Bridging Continuum and Molecular Models in a Beating Heart Simulation

By Takumi Washio

Numerical continuum modeling of muscle offers a very interesting research topic. In typical continuum materials, the strain uniquely determines the stress. However, this rule does not apply to muscles undergoing contraction. The actomyosin complexes in muscle [1] constantly change their conformations during the contraction phase of events, such as attachment, dissociation, and power strokes of the lever arm (see Figure 1a). Furthermore, the transition rates depend on the strain of the lever arm, which is affected by the continuum muscle contraction.

In computer simulations of a beating heart (see video below), modeling these molecular behaviors is necessary to correctly reproduce the muscle contraction. The contractile force in cardiac cells rises when the sarcoplasmic reticulum (SR) releases Ca²⁺ ions, and decreases when these ions are sequestered back into the SR. In many earlier studies, the temporal change of contractile force under a given Ca²⁺ transient was computed by a system of ordinary differential equations (ODEs). However, ODEs cannot easily model the stochastic and cooperative behaviors of actomyosin complexes. For example, in a healthy beating heart, the left ventricular pressure (LVP) (see Figure 1b, black line) falls to almost zero at the end-systole (or end of contraction), but nearly 10% of the peak Ca²⁺ concentration remains in the cytosol, that is, within heart cells (see Figure 1b, red line). Furthermore, the LVP falls much more rapidly than the Ca²⁺ concentration. These quick relaxation properties of cardiac muscle are assumed to originate from the cooperative behavior between the neighboring actomyosin complexes and the stochastic power-stroke, strain-dependent mechanism of the myosin lever arms. Unless we correctly model these molecular properties, we fail to reproduce the quick relaxation of the cardiac muscle. The consequence is insufficient blood-filling into the ventricular cavity during the diastolic phase.

Figure 1. Multiscale modeling. 1a. MC model of the actomyosin complex. The contraction force is contributed by three attached states (XB_PreR, XB_PostR1, and XB_PostR2). The transitions between N_XB and P_XB are influenced by the states of the T / T unit above the myosin head through the coefficients K_np,K_pn and the states of the neighboring myosins through the cooperative factors γⁿ, γ^-n (γ=40). 1b. Transients of the averaged Ca²⁺ concentration over the ventricle (red line) and the left ventricular pressure (LVP) (black line). 1c. Sarcomere model composed of actomyosin complexes and FEM ventricle model showing the twisted fiber orientations along the transmural line. Figure credit: Takumi Washio.

Multiscale Model

Failure to successfully simulate a beating heart with existing ODE models has motivated our method to directly mimic individual molecules in our beating heart model. To realize this idea, we couple macroscopic continuum dynamics by the finite element method (FEM) with microscopic molecular dynamics by the Monte Carlo (MC) method. These two approaches differ on both spatial and temporal scales. In our approach [2], we embed the sarcomere model of actomyosin complexes into each tetrahedral element of the FEM (see Figure 1c) along the fiber orientation. Here we assume that a single sarcomere force represents the contractile tension generated by all cardiac cells in the element. Usually we set the time step as \(\Delta T = 1\) ms in the FEM analysis, and subdivide this value into a finer time step \(\Delta t = 10 \mu s\) for the MC analysis.

Coupling of Molecules and Continuum

In our approach, the macroscopic active stress tensor \(S_{act}\) in the finite element time interval \([T, T, + \Delta T]\) is implicitly determined, ensuring compatibility between the molecular and continuum virtual works, as shown in Figure 2. To account for the state transitions of the actomyosin complexes during this time interval, we relate the strain rates on three scales (actomoysin complex, sarcomere, and continuum). In the FEM analysis, the second Piola-Kirchhoff stress tensor represents the stress under contractile tension \(T_f\):

