Introduction

Myocardial dysfunction is a considerable social and economic burden because it can lead to heart failure due to repeated stresses and injuries [1]. However, our understanding of cardiac dysfunction remains incomplete. Mann and Bristow [2] suggested that heart failure could be viewed as a biomechanical model orchestrated by different subsystems, and the downstream biological remodelling is one of the drivers leading to myocardial dysfunction progression. Multi-scale and multi-physics biomechanical modelling of heart dynamics provides a unique way to gain insights of myocardial function and holds the potential for patient risk stratification and clinical decision making [3, 4].

Developing physiologically realistic heart models is complicated by many challenges [5], such as the complex anatomical geometry, nonlinear material responses of the myocardium, fluid-structure interaction, and issues across different length/function scales. In the last several decades, many nonlinear models of cardiac mechanics have been developed using the finite element (FE) method [6–9]. Different constitutive laws are introduced to characterize myocardium, such as the “pole-zero” constitutive law [9] and “Holzapfel-Ogden” law [10]. A recent review on heart modelling is provided in Ref [5]

In addition to the aforementioned structure-only ventricular models, immersed boundary (IB) method has also been used in modelling ventricular dynamics since the 1970s [11]. Recently, Griffith and Luo [12] have developed a hybrid IB with finite element method (IB/FE) that describes the immersed structure using the finite elements elasticity. Gao et al. [13] successfully applied the IB/FE framework to an imaging-derived subject-specific human LV model in diastole, as well as a diseased heart with myocardial infarction [14]. Later, Cai et al. [15] developed a more automatic procedure to construct mathematical models from medical images, and simulated the end-diastolic (ED) and end-systolic (ES) behaviours of a healthy left ventricle. Their results demonstrated that it is possible to simulate realistic ventricular dynamics using the IB/FE approach in conjunction with structure-based constitutive models [13]. However, the electrophysiological modelling of the heart was oversimplified in their models.

Myocardial contraction is highly dependent on the dynamics of calcium ion (Ca²⁺) in individual myocytes [16]. Numerious myocyte models exist in the literature [17], such as the model developed by Luo et al. [18], the FK model by Fenton et al. [19], and so on. Human ventricular myocyte models are emerging in recent years. For example, the ten Tusscher model [20], the OVVR model [21], and the Grandi-Pasqualini-Bers (GPB) model [22]. The electrophysiology models in single myocyte have proven to be very useful to study the mechanisms underlying disturbances of Ca²⁺ due to functional dysfunction [23] and the effects of drugs [24].

The failing heart undergoes continuous electrical-physiological, metabolic and structural remodelling in order to compensate the loss of functioning myocytes [2]. Ca²⁺ is the most important ion regulating myocardial excitation-contraction coupling [16]. Continuous remodelling of ion channels in failing myocyte, such as a significant increase in the late Na⁺ current [25], will cause AP duration prolongation and altered intracellular Ca²⁺ homeostasis, the hallmark characteristics of failing myocytes [26]. Several studies [27–29] have used the GPB model to explore excitation–contraction coupling and repolarization abnormalities in single human myocyte, in particular the Ca²⁺ handling process.

The multiscale nature of the heart requires one to model myocyte dynamics at cellular level and the electrophysiology at tissue level [17]. Many cardiac mechanical models, our previous models included, have not incorporated the detailed electrophysiology model for the myocytes. These models assumed a homogeneous and simultaneous contraction in the LV wall, and used prescribed intracellular calcium transient (CaT), or specified contraction not explicitly based on electrophysiology [7, 14, 30]. In this work, we coupled a more detailed model of myocyte ion dynamics for excitation-contraction coupling into the tissue-level mechanics to study how Ca²⁺ dynamics at the cellular level affects the myocardial function at organ level, in particular in failing heart. This was achieved by extending our IB/FE LV model [13, 14] by embedding a modified GPB model together with a myofilament model. The new coupled model was then used to study the myocardial active contraction at both cellular- and tissue- levels in healthy and failing states.

