It is widely acknowledged that an imbalanced biomechanical environment can have significant effects on myocardial pathology, leading to adverse remodelling of cardiac function if it persists. Accurate stress prediction essentially depends on the strain energy function which should have competent descriptive and predictive capabilities. Previous studies have focused on myofibre dispersion, but not on fibres along other directions. In this study, we will investigate how fibre dispersion affects myocardial biomechanical behaviours by taking into account both the myofibre dispersion and the sheet fibre dispersion, with a focus on the sheet fibre dispersion. Fibre dispersion is incorporated into a widely-used myocardial strain energy function using the discrete fibre bundle approach. We first study how different dispersion affects the descriptive capability of the strain energy function when fitting to ex vivo experimental data, and then the predictive capability in a human left ventricle during diastole. Our results show that the chosen strain energy function can achieve the best goodness-of-fit to the experimental data by including both fibre dispersion. Furthermore, noticeable differences in stress can be found in the LV model. Our results may suggest that it is necessary to include both dispersion for myofibres and the sheet fibres for the improved descriptive capability to the ex vivo experimental data and potentially more accurate stress prediction in cardiac mechanics.
Citation: Debao Guan, Yuqian Mei, Lijian Xu, Li Cai, Xiaoyu Luo, Hao Gao. Effects of dispersed fibres in myocardial mechanics, Part I: passive response[J]. Mathematical Biosciences and Engineering, 2022, 19(4): 3972-3993. doi: 10.3934/mbe.2022183
Related Papers:
[1]
Cameron Browne .
Immune response in virus model structured by cell infection-age. Mathematical Biosciences and Engineering, 2016, 13(5): 887-909.
doi: 10.3934/mbe.2016022
[2]
A. M. Elaiw, N. H. AlShamrani, A. D. Hobiny .
Stability of an adaptive immunity delayed HIV infection model with active and silent cell-to-cell spread. Mathematical Biosciences and Engineering, 2020, 17(6): 6401-6458.
doi: 10.3934/mbe.2020337
[3]
Tinevimbo Shiri, Winston Garira, Senelani D. Musekwa .
A two-strain HIV-1 mathematical model to assess the effects of chemotherapy on disease parameters. Mathematical Biosciences and Engineering, 2005, 2(4): 811-832.
doi: 10.3934/mbe.2005.2.811
[4]
Yu Wu, Xiaopeng Zhao, Mingjun Zhang .
Dynamics of stochastic mutation to immunodominance. Mathematical Biosciences and Engineering, 2012, 9(4): 937-952.
doi: 10.3934/mbe.2012.9.937
[5]
Ting Guo, Zhipeng Qiu .
The effects of CTL immune response on HIV infection model with potent therapy, latently infected cells and cell-to-cell viral transmission. Mathematical Biosciences and Engineering, 2019, 16(6): 6822-6841.
doi: 10.3934/mbe.2019341
[6]
Abdessamad Tridane, Yang Kuang .
Modeling the interaction of cytotoxic T lymphocytes and influenza
virus infected epithelial cells. Mathematical Biosciences and Engineering, 2010, 7(1): 171-185.
doi: 10.3934/mbe.2010.7.171
[7]
A. M. Elaiw, A. S. Shflot, A. D. Hobiny .
Stability analysis of general delayed HTLV-I dynamics model with mitosis and CTL immunity. Mathematical Biosciences and Engineering, 2022, 19(12): 12693-12729.
doi: 10.3934/mbe.2022593
[8]
Khalid Hattaf, Noura Yousfi .
Dynamics of SARS-CoV-2 infection model with two modes of transmission and immune response. Mathematical Biosciences and Engineering, 2020, 17(5): 5326-5340.
doi: 10.3934/mbe.2020288
[9]
Jie Lou, Tommaso Ruggeri, Claudio Tebaldi .
Modeling Cancer in HIV-1 Infected Individuals: Equilibria, Cycles and Chaotic Behavior. Mathematical Biosciences and Engineering, 2006, 3(2): 313-324.
doi: 10.3934/mbe.2006.3.313
[10]
Xuejuan Lu, Lulu Hui, Shengqiang Liu, Jia Li .
A mathematical model of HTLV-I infection with two time delays. Mathematical Biosciences and Engineering, 2015, 12(3): 431-449.
doi: 10.3934/mbe.2015.12.431
Abstract
It is widely acknowledged that an imbalanced biomechanical environment can have significant effects on myocardial pathology, leading to adverse remodelling of cardiac function if it persists. Accurate stress prediction essentially depends on the strain energy function which should have competent descriptive and predictive capabilities. Previous studies have focused on myofibre dispersion, but not on fibres along other directions. In this study, we will investigate how fibre dispersion affects myocardial biomechanical behaviours by taking into account both the myofibre dispersion and the sheet fibre dispersion, with a focus on the sheet fibre dispersion. Fibre dispersion is incorporated into a widely-used myocardial strain energy function using the discrete fibre bundle approach. We first study how different dispersion affects the descriptive capability of the strain energy function when fitting to ex vivo experimental data, and then the predictive capability in a human left ventricle during diastole. Our results show that the chosen strain energy function can achieve the best goodness-of-fit to the experimental data by including both fibre dispersion. Furthermore, noticeable differences in stress can be found in the LV model. Our results may suggest that it is necessary to include both dispersion for myofibres and the sheet fibres for the improved descriptive capability to the ex vivo experimental data and potentially more accurate stress prediction in cardiac mechanics.
1.
Introduction
T cells are one of the components of the adaptive immune system. They are a subset of lymphocytes that mature in the thymus and are responsible for searching for pathogens. The presence of a pathogen will result in the presentation of its characteristic peptides by the antigen presenting cells (APC) [1]. When T cells of a given clonotype find their specific peptide, they will become activated and they will start secreting cytokines, namely interleukine 2 (IL-2), that signal to other components of the immune system and promote proliferation [2]. However, it may happen that a clonotype of T cells erroneously target self-antigens, thereby promoting an auto-immune response against its host. Regulatory T cells (Tregs) are a subset of T cells that have immune suppressive function – they are able to inhibit cytokine secretion by T cells [3,4,5,6,7,8]. At a genetic level, Tregs express Foxp3, a master regulator of the Treg phenotype inducing CD25, CTLA-4 and GITR expression, all correlating with a suppressive phenotype [9]. Auto-immune diseases may appear when Tregs are misregulated, e.g. IPEX [9].
Regarding laboratory data obtained in vivo or in vitro, there is a large literature on the subject, see, for instance, the reviews by Zhu et al.[2,10]. However, in vivo data on time evolution of the concentration of T cells with several measurements over a period of days or months are relatively scarce. One example is the data on immune response by T cells from mice infected withlymphocytic choriomeningitis virus (LCMV) by Homann et al.[11]. In vivo time series can be information is important to further understand immune response dynamics and to validate existing deterministic or stochastic immune response models.
Modelling of immune responses by T cells can be made using different mathematical tools. See Callard et al.[12] and Lythe et al.[13] for reviews. De Boer et al.[14], Burroughs et al.[3], Pinto et al.[15], Blyuss et al.[16] and Khailaie et al.[17] studied systems of ordinary differential equations (ODE), while León et al.[18] used a hypergeometric distribution in a discrete model, and Bianca and Brézin [19] used thermostatted kinetic theory methods. In [20] the authors studied bifurcations and employ stochastic modelling with the goal of addressing individual variability as well as variability in the immune system as a whole by means of the Doob-Gillespie algorithm.
In this paper we use the ODE model with two clonotypes in Pinto et al.[15], based on the model from Burroughs et al.[3]. The model with two clonotypes was previously used to qualitatively describe the appearance of autoimmune diseases after an infection due to the bystanter proliferation of the clonotype of autoimmune T cells [21] and the suppression of autoimmune responses as a consequence of decrease of the clonotype of autoimmune T cells during an infection with sufficiently long duration [22]. Previous research on the model with one clonotype [3,23] showed that T cells' growth due to cytokine has a quorum population threshold as suggested in de Boer et al.[14]. For the reported values of the parameters, below a given threshold of antigen stimulation of T cells b<bL, only one stable equilibrium is found - a controlled state with a low concentration of T cells. Above a higher threshold of antigenic stimulation of T cells b>bH, only one stable equilibrium is present - an immune response state, characterized by a high concentration of T cells. Between the two antigenic stimulation thresholds, for b in [bL,bH], there is a bi-stability region, where both an immune response state and a controlled state are present. The thresholds bL and bH are saddle-node bifurcation points that bound the hysteresis, see e.g. [3,22]. Moreover, some parameters can unfold the hysteresis at a cusp bifurcation point [23,24,25]. When a linear tuning between the antigenic stimuli of T cells and Tregs is considered, we observe that the thresholds of antigenic stimulation of T cells can change substantially with the slope parameter and we can also observe a transcritical bifurcation and other saddle-node bifurcations [26,27]. See also the reviews by Burroughs et al.[28] and by Pinto et al.[15].
In this paper, we aim to calibrate the ODE model by obtaining a quantitative fit to data, thereby complementing the qualitative description of the properties of the model from previous publications. Thus, we fit the time dynamics of the ODE model to the data from Homann et al.[11] regarding the CD4+ T cell numbers responding to LCMV epitopes gp61 or NP309. Previously, de Boer et al.[29] fitted the data on the same epitopes considering two or three separate phases. The first phase of their model is an exponential growth phase, when the immune system is simulated. It is followed by one or two exponential decay phases, when the stimulation is small. De Boer et al. used the same methodology to study data from other epitopes in a different publication [30]. We use a different model from de Boer et al. since we consider simultaneously two clonotypes of T cells, Tregs and IL-2; and our ODE system includes non-linear terms. The novelty of our approach is to fit a unique ODE model to the growth and decay phases, where only the antigenic stimulation of T cells may have different values between the immune response phase and the contraction phase. This article is structured as follows. We start by describing the adopted immune response model to fit the data in Section 2. Next, in Section 3 we fit the model to the data. Finally, in Section 4 we present a discussion and an analysis of the results obtained in the fits. The parameters of the model are presented in the Appendix 7.
2.
Immune response model
We use the model from Section 3 in Pinto et al.[15] for CD4+ T cells and regulatory T cells described by a set of ordinary differential equations. We consider one population of Tregs (R,R∗) and two populations of conventional T cells (Tb,Tb∗) and (Tc,Tc∗), responding to different antigens. Both populations of T cells and Tregs require specific antigenic stimulation for activation. Levels of antigenic stimulation are denoted by a for Tregs and by b and c for conventional T cells. T cells are activated by their specific antigen, from a non-secreting state, denoted Tb or Tc, to a IL-2 secreting state Tb∗ or Tc∗. In this paper, the values of b and c can change along time as a consequence of pathogen presence. Similarly, Tregs are activated by self antigens at a constant level a from an inactive state R, to an active state R∗. Tregs do not secrete IL-2. Conventional T cells acquire proliferative capacity in the presence of IL-2. Some authors suggest that Tregs proliferate in vivo with IL-2 [4,5,8], while the opposite may happen in vitro[7]. We will consider that Tregs also proliferate in the presence of IL-2, although less efficiently than normal T cells [3,15]. Secreting T cells T∗ can revert to the non-secreting state T at a rate k or due to the non-specific suppression by active Tregs R∗, at a rate proportional to γ. Nevertheless, non-secreting T cells are still able to proliferate in the presence of IL-2 [6,9]. Moreover, we consider two types of death terms: linear, the d terms, and quadratic, the β terms. The latter acts as a saturation mechanism and can be related to Fas-FasL death by apoptosis [31,32]. Finally, we include an influx of (auto) immune T cells into the tissue (Tinb and Tinc) and Tregs (Rin), which can represent previously activated T cells from circulation or naïve T cell input from the thymus. This model does not include memory cells explicitly, hence it may only be suited to simulate time evolutions over a few weeks.
The model consists of a set of seven ordinary differential equations. We have a compartment for each T cell population (inactive Tregs R, active Tregs R∗, non-secreting T cells Tb and Tc, and secreting T cells Tb∗ and Tc∗) and interleukine 2 density I.
with N=R+R∗+Tb+Tb∗+Tc+Tc∗. The parameters for this model and the ranges considered in the optimization procedure for the fits are given in the appendix (7) in Tables 1 and 2. We will consider the death rates to be equal dT=dR and dR∗=dT∗. We also assume equal relaxation rates k=ˆk, and different. Furthermore, we allow for different inflows of T cells Tinb and Tinc and Tregs Rin.
3.
Simulations and fits
In this section we do the fitting of the immune response model described in the previous section to the data from Homann et al.[11]. The data consist of two time series that characterize the immune response by CD4+ T cells to lymphocytic choriomeningitis virus - LCMV in mice obtained from laboratory experiments. Each time series corresponds to the total number of T cells that respond to one of the LCMV epitopes studied here: gp61 and NP309.
Most of the parameters of the model were obtained from the literature [3,15,23,33,34,35,36,37] and we assumed those presented in Table 1 to be fixed a priori. In Table 2 we list the parameters that were part of the optimization procedure to obtain the fit. These parameters are the T cell maximum growth rate ρ/α, the capacity of T cells Tcap, the level of homeostatic Tb cells that are triggered by the gp61 epitope, and the ratio between the homeostatic levels of T cells that are triggered by the NP309 epitope with those triggered by the gp61 epitope.
The pathogen was considered to be present within a limited period of time. To model this we used step functions for the antigen simulation parameters b and c, with specific parameters for each clonotype of T cells responding to each epitope: bmax for the gp61 epitope and cmax for the NP309 epitope. More precisely, we assumed b and c to be functions of time having a certain high antigen stimulation intensity b=bmax and c=cmax within a given time interval t∈[tini,tend], and very small values, b=10−5 and c=10−5 outside that interval.
The values of the six parameters of the fit were obtained from the optimization procedure. The software used for the simulations of the model was GNU Octave - version 4.2.2 and the routine that was used for the numerical integration of the ODE system was lsode. The initial condition of the system (at time t=0) for each set of parameters was obtained by simulating the system for a long time with small values of antigenic stimulations of T cells, b=c=10−5. This allowed us to obtain an initial condition close to the homeostatic controlled state, with a low total number of T cells.
The optimization procedure aimed to minimize the sum of squares of the residuals in a log scale, i.e. to minimize
where the values xb(ti) and xc(ti) are the observations of total T cell numbers from Homann et al.[11] for the gp61 and NP309 epitopes respectively, ti represents the time of each observation and Tb(ti)+Tb∗(ti) and Tc(ti)+Tc∗(ti) the total T cell numbers as obtained from the ODE fit at time ti. The total residual sum of squares has two components, one corresponding to the gp61 epitope time series res2b and another corresponding to the NP309 epitope time series res2c.
The optimum parameter estimates were obtained by using Octave's routine fminsearch implementing the Nelder-Mead algorithm/downhill simplex method. We started by running the algorithm on 500 random initial points in the six-dimensional parameter space. Due to numerical accuracy issues, the algorithm could stop at local sup-optimal minimum points. Hence, we had the necessity of running the algorithm repeatedly from distinct initial points. The best fit estimates of the parameters to the data, as well as 95% confidence intervals are presented in Table 2. The confidence intervals were obtained by deviating each parameter from the optimal value until the residual sum of squares was significantly larger - according to the F statistic - than the minimum residual [38]. The fit is presented in Figure 1 together with the residual mean square MNSQ=res2/df, where the number of degrees of freedom df is the difference between the number of data points and the number of free parameters. For epitope gp61 we have that res2b=0.565 and for epitope NP309 we have res2c=0.995. Hence, for 14+8−6=16 degrees of freedom, we have that MNSQ=0.0975.
Figure 1.
Dynamics for the best fit of the ODE model to the laboratory measurements of CD4+ T cells numbers from Homman et al.[11]. MNSQ=0.0975. Circles: gp61 epitope data. Crosses: NP309 epitope data. The vertical axis is in logarithmic scale. The vertical dashed lines mark the time of the beginning tini=3.31 and the end tend=8.74; of the immune response phase. The parameters of the fit are shown in Tables 1 and 2. Top: Total Tb+Tb∗ (green line); non-secreting Tb (blue dots); secreting Tb∗ (red dashes). Middle: Total Tc+Tc∗ (green line); non-secreting Tc (blue dots); secreting Tc∗ (red dashes). Bottom: Total R+R∗ (green line); inactive R (blue dots); active R∗ (red dashes).
We observe in Figure 1 that the total number of T cells responding to either epitope follows the data from Homann et al.[11]. Before tini the system has not detected the presence of a pathogen, thus being in homeostasis. During this period, the numbers of non-secreting T cells Tb and Tc are very close to their respective total. Between tini and tend, the numbers of secreting T cells Tb∗ and Tc∗ initially dominate, and that after day 7, approximately, the numbers of the non-secreting T cells become larger than the numbers of secreting T cells. This happens because the number of Tregs (R and R∗) increases due to the increased presence of IL-2. When the number of active Tregs R∗ is sufficiently high, the suppressive action of Tregs is stronger than the antigenic stimulation of T cells. Hence, as the number of active Tregs increases, the secretion of IL-2 cytokines becomes increasingly inhibited, thus reducing the proportion of secreting T cells, even if the total number of T cells is still increasing. Moreover, the overall larger numbers of cells result in an increased Fas-FasL (quadratic) death rate, thereby decreasing the balance between growth and death. After tend, growth halts, since the antigenic stimulation of T cells nearly vanishes. Hence, almost all the secreting T cells relax to their non-secreting state leading to a very low secretion rate of IL-2 and to its quick decay. Therefore, in this phase, the dynamics is governed by the death terms. Initially, the Fas-FasL quadratic death term dominates the contraction phase of all cells. As their number is decreasing, the rate of change is regulated by the balance between the linear death and the inflow of cells into the system.
Besides the fit for the first 4 weeks we also obtained the fit for a time span of 3 weeks (13 data points for epitope gp61 and 7 data points for epitope NP309) and 8 weeks (13 data points for epitope gp61 and 7 data points for epitope NP309). In both cases, the fitted values of the parameters belong to the confidence intervals in Table 2: for the 3 weeks fit we obtained: MNSQ=0.0745, ρ/α=2.06, Tcap=1.22×107, Thomb=2480, Thomc/Thomb=0.279, tini=3.23, and tend=9.36; while for the 8 weeks fit we got: MNSQ=0.1207, ρ/α=2.15, Tcap=2.54×107, Thomb=1732, Thomc/Thomb=0.363, tini=3.22, and tend=8.50.
4.
Discussion
Our simulations can fit the data both in the period of high antigenic stimulation around the first week, and afterwards until day 28, during the contraction phase. De Boer et al.[29] analysed previously the data from Homann et al.[11]. Regarding epitope gp61, their model with a biphasic contraction model had a residual mean square MNSQgp61=0.06, while the one with a single contraction phase had MNSQgp61=0.17. Regarding the NP309 epitope, their fits had similar residual mean squares: both the biphasic contraction model MNSQNP309=0.10 and the single contraction model MNSQNP309=0.12. We note that our model is different since we are fitting simultaneously the data from gp61 and NP309 epitopes using 7 non-linear ordinary differential equations, we use two death terms and we change the antigenic stimulations of T cells b and c in time, instead of having a distinct linear differential equation for each phase, as in de Boer et al.[29]. Using step functions for the antigenic stimulation of T cells can be interpreted biologically as the pathogen being either present or not at each time period, which is a simplification of reality in itself. However, when considering the transition from the immune response phase to the contraction phase, using step functions for parameters b and c should not differ substantially from using a distinct equation in each phase, since the transitions between the non-secreting and the secreting states occur in time scales of hours, while data are daily. However, the non-linearities of the model we used, allowed for a smooth fit during the contraction phase.
The assumption of the initial condition being the controlled state is consistent with the assumption that the mice T cells were in a homeostatic state, which is also consistent with the laboratory methodology described in [11] that reports the mice being infected at 6-8 weeks of age. In our fit, a larger homeostatic population of T cells responding to epitope gp61, Thomb>Thomc, and similar growth rates for both populations of T cells, yields a similar outcome at day 8-9 to the model presented in de Boer et al.[29] in which at around day 3, the concentration of T cells is similar for both clonotypes, but the growth rate of the T cells responding to epitope gp61 is larger. We note that in our model, the growth rate of T cells is not constant since it is a function of the concentration of the IL-2 cytokine, which in turn depends on the concentration of secreting T cells. Thus, we were not surprised by finding that the fitted value of the maximum growth rate was higher than the average growth rate reported by de Boer et al.[29].
The chosen value of the death rate of T cells is lower than the death rate of the first decay phase, the apoptosis phase, in de Boer et al.[29], being closer to the death rate of memory T cells from the second phase. Other small values of the death rate resulted in fits with similar mean square residuals. Thus, we opted to keep the death rate parameter fixed. We note that our model did not consider the differentiation to memory cells which led us to consider only the first 28 days of the data in Homann et al.[11]. This may represent a limitation that might be addressed in the future by considering a different model with a specific compartment for memory T cells. The asymmetry in the death rates between secreting and non-secreting T cells, from Pinto et al.[15], was present since we considered more than one clonotype of T cells (see also e.g. Burroughs et al.[21] and Oliveira et al.[22]). As for the fitted values for time spans of 3 weeks and time spans of 8 weeks, we observed that they were within the confidence intervals for the fits for 4 weeks, hence exhibiting no sensitivity at least with respect to small changes in the considered time span. Regarding the Fas-FasL death term, it allows a continuous deceleration of the death rate during the contraction period. We can observe this effect more visibly in Figure 4 (A) of Pinto et al.[15]. This deceleration can be explained by the different relative importance of the death terms along time. After the end of the immune activation phase, while the population size of T cells is still large, the Fas-FasL death rate, proportional to the square of the concentration of T cells, is also high. Hence, we observe an apoptosis phase, characterized by a fast decrease in T cell numbers. As the T cells die, this quadratic term looses relevance and the death rate approaches the linear death term. This smooth behaviour is qualitatively similar to the biphasic contraction phase described in de Boer et al.[29].
Another limitation is to consider that the antigenic stimulations of T cells b and c, and thus pathogen activation intensity, are step functions of time instead of smooth functions, which might be more realistic. However, this would add more parameters to the model. We recall that the parameters that were kept constant had their values obtained from the literature, and we opted to have the same values for both epitopes due to the small number of data points available. Biologically, having the same values may be interpreted as meaning that the respective characteristics of the T cells are similar for the two epitopes.
We used a model that has processes occurring in different time scales, from minutes, to hours to days. The laboratory restrictions, that led to a small number of data points for a continuous time fit, allow for other models to obtain good fits to the data from Homann et al.[11]. Furthermore, to better understand the time dynamics of Tregs in immune responses, we would require that their numbers were also measured. Hence, the homeostatic number of Tregs may be different. Other parameter values had fits with similar residual mean squares, provided the secretion inhibition was adjusted accordingly (data not shown). Nevertheless, in our fit, the presence of Tregs is felt on the decrease of the growth rate, observed at the end of the immune activation phase.
Our model has two main strengths. Firstly, it fits simultaneously the data for the gp61 and NP309 epitopes using the same equations in all phases, the immune activation phase and the contraction phase, and without the need to consider biphasic regimes. Secondly, with respect to the gp61 and NP309 epitopes, only the T cells homeostatic values were different.
For the future, it would be interesting to study human data on the concentration of T cells and Tregs, over large periods of time, ideally measured with a small time resolution. With longitudinal data, we could observe the natural occurrence of infections and the following recovery. This data would allow the researchers to use mathematical tools to create and to validate dynamical models with the inclusion of memory cells, that are key to develop long term immunity to some infections. Moreover, dynamics of auto-immune diseases could also be studied and modelled, and in the process we could estimate parameters and their confidence intervals. However, ethical aspects of such studies should be carefully evaluated, weighting the benefits and the costs.
In conclusion, the model we used allows to be calibrated to obtain a good fit to the data from Homann et al.[11] keeping most of the parameters fixed at the values obtained from the literature, and optimizing a set of 6 parameters that are related with the beginning and ending times of the immune response, the maximum growth rate of T cells, the T cells capacity, and the homeostatic numbers of T cells responding to gp61 and NP309 epitopes. Therefore, besides the qualitative features of the model described in other publications, we showed that this model is suited to describe quantitatively immune responses by CD4+ T cells.
Acknowledgements
The authors would like to thank the financial support by the ERDF – European Regional Development Fund through the Operational Programme for Competitiveness and Internationalization – COMPETE 2020 Programme within project "POCI-01-0145-FEDER-006961" and by national Funds through the FCT – Fundação para a Ciência e a Tecnologia (Portuguese Foundation for Science and Technology) as part of project UID/ EEA/ 50014/ 2013, and within project "Dynamics, Optimization and Modelling" with reference PTDC/ MAT-NAN/ 6890/ 2014 as part of project "Modelling, Dynamics and Games" with reference PTDC/MAT-APL/31753/2017 and project "NanoSTIMA – Macro-to-Nano Human Sensing: Towards Integrated Multimodal Health Monitoring and Analytics" with reference NORTE-01-0145-FEDER-000016, financed by the North Portugal Regional Operational Programme (NORTE 2020), under the PORTUGAL 2020 Partnership Agreement. Atefeh Afsar would like to thank the financial support of FCT - Fundação para a Ciência e a Tecnologia - through a PhD. grant of the MAP-PDMA program with reference PD/BD/142886/2018.
Conflict of interest
The authors declare no conflict of interest.
A.
Appendix
The parameters of our model and their default values are presented in the following two tables.
Table 1.
Fixed parameter values for the model of T cells and Tregs from [3,15,23,33,34,35,36,37]. The formula for the maximum cytokine concentration is obtained from [3,23].
Fixed Parameters
Symbol
Value
T cells Tb,Tb∗,Tc,Tc∗
Death rate of inactive T cells (day−1)
dT
10−3
Death rate ratio of active : inactive T cells
dT∗/dT
0.1
Secretion reversion2 (day−1)
k
2.4
Maximum antigen stimulation level of Tb cells
bmax
100
Maximum antigen stimulation level of Tc cells
cmax
100
Tregs R, R∗
Growth rate ratio Tregs : T cells
ϵ
0.6
Relaxation rate (day−1)
ˆk
2.4
Death rate ratio of inactive Tregs : inactive T cells
1 Maximum T cell density for severe infections (based on LCMV).
2 This is in the absence of Tregs.
3 Tregs input level is given by Rin=Rhom(dR−ˆka(dR−dR∗)dR∗+ˆk(1+a)).
4 This is taken as 20 times the receptor affinity (10pM).
Table 2.
Fitted parameter values and 95% confidence intervals for the model of T cells and Tregs. MNSQ=0.0975. The formulas for the T cell maximum growth rate and the capacity of T cells were obtained from [3,23].
3Tc cell input level is given by Tinc=TinbThomc/Thomb.
References
[1]
M. R. Zile, C. F. Baicu, W. H. Gaasch, Diastolic heart failure—abnormalities in active relaxation and passive stiffness of the left ventricle, N. Engl. J. Med., 350 (2004), 1953–1959. https://doi.org/10.1056/NEJMoa032566 doi: 10.1056/NEJMoa032566
[2]
K. Mangion, H. Gao, D. Husmeier, X. Luo, C. Berry, Advances in computational modelling for personalised medicine after myocardial infarction, Heart, 104 (2018), 550–557. http://dx.doi.org/10.1136/heartjnl-2017-311449 doi: 10.1136/heartjnl-2017-311449
[3]
G. A. Holzapfel, R. W. Ogden, Constitutive modelling of passive myocardium: a structurally based framework for material characterization, Philos. Trans. R. Soc., A, 367 (2009), 3445–3475. https://doi.org/10.1098/rsta.2009.0091 doi: 10.1098/rsta.2009.0091
[4]
H. Gao, H. Wang, C. Berry, X. Luo, B. E. Griffith, Quasi-static image-based immersed boundary-finite element model of left ventricle under diastolic loading, Int. J. Numer. Method. Biomed. Eng., 30 (2014), 1199–1222. https://doi.org/10.1002/cnm.2652 doi: 10.1002/cnm.2652
[5]
H. Gao, A. Aderhold, K. Mangion, X. Luo, D. Husmeier, C. Berry, Changes and classification in myocardial contractile function in the left ventricle following acute myocardial infarction, J. R. Soc. Interface, 14 (2017), 20170203. https://doi.org/10.1098/rsif.2017.0203 doi: 10.1098/rsif.2017.0203
[6]
B. Baillargeon, N. Rebelo, D. D. Fox, R. L. Taylor, E. Kuhl, The living heart project: a robust and integrative simulator for human heart function, Eur. J. Mech. A/Solids, 48 (2014), 38–47. https://doi.org/10.1016/j.euromechsol.2014.04.001 doi: 10.1016/j.euromechsol.2014.04.001
[7]
K. L. Sack, E. Aliotta, D. B. Ennis, J. S. Choy, G. S. Kassab, J. M. Guccione, et al., Construction and validation of subject-specific biventricular finite-element models of healthy and failing swine hearts from high-resolution dt-mri, Front. Physiol., 9 (2018). https://doi.org/10.3389/fphys.2018.00539
[8]
H. Gao, K. Mangion, D. Carrick, D. Husmeier, X. Luo, C. Berry, Estimating prognosis in patients with acute myocardial infarction using personalized computational heart models, Sci. Rep., 7 (2017), 13527. https://doi.org/10.1038/s41598-017-13635-2 doi: 10.1038/s41598-017-13635-2
[9]
S. I. H. Richardson, H. Gao, J. Cox, R. Janiczek, B. E. Griffith, C. Berry, et al., A poroelastic immersed finite element framework for modelling cardiac perfusion and fluid–structure interaction, Int. J. Numer. Method. Biomed. Eng., 37 (2021), e3446. https://doi.org/10.1002/cnm.3446 doi: 10.1002/cnm.3446
[10]
H. M. Wang, H. Gao, X. Y. Luo, C. Berry, B. E. Griffith, R. W. Ogden, et al., Structure-based finite strain modelling of the human left ventricle in diastole, Int. J. Numer. Method. Biomed. Eng., 29 (2013), 83–103. https://doi.org/10.1002/cnm.2497 doi: 10.1002/cnm.2497
[11]
D. Guan, X. Luo, H. Gao, Constitutive modelling of soft biological tissue from ex vivo to in vivo: myocardium as an example, in International Conference by Center for Mathematical Modeling and Data Science, Osaka University, Springer, (2020), 3–14. https://doi.org/10.1007/978-981-16-4866-3_1
[12]
H. Gao, W. G. Li, L. Cai, C. Berry, X. Y. Luo, Parameter estimation in a holzapfel–ogden law for healthy myocardium, J. Eng. Math., 95 (2015), 231–248. https://doi.org/10.1007/s10665-014-9740-3 doi: 10.1007/s10665-014-9740-3
[13]
U. Noè, A. Lazarus, H. Gao, V. Davies, B. Macdonald, K. Mangion, et al., Gaussian process emulation to accelerate parameter estimation in a mechanical model of the left ventricle: a critical step towards clinical end-user relevance, J. R. Soc. Interface, 16 (2019), 20190114. https://doi.org/10.1098/rsif.2019.0114 doi: 10.1098/rsif.2019.0114
[14]
V. Davies, U. Noè, A. Lazarus, H. Gao, B. Macdonald, C. Berry, et al., Fast parameter inference in a biomechanical model of the left ventricle by using statistical emulation, J. R. Stat. Soc. Ser. C Appl. Stat., 68 (2019), 1555–1576. https://doi.org/10.1111/rssc.12374 doi: 10.1111/rssc.12374
[15]
G. Sommer, A. J. Schriefl, M. Andrä, M. Sacherer, C. Viertler, H. Wolinski, et al., Biomechanical properties and microstructure of human ventricular myocardium, Acta Biomater., 24 (2015), 172–192. https://doi.org/10.1016/j.actbio.2015.06.031 doi: 10.1016/j.actbio.2015.06.031
[16]
F. Ahmad, S. Soe, N. White, R. Johnston, I. Khan, J. Liao, et al., Region-specific microstructure in the neonatal ventricles of a porcine model, Ann. Biomed. Eng., 46 (2018), 2162–2176. https://doi.org/10.1007/s10439-018-2089-4 doi: 10.1007/s10439-018-2089-4
[17]
Y. Lanir, Multi-scale structural modeling of soft tissues mechanics and mechanobiology, J. Elast., 129 (2017), 7–48. https://doi.org/10.1007/s10659-016-9607-0 doi: 10.1007/s10659-016-9607-0
[18]
T. C. Gasser, R. W. Ogden, G. A. Holzapfel, Hyperelastic modelling of arterial layers with distributed collagen fibre orientations, J. R. Soc. Interface, 3 (2006), 15–35. https://doi.org/10.1098/rsif.2005.0073 doi: 10.1098/rsif.2005.0073
[19]
T. S. E. Eriksson, A. J. Prassl, G. Plank, G. A. Holzapfel, Modeling the dispersion in electromechanically coupled myocardium, Int. J. Numer. Method. Biomed. Eng., 29 (2013), 1267–1284. https://doi.org/10.1002/cnm.2575 doi: 10.1002/cnm.2575
[20]
G. A. Holzapfel, R. W. Ogden, S. Sherifova, On fibre dispersion modelling of soft biological tissues: a review, Proc. Math. Phys. Eng. Sci., 475 (2019). https://doi.org/10.1098/rspa.2018.0736
[21]
A. V. Melnik, X. Luo, R. W. Ogden, A generalised structure tensor model for the mixed invariant i8, Int. J. Non-Linear Mech., 107 (2018), 137–148. https://doi.org/10.1016/j.ijnonlinmec.2018.08.018 doi: 10.1016/j.ijnonlinmec.2018.08.018
[22]
A. Pandolfi, A. Gizzi, M. Vasta, Coupled electro-mechanical models of fiber-distributed active tissues, J. Biomech., 49 (2016), 2436–2444. https://doi.org/10.1016/j.jbiomech.2016.01.038 doi: 10.1016/j.jbiomech.2016.01.038
[23]
A. Gizzi, A. Pandolfi, M. Vasta, Statistical characterization of the anisotropic strain energy in soft materials with distributed fibers, Mech. Mater., 92 (2016), 119–138. https://doi.org/10.1016/j.mechmat.2015.09.008 doi: 10.1016/j.mechmat.2015.09.008
[24]
D. Guan, X. Zhuan, W. Holmes, X. Luo, H. Gao, Modelling of fibre dispersion and its effects on cardiac mechanics from diastole to systole, J. Eng. Math., 128 (2021), 1–24. https://doi.org/10.1007/s10665-021-10102-w doi: 10.1007/s10665-021-10102-w
[25]
G. A. Holzapfel, R. W. Ogden, On fiber dispersion models: exclusion of compressed fibers and spurious model comparisons, J. Elast., 129 (2017), 49–68. https://doi.org/10.1007/s10659-016-9605-2 doi: 10.1007/s10659-016-9605-2
[26]
K. Li, R. W. Ogden, G. A. Holzapfel, A discrete fibre dispersion method for excluding fibres under compression in the modelling of fibrous tissues, J. R. Soc. Interface, 15 (2018), 20170766. https://doi.org/10.1098/rsif.2017.0766 doi: 10.1098/rsif.2017.0766
[27]
D. Guan, J. Yao, X. Luo, H. Gao, Effect of myofibre architecture on ventricular pump function by using a neonatal porcine heart model: from dt-mri to rule-based methods, R. Soc. Open Sci., 7 (2020), 191655. https://doi.org/10.1098/rsos.191655 doi: 10.1098/rsos.191655
[28]
K. Li, R. W. Ogden, G. A. Holzapfel, Computational method for excluding fibers under compression in modeling soft fibrous solids, Eur. J. Mech. A/Solids, 57 (2016), 178–193. https://doi.org/10.1016/j.euromechsol.2015.11.003 doi: 10.1016/j.euromechsol.2015.11.003
[29]
M. Vasta, A. Gizzi, A. Pandolfi, On three-and two-dimensional fiber distributed models of biological tissues, Probab. Eng. Mech., 37 (2014), 170–179. https://doi.org/10.1016/j.probengmech.2014.05.003 doi: 10.1016/j.probengmech.2014.05.003
[30]
G. A. Holzapfel, J. A. Niestrawska, R. W. Ogden, A. J. Reinisch, A. J. Schriefl, Modelling non-symmetric collagen fibre dispersion in arterial walls, J. R. Soc. Interface, 12 (2015), 20150188. https://doi.org/10.1098/rsif.2015.0188 doi: 10.1098/rsif.2015.0188
[31]
F. Ahmad, J. Liao, S. Soe, M. D. Jones, J. Miller, P. Berthelson, et al., Biomechanical properties and microstructure of neonatal porcine ventricles, J. Mech. Behav. Biomed. Mater., 88 (2018), 18–28. https://doi.org/10.1016/j.jmbbm.2018.07.038 doi: 10.1016/j.jmbbm.2018.07.038
[32]
G. A. Holzapfel, T. C. Gasser, R. W. Ogden, Comparison of a multi-layer structural model for arterial walls with a fung-type model, and issues of material stability, J. Biomech. Eng., 126 (2004), 264–275. https://doi.org/10.1115/1.1695572 doi: 10.1115/1.1695572
[33]
D. Guan, F. Ahmad, P. Theobald, S. Soe, X. Luo, H. Gao, On the aic-based model reduction for the general holzapfel–ogden myocardial constitutive law, Biomech. Model. Mechanobiol., 18 (2019), 1213–1232. https://doi.org/10.1007/s10237-019-01140-6 doi: 10.1007/s10237-019-01140-6
[34]
S. Klotz, I. Hay, M. L. Dickstein, G. Yi, J. Wang, M. S. Maurer, et al., Single-beat estimation of end-diastolic pressure-volume relationship: a novel method with potential for noninvasive application, Am. J. Physiol. Heart Circ. Physiol., 291 (2006), H403–H412. https://doi.org/10.1152/ajpheart.01240.2005 doi: 10.1152/ajpheart.01240.2005
[35]
A. V. Melnik, H. B. Da Rocha, A. Goriely, On the modeling of fiber dispersion in fiber-reinforced elastic materials, Int. J. Non-Linear Mech., 75 (2015), 92–106. https://doi.org/10.1016/j.ijnonlinmec.2014.10.006 doi: 10.1016/j.ijnonlinmec.2014.10.006
[36]
D. H. Cortes, S. P. Lake, J. A. Kadlowec, L. J. Soslowsky, D. M. Elliott, Characterizing the mechanical contribution of fiber angular distribution in connective tissue: comparison of two modeling approaches, Biomech. Model. Mechanobiol., 9 (2010), 651–658. https://doi.org/10.1007/s10237-010-0194-x doi: 10.1007/s10237-010-0194-x
[37]
X. Zhuan, X. Luo, H. Gao, R. W. Ogden, Coupled agent-based and hyperelastic modelling of the left ventricle post-myocardial infarction, Int. J. Numer. Method. Biomed. Eng., 35 (2019), e3155. https://doi.org/10.1002/cnm.3155 doi: 10.1002/cnm.3155
[38]
G. A. Holzapfel, R. W. Ogden, An arterial constitutive model accounting for collagen content and cross-linking, J. Mech. Phys. Solids, 136 (2020), 103682. https://doi.org/10.1016/j.jmps.2019.103682 doi: 10.1016/j.jmps.2019.103682
[39]
J. Xi, P. Lamata, S. Niederer, S. Land, W. Shi, X. Zhuang, et al., The estimation of patient-specific cardiac diastolic functions from clinical measurements, Med. Image Anal., 17 (2013), 133–146. https://doi.org/10.1016/j.media.2012.08.001 doi: 10.1016/j.media.2012.08.001
[40]
M. Strocchi, M. A. F. Gsell, C. M. Augustin, O. Razeghi, C. H. Roney, A. J. Prassl, et al., Simulating ventricular systolic motion in a four-chamber heart model with spatially varying robin boundary conditions to model the effect of the pericardium, J. Biomech., 101 (2020), 109645. https://doi.org/10.1016/j.jbiomech.2020.109645 doi: 10.1016/j.jbiomech.2020.109645
This article has been cited by:
1.
Aliyu Yusuf, Isabel Figueiredo, Atefeh Afsar, Nigel Burroughs, Alberto Pinto, Bruno Oliveira,
The Effect of a Linear Tuning between the Antigenic Stimulations of CD4+ T Cells and CD4+ Tregs,
2020,
8,
2227-7390,
293,
10.3390/math8020293
2.
Atefeh Afsar, Filipe Martins, Bruno M. P. M. Oliveira, Alberto A. Pinto,
2021,
Chapter 1,
978-3-030-78162-0,
1,
10.1007/978-3-030-78163-7_1
3.
Atefeh Afsar, Maria del Carmen Chacon Castro, Adedamola Saidi Soladogun, Li Zhang,
Recent Development in the Understanding of Molecular and Cellular Mechanisms Underlying the Etiopathogenesis of Alzheimer’s Disease,
2023,
24,
1422-0067,
7258,
10.3390/ijms24087258
4.
Atefeh Afsar, Li Zhang,
Putative Molecular Mechanisms Underpinning the Inverse Roles of Mitochondrial Respiration and Heme Function in Lung Cancer and Alzheimer’s Disease,
2024,
13,
2079-7737,
185,
10.3390/biology13030185
5.
Atefeh Afsar, Min Chen, Zhenyu Xuan, Li Zhang,
A glance through the effects of CD4+ T cells, CD8+ T cells, and cytokines on Alzheimer's disease,
2023,
21,
20010370,
5662,
10.1016/j.csbj.2023.10.058
6.
Mingran Zhang, Junling Ma, Meili Li,
Original Antigenic Sin in CD4+ T Cells,
2025,
0019-2805,
10.1111/imm.13916
Table 1.
Fixed parameter values for the model of T cells and Tregs from [3,15,23,33,34,35,36,37]. The formula for the maximum cytokine concentration is obtained from [3,23].
Fixed Parameters
Symbol
Value
T cells Tb,Tb∗,Tc,Tc∗
Death rate of inactive T cells (day−1)
dT
10−3
Death rate ratio of active : inactive T cells
dT∗/dT
0.1
Secretion reversion2 (day−1)
k
2.4
Maximum antigen stimulation level of Tb cells
bmax
100
Maximum antigen stimulation level of Tc cells
cmax
100
Tregs R, R∗
Growth rate ratio Tregs : T cells
ϵ
0.6
Relaxation rate (day−1)
ˆk
2.4
Death rate ratio of inactive Tregs : inactive T cells
Table 2.
Fitted parameter values and 95% confidence intervals for the model of T cells and Tregs. MNSQ=0.0975. The formulas for the T cell maximum growth rate and the capacity of T cells were obtained from [3,23].
Death rate ratio of inactive Tregs : inactive T cells
dR/dT
1
Death rate relative ratio of Tregs : T cells
dR∗dR/dT∗dT
1
Tregs antigen stimulation level (day−1)
aˆk
1
Homeostatic Treg level3 (cells day−1)
Rhom
100
Secretion inhibition
γ
10/Rhom
Interleukin-2 (IL-2) Cytokine
Max. cytokine concentration4 (pM)
1/α
200
IL2 secretion rate (pM day−1)
σ
144
Cytokine decay rate (day−1)
σδ
36
Parameter
Symbol
Fitted value
95% CI1
T cell maximum growth rate (day−1)
ρ/α
2.06
[1.82,2.35]
Capacity of T cells1 (107 cells)
Tcap=ρ/(αβ)
1.76
[1.10,2.99]
Homeostatic Tb cells level2 (cells day−1)
Thomb
2409
[1637,3677]
Homeostatic levels ratio3Thomc : Thomb
Thomc/Thomb
0.327
[0.204,0.527]
Beginning of the immune activation phase (day)
tini
3.31
[3.00,3.59]
End of the immune activation phase (day)
tend
8.74
[7.94,9.91]
Figure 1. (a) Schematic of heart from Wikipedia (left). Microstructural arrangement of myofibres and sheet fibres in the myocardium (middle). The blue network describes sheet fibres that connect myofibres (red columns). Dispersed fibre field is drawn on the right. An unit vector fn (red) representing the myofibre direction defined by Θ and Φ with respect to the mean myofibre direction f0 in the (f0,s0,n0) fibre system. Similarly, the greed arrow is the dispersed sheet with θ and ϕ with respect to s0. (b) Illustration of the unit hemisphere domain centralised with the mean myofibre direction f0 (the red arrow). In the DFD method, it is divided into N discrete triangular elements with representative myofibre directions fn (green arrows) at the centroid of each triangular surface
Figure 2. Illustrations of ρ with η=4.5, η=3.9, η=0.0 and τ∈[−π/2,π/2]. When η=0.0, it is a uniform distribution
Figure 3. The left is the sketch of a left ventricle with inside myofibres (red lines) and a cubic sample cut from the ventricular wall. The right is a sketch of all six possible shear modes where f0, s0 and n0 denote the myofibre, sheet and sheet-normal direction, respectively. (ij) refers to shear in the j0 direction within the i0j0 plane, where i≠j∈{f,s,n}
Figure 4. Fibre dispersion distributions of the six cases. Cases 1–3 include fibre dispersion both along f0 and s0, whilst Cases 4–6 only include fibre dispersion along f0
Figure 5. The human LV model with 133,042 linear tetrahedral elements and 26,010 nodes (a), and the rule-based myofibre structure (b) that the myofibre angle varies from −60∘ at the epicardium to 60∘ at the endocardium. The two black lines (d and l) are used to determine radial expansion and longitudinal elongation ratios, respectively
Figure 6. Goodness-of-fit of the H-O model with the six cases as listed in Table 1
Figure 7. Myofibre stress distributions at end of diastole in the finite-element simulations when including the six cases of fibre dispersion. σσff and σσss denote the stress component along the mean myofibre direction f0 and the mean sheet direction s0, respectively. δσσff and δσσss are the differences comparing to Case 1
Figure 8. Comparisons of end-diastolic pressure-volume relationship computed by our models and experimental data from ex vivo human hearts [34]. Normal, healthy heart; ICM, ischemic cardiomyopathy; DCM, diopathic dilated cardiomyopathy; and LVAD, hearts supported by a left ventricular assist device
Figure 9. Differences of end diastolic stresses from Test 1 to Test 4 with varied myofibre and sheet rotation compared to Case 1. (left) myofibre stress differences; (right) sheet stress differences. Stress distributions for Case 1 can be found in Figure 7