\[S_{act}= \frac{T_f}{\lambda} \mathbf{f} \bigotimes \mathbf{f} = \displaystyle {\sum\limits_{i,j}} \frac {T_f}{\lambda} f_i f_j \mathbf{\mathit{e}}_i \bigotimes \mathbf{\mathit{e}}_j, \tag1 \]

where the unit vector \(\mathbf{f}= \sum\limits_{i=1}^3 f_i \mathbf{e}_i\) denotes the fiber orientation at the material point \(\mathbf{\mathit{X}}\), and \(\lambda\) denotes the stretch along the fiber orientation \(f\), given by

\[\lambda = \Bigg{\|} \frac{\partial \mathbf{x}}{\partial \mathbf{X}} \mathbf{\mathit{f}} \enspace \Bigg{\|} = \sqrt {\displaystyle {\sum\limits_{i, j, k-1}^3} \frac{\partial x_k}{\partial X_i}  \frac{\partial x_k}{\partial X_j} f_i f_j}. \]

The insertion of \(\lambda\) into the denominator of \((1)\) is justified by the infinitesimal relationship \(T_f \delta \lambda = \mathbf{S}_{act} \cdot \mathbf{\delta E}\), where \(\mathbf{E}\) is the Green-Lagrange strain tensor. The contractile tension \(T_f\) is computed by summing up the molecular forces produced by the actomyosin complexes and arranged along the actin filaments in the sarcomere model:

\[T_f = \frac{R_S}{SA_0} \frac{2}{ns} \sum\limits_{is=1}^{ns} \sum\limits_{im=1}^{nm} F_M(im, is).\]

Here, \(SA_0 (=1000 \enspace \textrm{nm}^2)\) is the cross-sectional area of a single actin filament; \(R_S\) denotes the volume ratio of the sarcomere \((=0.5)\); \(nm(=38)\) is the number of myosin molecules surrounding the binding sites, arranged along one of the two spirals in the actin filament; and \(ns\) is the number of actin filament samples in the sarcomere model. The force generated by an individual myosin molecule is given by

\[F_M(im, is) = \frac{\Delta t}{\Delta T} \sum\limits_{k=1}^{nt} \delta_{A,k} (im, is) K_M^{T + k\Delta t} x (im, is),\]

where the FEM time step interval \([T, T+\Delta T]\) is divided into \(nt\) MC time steps. \(\delta_{A,k}\) takes \(1\) in the attached case and \(0\) in the detached case for \(k=1,..,nt\), while \(k_M\) denotes the spring constant of the myosin lever arm and \(^{T+k\Delta t}x\) is the strain of the myosin lever arm at time \(T+k\Delta t\). This strain is a function of the initial strain at the attachment \((x_{init})\) under the thermal fluctuation, the power stroke from the latest attachment from \(t_A\) to \(T+k\Delta T (^{T+k\Delta t}x_{PS})\), and the length change due to the filament sliding:

\[^{T+k\Delta t}x(im, is)=x_{init}+^{T+k\Delta t}x_{PS}+\frac{SL_0}{2} \left \{ {\int_{\textrm{min}(t_A,T)}^{T}} ^{t}  \dot{\lambda}dt+ ^{T+\Delta T} \dot{\lambda} (T+k\Delta t - \textrm{max}(t_A, T)) \right \}. \tag2 \]

Here, \(SL_0 / 2 \) is the unloaded half-sarcomere length that relates the muscle stretch rate along the fiber orientation \((\dot{\lambda})\) to the shortening velocity of the sarcomere \((-SL_0\dot{\lambda}/2)\).

Figure 2. Strain rates in the model (top) and the virtual works (bottom) on the three scales: the actomoysin complex, the sarcomere, and the continuum. Figure credit: Takumi Washio.

The working stroke size \(s\) increments \(x_{PS}\) during the forward transition from the pre-power stroke to the post-stroke state, and decrements \(x_{PS}\) in the reversal stroke from the post-stroke to the pre-stroke state. The rate constants of the forward \((f)\) and backward \((b)\) transitions between the pre and post-stroke states are determined as follows, with the relationship given by the statistical equilibrium:

\[\frac{f(x)}{b(x)}= \textrm{exp} \Bigg ( - \frac {E_{Post}-E_{\textrm{Pr}e}+k_M\Big((x+s)^2 -x^2 \Big) /2}{kT} \Bigg), \tag3 \]

where \(k\) and \(T\) denote the Boltzmann constant and the temperature respectively, and \(E_{\textrm{Pr}e}\) and \(E_{Post}\) are the free energies of the myosin head in the respective pre and post-stroke states. The difference \(E_{Post}-E_{\textrm{Pr}e}\) corresponds to the partial transfer of the chemical energy obtained by ATP hydrolysis to the mechanical stress energy generated by the strain increment \(s\). At the current FEM step in \([T,T+\Delta T]\), the filament sliding contribution (third term on the right-hand side of \((2)\)) is computed by implicitly assuming the stretch rate \(\dot{\lambda}\) at \(T+\Delta T\). Through this implicit scheme, we can correctly incorporate the temporal stiffness into the Jacobian matrix and stabilize the Newton-Raphson (NR) iteration in the FEM analysis within a reasonable time step \(\Delta T\). Note that the NR steps iteratively reuse the computational results obtained in the MC steps, in which \(\dot{\lambda}\) at \(T\), not at \(T+ \Delta T\), is referred to compute the strain \(x\) of \((2)\).

Beating Heart Simulation Results 

Figure 3b shows the behavior of the sarcomere model embedded in the inner layer of the left ventricular free wall. At the end-systolic phase, or end of contraction  (range surrounded by the broken lines in Figures 3a and 3b), the relative current of the backward to forward transitions increases as the sarcomere shortening decelerates and the sarcomere contraction finally ends. A quick stretch immediately follows the increased frequency of backward transitions. In cooperation with neighboring myosin molecules, these backward transitions swiftly reduce the blood pressure (see Figure 3a, red line), facilitating the quick filling of blood into the left ventricle (see Figure 3a, black line). Figure 3c presents the contours of the arm strain distribution. In the post-stroke state \((\textrm{XB}_{\textrm{PostR}2})\), the center of the distribution shifts to larger strains towards the end-systolic phase. Such a distribution shift has a major impact on balancing the frequencies of the forward and backward transitions determined by \((3)\). 

Figure 3. Simulation results of a beating heart. 3a. Time transients of the left ventricular pressure and volume. 3b. Sarcomere force and length, and the currents of the forward and backward transitions. 3c. Stain distributions in the lever arms of the sarcomere model embedded in the inner layer. Figure credit: Takumi Washio.

We are currently developing more detailed molecular models under the post-K supercomputer project.¹ This study will provide insights into the molecular-level mechanisms that control the state transitions, which are necessary for linking various mutants of contractile proteins with heart failures.

¹ This work is supported in part by MEXT as Strategic Programs for Innovative Research Field 1 Supercomputational Life Science and a social and scientific priority issue (integrated computational life science to support personalized and preventive medicine), to be tackled using a post-K computer.

References 
[1] Sweeney H.L., & Houdusse, A. (2010). Structural and Functional Insights into the Myosin Motor Mechanism. Annual Review of Biophysics, 39, 539-57.
 [2] Washio, T., Yoneda, K., Okada, J., Kariya, T., Sugiura, S., & Hisada, T. (2016). Ventricular fiber optimization utilizing the branching structure. International Journal for Numerical Methods in Biomedical Engineering, 32(7), e02753.

Takumi Washio has worked at the Computational Biomechanics Laboratory School of Frontier Sciences at the University of Tokyo since 2004. He has engaged in development of the multiscale and multiphysics heart simulator. In 2015, his group started a company, UT-Heart Inc., to contribute to medical and clinical communities by applying their heart simulator.