The remaining of the paper is organised as follows. In Methodology, we report how we use a modified GPB model to model the myocyte electrophysiology [29], and the filament model by Niederer et al. [31] (henceforth referred to as the NHS model) for the active contraction. This is followed by studying four different cases of myocyte contraction from healthy to failing states. In Results, we describe how we embed the active tension generation using the NHS and GPB model into a IB/FE LV model, to study the pump function in the healthy case and the worst failing case. We subsequently show how changes in cellular electrophysiology may affect the LV pump function at the organ level. The paper ends with Discussion and Conclusions.

Methodology

Myocyte dynamics at cellular level

Ion dynamics: The GPB model.

In this study, a modified GPB model by Cardona et al. [29] was used to simulate the Ca²⁺ dynamics in healthy and failing states. The reasons for choosing the GPB model were that (1) it matches experimental data well [22]; (2) it is adequate to analyse action potential with detailed Ca²⁺ dynamics, which plays a crucial role in excitation-contraction in myocardium [16]; (3) other studies have demonstrated that it can reproduce the electrical changes of failing myocytes caused by remodelling in ion channels and transporters [28, 29]. The original GPB model consists of 38 ordinary differential equations, which describe the kinetics of different ion channels. Details of the GPB model and parameters are given in S1 and S2 Appendices.

In the failing myocytes, the intracellular sodium concentration and Ca²⁺ handing are closely related. For example, the intracellular sodium concentration, especially I_NaL, is increased in failing myocytes, which contributes to cellular Ca²⁺ accumulation. The introduced I_NaL in the modified GBP model was modelled following the Hodgkin-Huxley formula [32], (1) (2) (3) and where g_Nal is the maximum conductance, m_L is the activation gate, h_L is the steady-state inactivation gate, τ_hl (233 ms) is the time constant of inactivation, E_junc and E_sl are, respectively, the Nernst potential for the sodium ion in the junctional and subsarcolemmal compartments (see S1 and S2 Appendices for details). Cardona et al. [29] found that the most impactful parameters on Ca²⁺ handing are the sarcoplasmic reticulum Ca-ATPase (SERCA) function, I_NaL, the Na/K pump current (I_Nak), the sarcoplasmic reticulum leak current (I_leak), the background Ca current (I_Cab) and the Na-Ca exchanger current (I_ncx). The modified GPB model was further extended by Trenor et al. [28], who introduced the impact factors to simulate the transition from healthy to failing states. The impact factors are coefficients that can change the parameters on Ca²⁺ handling in the modified GPB model. The final formulations of the modified GPB model employed in this study are given in S3 Appendix.

Based on the work of Trenor et al. [28], four cases were designed to represent different failing states; Case A: the healthy state, Cases B and C: two intermediate failing states, and Case D: the end-stage failure state. The impact factors of Case A and Case D were adopted from [28], and the impact factors of Case B and Case C were obtained by reducing the values of Case A by one-third and two-thirds of the difference between Cases A and D, as shown in Table 1. In healthy state, the peak value of the intracellular CaT obtained by this modified GPB model were between 0.3 μM and 0.4 μM. Thus, in order to make the peak CaT to be around 1 μM which was needed by the later active contraction model [31], two diffusion constants were changed (needing to note that changes of and may not be the ideal way to obtain an intracellular CaT profile with a peak value of 1 μM), they were (1) the diffusion constant between the junctional cleft and subscarcolemmal: L/ms, and (2) the diffusion constant between the subsarcolemmal and the bulk cytosol space: L/ms.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Impact factors in the four different cases based on the modified GPB model.

https://doi.org/10.1371/journal.pone.0174834.t001

Myofilament active contraction: The NHS model.

The active tension in each myocyte, T(t), was modelled as a function of the myocyte stretch (λ_f), the stretch rate (), and the intracellular CaT [31]. Details of the NHS model and the parameters are given in S4 Appendix. To study the influence of CaT on the active tension generation in individual myocytes, the intracellular CaT from the modified GPB model was coupled to the NHS model to trigger active contraction. We further assumed that in the initial time T = 0 kPa, and λ_f = 1, during myocyte contraction. Therefore, the simulated active tension in one myocyte corresponded to an isometric tension experiment at resting length.

Myocardial dynamics at organ level: The LV model

IB/FE formulation.

The in-house developed IB/FE framework was employed to simulate the LV contraction [12], in which the Lagrangian equations of the LV dynamics were approximated on an FE mesh, and the Eulerian equations of the fluid were approximated on a Cartesian grid [11]. A regularised delta function δ was used to describe the fluid-structure interaction, which implies that the IB/FE approach permits nonconforming discretization of the fluid and structure domains. In brief, the governing equations of the coupled fluid-structure system are (4) (5) (6) (7) where Ω⊂R³ denotes the physical region occupied by the fluid-structure system, and E⊂R³ denotes the region occupied by the immersed solid in the reference configuration. ρ is the fluid density, p is the pressure, and μ is the viscosity. x = (x₁, x₂, x₃) ∈ Ω denotes the Cartesian (Eulerian) coordinates, X = (X₁, X₂, X₃) ∈ E denotes the material (Lagrangian) coordinates in the reference configuration. χ(X, t) ∈ Ω gives the physical position of material point X at time t. Therefore, the physical region occupied by the structure at time t is Ω^e(t) = χ(E, t), and the physical domain occupied by the fluid at time t is Ω^f (t) = Ω − Ω^e(t). The IB/FE LV model is shown in Fig 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Schematic illustration of the LV model.

Ω is the total computational domain, Ω^e is the LV structure, f denotes the myofibre direction, s denotes the sheet direction, and n denotes the normal orientation of the plane spanned by the myofibre and sheet directions.

https://doi.org/10.1371/journal.pone.0174834.g001

The total Cauchy stress tensor of the coupled fluid-structure system was defined as (8) where σ^f = −pI + μ[∇u + ∇u^T] is the stress tensor of a viscous incompressible fluid, existing in both the solid and fluid regions. I is the identity matrix, and σ^e is the elastic stress tensor determined by the deformation of the immersed structure, defined as (9) where is the structural deformation gradient and J = det(F), P^e is the first Piola-Kirchoff stress tensor, consisting of the passive elastic response of the myocardium P^p and the active stress P^a, that is (10) P^p was calculated from a strain-invariant based strain-energy function (W) introduced by Holzapfel and Ogden [10], (11) and (12) where a, b, a_i and b_i (i = f, s, fs) are eight non-negative material parameters. I₁ = tr(C), and C = F^T F is the right Cauchy-Green deformation tensor. (i = f, s), , f = Ff₀, and , s = Fs₀. f₀ and s₀ are the fibre and sheet directions in the reference configuration, as shown in Fig 1. I_4f and I_4s are the squared stretches along the fibre and sheet directions, respectively. ensures that the myofibres only exert stresses when stretched. The coupling effects of the fibre and sheet directions are represented by . The last term in Eq (12) is to ensure that when F = I, P^p = 0.

The active contraction stress tensor in the current configuration was (13) where T(x, t) is the active tension based on the NHS model, triggered by the intracellular CaT from the modified GPB model, is the unit vector of the current myofibre direction. The NHS model parameters were initially obtained from animal experiments. To apply the NHS model to human heart, we introduced an extra scaling parameter T_scale in Eq (13) as in our previous study [14], (14)

LV model construction and implementation.

The LV model was constructed from an in vivo magnetic resonance imaging study of a healthy volunteer [15, 33], and a rule-based method was used to generate the myocardial fibre and sheet directions [8]. In the simulations, the physical domain occupied by the fluid-structure system (Ω) was taken to be a 15 cm×15 cm×20 cm box that was discretized with 96 × 96 × 128 Cartesian grids. An explicit version of crank Nicolson-Adams Bashforth scheme was used for time stepping, which required a relatively small time step size (1.22 × 10⁻⁴s). The LV model was implemented within the open-source IBAMR software framework (https://github.com/IBAMR/IBAMR).

The boundary conditions used were such that the longitudinal and circumferential displacements of the basal plane were constrained, whereas the radial displacements of the basal plane were set free. The no-slip condition was imposed on the left ventricular wall, and a zero pressure condition was applied on the physical domain of the fluid boundary. A spatially uniform pressure was loaded on the endocardial surface. Because it is difficult to measure the LV pressure in vivo, a population-based ED pressure (8 mmHg) was assumed [15]. In diastolic filling, the endocardial pressure was linearly increased from 0 to the assumed ED pressure. Then the myocardium started to contract with increased intracellular CaT and increased systolic pressure. To simulate the ES state, both the ES pressure (150 mmHg, obtained from the cuff measurements of the subject) and the peak value of intracellular CaT were maintained until a steady-state LV dynamics was reached.

The parameters used in the IB/FE LV model were: a = 0.24 kPa, b = 5.08, a_f = 1.46 kPa, b_f = 4.15, a_s = 0.87 kPa, b_s = 1.6, a_fs = 0.3 kPa, b_fs = 1.3 and T_scale = 3.0. These were obtained inversely by ensuring that the LV volume differences between the clinical measurements and the simulated values were less than 5% at end-diastole and end-systole [14].

Results

Myocyte contraction

The simulated AP profiles in individual myocytes in one cardiac cycle for the fours cases are showed in Fig 2(A), in which the spike notch dome morphology can be found. The characteristics of the AP profiles are summarised in Table 2. In Case A, the AP amplitude is 122.33 mV, which is comparable to published experimental measurements (100 mV [34], 132 mV [35], and 135 mV [36]). The maximal upstroke velocity in Case A is also within the range of experimental measurements (228 ± 11 V/s [37] to 446 ± 46 V/s [38]). The APD in Case D is much longer than that in Case A, similarly for the effective refractory period (ERP, measured as the interval between the AP upstroke and the repolarization level). The resting potential, the amplitude of the AP, and the maximal upstroke velocity in Cases B, C and D are much lower than the values in Case A.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Numerical results of myocyte contraction in a single myocyte.

Profiles of (A) action potential, (B) Intracellular CaT, and (C) Active tension.

https://doi.org/10.1371/journal.pone.0174834.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Characteristics of the active potential profiles, the intracellular CaT, and the active tension generation in a single myocyte.

https://doi.org/10.1371/journal.pone.0174834.t002

The comparison of intracellular CaT in the four cases are shown in Fig 2(B), and their characteristics in Table 2. The intracellular CaT in Case A peaks rapidly with the highest peak value compared to other cases. In contrast, the CaT in Case D peaks tardily at a much lower value, almost 63% less than the value in Case A. The intracellular CaT from the two intermediate cases (Case B and Case C) lie between Case A and Case D. After 0.2s until 0.6s, the CaT in Case D is slightly higher than others (Fig 2(B)). The maximal upstroke velocity is also highest in Case A and lowest in Case D.

Fig 2(C) shows the isometric active tension generation for the four cases, and the characteristics of the generated active tension are also summarized in Table 2. It can be found that the generated active tension and the maximal upstroke velocity decrease with the decrease of the intracellular CaT. In Cases C and D, the peak active tension is about one-third of the value in Case A, with a much lower upstroke velocity. Interestingly, the developed active tension in Case D after 0.5s is slightly higher than other three cases, which could partially be explained by the higher CaT in Case D after 0.2s as shown in Fig 2(B), however, the difference in active tension is not much.

LV contraction

The deformed end-diastolic LV geometries are same in Cases A and D because of the same passive material properties and boundary conditions used, shown in Fig 3(A). The deformed end-systolic geometry in Case A is shown in Fig 3(B). In Case A, the LV cavity volume is 145 mL at end-diastole and 61 mL at end-systole, with an ejection fraction (EF) of 57.55%, which are in agreement with the measured end-diastolic and end-systolic volumes (ED: 143 mL, ES: 64 mL). The deformed end-systolic geometry in Case D is showed in Fig 3(C), the corresponding LV cavity volume is 189 mL, and the mean myofibre strain is 0.07 ± 0.06. The larger end-systolic LV cavity volume and the positive systolic myofibre strain in Case D suggest that the LV model does not contract. This is largely because of the much lower level of intracellular CaT in Case D, which results in a much lower systolic active tension. The average systolic active tension in Case D is 50.2 ± 4.3 kPa, a 38% less than that of Case A (81.53 ± 21.0 kPa). The lower myocardial active tension in Case D agrees with the cellular level results: the isometric tension in Case D is 66% less compared to the peak value in Case A.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Deformed LV structures.

A: Deformed LV structures at end-diastole; B: the end-systolic state for Case A; C: the end-systolic state for Case D.

https://doi.org/10.1371/journal.pone.0174834.g003

Remodelling in myocytes, such as increasing myocardial contractility (i.e T_scale), decreasing the afterload (i.e. the end-systolic pressure), or increasing diastolic filling pressure (i.e. the ED pressure), is needed to maintain the minimum required pump function in a failing heart [2, 39]. We further investigated the effects of decreasing the end-systolic pressure and increasing T_scale on LV pump function in Case D, and the results are shown in Fig 4(A). As expected, the LV end-systolic volume decreases with the decrease of the end-systolic pressure or the increase of T_scale. Fig 4(B) shows the corresponding changes of LV stroke volume, defined as ED volume—ES volume. When the end-systolic pressure is 100 mmHg and T_scale is 6.0, the end-systolic volume is 116 mL and an EF of 19% is achieved in Case D, which still indicates a severe heart failure. The LV stroke volume in case A is 83.60 mL, almost four times greater. Fig 4(C) is the average systolic active tension in Case D. For a given end-systolic pressure, a higher T_scale is associated with a higher active tension and a stronger contraction. On the other hand, for a given T_scale, although a lower end-systolic pressure reduces the active tension generation, the LV stroke volume actually increases. This is because the active tension is length and strain-rate dependent, therefore a lower end-systolic pressure means that the myocardium is able to contract early and quickly, leading to a larger LV stroke volume. Fig 4 shows that only when the end-systolic pressure is below 120 mmHg, the LV can pump the blood into the systemic circulation. These results may suggest that it is insufficient to improve the pump performance by simply increasing the myocardial contraction, or reducing the end-systolic pressure.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. LV pump function in case D under different end-systolic pressures and T_scale.

A: The end-systolic volume; B: LV stroke volume; C: Average systolic active tension.

https://doi.org/10.1371/journal.pone.0174834.g004

Fig 5(A) shows the end-diastolic and end-systolic volumes under different end-diastolic pressures in case D for a given end-systolic pressure (= 100 mmHg) and T_scale(= 6.0). The end-diastolic volume increases with the increase of the end-diastolic pressure, which is consistent with the Frank-Starling law of the heart. On the other hand, the end-systolic volume does not change much. The LV stroke volume versus end-diastolic pressure is shown in Fig 5(B). The LV stroke volume and EF reach their maximum values (47mL, 29.19%) at end-diastolic pressure = 16 mmHg, but still less than the value in Case A. The results suggest that with increased end-diastolic pressure, the left ventricle can pump more blood, which is beneficial in terms of meeting body’s blood demand. However, studies have found that elevated diastolic filling pressure may associate with terminal heart failure in the longer run [40, 41].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. LV pump function under different end-diastolic pressures in case D.

A: The end-diastolic and end-systolic volumes; B: LV stroke volume.

https://doi.org/10.1371/journal.pone.0174834.g005

Discussion

From biomechanical perspective, myocaridal dysfunction could be considered as biomechanical model with interrelated and intricate changes and remodelling in cardiac structure and function as the result of downstream biological abnormalities [2]. Advanced biomechanical modelling of heart mechanics with multi-scale and multi-physics may hold the potential to improve our understanding of failing hearts and shed lights on effective treatments [3]. This computational work, based on an advance human myocyte electrophysiological model, a myofilament model for excitation-contraction coupling, and a multi-physics organ-level LV mechanics model, demonstrates that a multi-scale biomechanical model can be used to investigate the effects of downstream biological abnormalities on ventricular pump function. Once validated, it will provide a platform to understand the functional and structural remodelling in failing hearts.

Using this multi-scale LV model, we found that (1) profiles of AP, intracellular CaT and active tension are very different at cellular level from healthy to failing states; (2) a lower CaT can lead to a much less active tension generation; and (3) a lower intracellular CaT in individual myocytes can result in LV pump function failure even with increased myocardial contractility, decreased afterload, and increased diastolic filling pressure. Our modelling results at both the cellular and tissue levels for one healthy and diseased LV cases showed good agreement with experimental measurements and earlier studies, additional comparison is provided in S1 Table.

Specifically, the APD in the failing case (D) is prolonged compared to the value in the healthy case (A) because of the much slower rate of repolarization in Case D, this is consistent with experimental findings [42, 43]. As the myocardial dysfunction worsens (from Case A to Case D), the peak and the maximal upstroke velocity of intracellular CaT decreases, and the curve gradually flattens. With decreased intracellular CaT and its upstroke velocity, the generated active tension also decreases, to almost 66% in Case D. Although the CaT in Case D is higher than other cases after 0.2s, the active tension in Case D is generally lower than in Case A. This suggests that in a failing state, the remodelling in ion dynamics (by prolonging APD and increasing CaT in the later contraction phase) has limited effect on active tension generation. Therefore, our study suggests that a well-balanced physiological Ca²⁺ handling mechanism is an essential component in maintaining normal myocyte contractile function, which agrees with [44].

By simulating a pseudo-isometric experiment in a myocyte, we found that the active tension from Case D is much lower compared to Case A. To see how this affects the LV performance, we embedded the modified GPB model into a dynamic LV model [14]. Our results show that the reduction in average systolic active tension across the whole LV in Case D is 39% compared to the Case A, and no blood can be pumped out (the LV cavity volume at end-systole is greater than its end-diastolic volume). This is when sudden death occurs.

A decline in LV pump function will activate various compensatory mechanisms to restore a normal homeostatic cardiovascular function [2, 39]. These include a higher heart rate, myocyte hypertrophy, and increased LV end-diastolic volume (i.e. through elevated diastolic filling pressure). Treatments of heart failure typically aims to preserve the cardiac output by reducing the blood pressure (the afterload) using intravenous vasodilators, or increasing the myocardial contractility using positive iontropes [45, 46]. Indeed, our results show that decreased end-systolic pressure and increased myocardial contractility are not enough to maintain normal LV pump function. For instance, with the end-systolic pressure of 100 mmHg, the LV model can pump the blood out only when T_scale is above 4. However, the pump function is compromised even when T_scale is increased to 6, the LV stroke volume is only 29 mL, which is far less compared to the healthy Case A (83 mL). Moreover, continuous increase of the myocardial contractility can deteriorate the LV function in the longer term leading to cardiac de-compensation. It is perhaps for this reason that some inotropic therapies increased, rather than decreased, mortality [47].

In addition, the capacity of decreasing afterload or increasing myocardial contractility is not without limitation, therefore, other compensation mechanisms will be activated, such as elevated diastolic filling pressure through the Frank-Starling law of the heart [39, 40]. In case D, when the end-systolic pressure is 100 mmHg and T_scale is 6, the LV stroke volume increases with the increase of end-diastolic pressure. The LV stroke volume reaches 48 mL when the end-diastolic pressure is 16 mmHg. It is expected that with further increased end-diastolic pressure, LV stroke volume will increase further. However, clinical observations suggest that the elevated diastolic filling pressure would reflect patients at increased risk of developing late clinical symptoms of heart failure in the long term [41], even though it can increase the LV stroke volume in the short term as suggested by our models.

There are two recent mathematical descriptions of human ventricle myocytes, the GPB model [22] and the OVVR model [21]. Both models can match experimental data well though discrepancies exist in cellular and tissue levels [27, 48]. In this study the GPB model was chosen because it is adequate to analyse action potential, Ca²⁺ dynamics. Furthermore, the GPB model has been demonstrated to reproduce the electrical changes of failing myocytes due to ion channel and transporter remodelling [28, 48]. However, the issue with linking the GPB model to the NHS model is that the peak CaT in the GPB model is less than 1μM, which is not compatible with the active tension generation model in the multi-scale LV model. Thus, we had to increase two diffusion constants so that the peak CaT is around 1μM. OVVR model, on the other hand, has a CaT amplitude of about 1μM. The adjustments of J_Cajuncsl and J_Caslmyo may not be the ideal way to tune the GPB model to obtain the required CaT profile, a carefully designed calibration procedure is required. Alternatively, one can use a different myofilament model to make the two models compatible since the NHS model was developed based on animal experiments. Despite of this, our modelling results from the adopted GPB model agree well with prior studies and experiments, for example, in the failing state, the CaT amplitude decreases, its rise and decay rates are slowed down [49], additional comparisons with other studies are provided in a supplementary table. We expect it will work equally well by replacing the GPB model with the OVVR model in our model, and it will be interesting to compare the behaviour of both models in a multi-scale biomechanical myocardial model, and explore other myofilament models developed for human.

The remodelling of the failing heart occurs in many aspects at both cellular and organ levels, in terms of structure [50], electrophysiology [51], metabolism [52], and couplings among them [16]. Cardiomyocytes in healthy and failing hearts differ in numerous ways, such as prolongation of AP, altered Ca²⁺ handing, intracellular Na⁺ accumulation, and impaired electro-mechanical coupling [16, 44, 51]. It is well known that failing myocytes have abnormal CaT, with decreased amplitude and slowed kinetics [49]. Except for the functional remodelling of ion channels and pumps as described in the adopted GPB model, cellular structural remodelling is also recognized as one of the underlying causes for altered intracellular Ca²⁺ homeostasis. These include alternations in the T-tubule structure [53], micro-architecture changes in sarcoplasmic reticulum [54], and alterations in molecular and biochemical structure of myofibrils [55]. Furthermore, metabolic remodelling also contributes to the impaired excitation-contraction coupling due to ATP deficiency [52]. In this study, we only studied the effects of CaT from different failing stages on myocardial active contraction that is caused by ion channel and transporter remodelling [28, 48]. However, we expect that other types of remodelling mechanism can be ready implemented in our model by including more detailed descriptions of related remodelling process.

Similar as the GPB model used here, the structural remodelling of the myocyte is not modelled in the NHS model [31]. An anatomical multi-scale myocardial description from sub-cellular to organ levels is needed in order to capture the spatially varied structural remodelling. However, there are many challenges in following this approach [56]. Even if one can develop detailed mathematical description based on proper upscaling, to solve this a model also requires efficient numerical methods and high-performance computers.

Other limitations of the work include: (1) we have not thoroughly analysed the effects of individual ion channel remodelling on myocardial active contraction, but focused on how CaT profiles from different failing stages affect LV pump function; (2) we assumed that the LV contraction was spatially homogeneous and temporally simultaneous, which has been widely adopted in heart simulations [6]. This treatment has limitations that the LV model will not be applicable to ischemic myocardial dysfunction, i.e. myocardial infarction, but more relevant to nonischemics myocardial dysfunction, such as ventricular fibrosis; (3) we also only simulated the end-diastolic and end-systolic phases of LV dynamics, not the whole cardiac cycle which involves interactions of the heart valves. Therefore, the developed multi-scale myocardial contraction model should be considered as a phenotype representation of failing hearts, rather than a personalized model. Nevertheless, this computational work provides a multi-scale platform for studying the effects of myocardial remodelling on LV pump function in a biomechanical perspective.

Conclusion

In this study, we have modelled myocardial excitation-contraction coupling both at cellular and organ levels by incorporating a modified GPB model and the NHS model at healthy and failing states. The LV model was implemented in an in-house developed IB/FE framework. Results show that the active tension decreases with the decrease of intracellular CaT both at cellular and organ levels. In a failing heart, decreased systolic blood pressure, enhanced contractility and elevated diastolic filling pressure can improve heart pump function in the short term (i.e. increased stroke volume), however, they may lead to terminal heart failure in a longer term if a balanced physiological calcium ion handling mechanism is not restored. This study forms part of the continuous effort in multi-scale electro-mechanics modelling for clinical problems, and provides a platform to study remodelling in a failing heart from a biomechanical perspective.

Supporting information