Coupling Langevin Dynamics With Continuum Mechanics: Exposing the Role of Sarcomere Stretch Activation Mechanisms to Cardiac Function

High-performance computing approaches that combine molecular-scale and macroscale continuum mechanics have long been anticipated in various fields. Such approaches may enrich our understanding of the links between microscale molecular mechanisms and macroscopic properties in the continuum. However, there have been few successful examples to date owing to various difficulties associated with overcoming the large spatial (from 1 nm to 10 cm) and temporal (from 1 ns to 1 ms) gaps between the two scales. In this paper, we propose an efficient parallel scheme to couple a microscopic model using Langevin dynamics for a protein motor with a finite element continuum model of a beating heart. The proposed scheme allows us to use a macroscale time step that is an order of magnitude longer than the microscale time step of the Langevin model, without loss of stability or accuracy. This reduces the overhead required by the imbalanced loads of the microscale computations and the communication required when switching between scales. An example of the Langevin dynamics model that demonstrates the usefulness of the coupling approach is the molecular mechanism of the actomyosin system, in which the stretch-activation phenomenon can be successfully reproduced. This microscopic Langevin model is coupled with a macroscopic finite element ventricle model. In the numerical simulations, the Langevin dynamics model reveals that a single sarcomere can undergo spontaneous oscillation (15 Hz) accompanied by quick lengthening due to cooperative movements of the myosin molecules pulling on the common Z-line. Also, the coupled simulations using the ventricle model show that the stretch-activation mechanism contributes to the synchronization of the quick lengthening of the sarcomeres at the end of the systolic phase. By comparing the simulation results given by the molecular model with and without the stretch-activation mechanism, we see that this synchronization contributes to maintaining the systolic blood pressure by providing sufficient blood volume without slowing the diastolic process.

Introduction

With the advances in computational science made possible by improvements in hardware technology, it is now possible to create multi-scale simulation models of the heart in which the macroscopic behaviors of the beating heart can be reproduced and analyzed based on molecular mechanisms of the excitation-contraction coupling process (Kerckhoffs et al., 2007; Gurev et al., 2011; Sugiura et al., 2012). These models are based on many studies of cell models of cardiac electrophysiology (Luo and Rudy, 1994; ten Tusscher et al., 2004; Grandi et al., 2010). We also note that tissue modeling has provided deep insights into the nature of coupling and other interactions among cells in the heart wall (Clayton et al., 2011). Central to these in silico heart studies is an accurate model of crossbridge kinetics, which not only forms the basis of cardiac mechanics, but also has clinical relevance in the light of the many reports showing the involvement of sarcomeric proteins in the pathogenesis of cardiomyopathies (Cahill et al., 2013).

Ideally, a molecular dynamics simulation of actomyosin should be coupled with a macroscopic finite element model of the heart because with such a model the impact of mutations in the myosin molecules on cardiac function can be directly assessed. However, it is not possible to perform such simulations even with the best available high-performance computers, and current multi-scale heart simulators usually adopt state-transition models of crossbridge cycling. In these models, the rate constants for transitions between states are governed by the energy of each state (Huxley and Simmons, 1971), but the minimum in the energy landscape corresponding to each state ignores its width in the infinitely-sharp minimum approximation, in which the angle of each lever arm is fixed in the most stable configuration. Obviously, this is a simplification of the behavior of real myosin molecules experiencing thermal fluctuations, and we have recently reported that a model with an energy landscape possessing wide minima can reproduce experimental findings with higher accuracy (Marcucci et al., 2016). However, in that paper, we only examined simple Langevin dynamics with a single variable representing the free energy potential during the power stroke, and solved it using a Monte Carle (MC) simulation. In that case, the Kramers-Smoluchovski approximation (Gardiner, 2004) was used to obtain the rate constants of the transitions between the multiple states, which were given by discretizing the one-dimensional range of the angles of the lever arms. If we try to formulate a more realistic free energy potential as a function of multiple variables, the number of MC states increases explosively, and it is no longer possible to find the rate constants between the MC states theoretically. Therefore, it is desirable to establish a numerical scheme that directly couples the Langevin dynamics of the molecules with the macroscopic continuum dynamics.

Here, we report a novel numerical method to couple the microscale simulation of crossbridge kinetics described by the Langevin equation with the macroscopic mechanics simulations using the finite element method, even though the time scales differ considerably. In this method, the time step of the macroscopic model is set at a multiple of that from the microscopic model to reduce computational overhead. The validity of the method was confirmed with a comparison of the simulation results with the recently reported experimental findings on the spontaneous oscillation of cardiac sarcomeres (Ishiwata et al., 2011), which can be reproduced only by correctly handling the coupling of the motion of the sarcomeres with the actomyosin dynamics. By applying this method, we also show that a trapped crossbridge mechanism greatly facilitates ventricular function through the stretch-activation of the cardiac muscle (Stelzer et al., 2006). A notable feature of the stretch-activation is a long-lasting increase in the contractile tension after a small, rapid stretch is applied during activation. In the usual stretch-activation experiments, the stretch is 1% of the muscle length, which closely corresponds to the microscale size of lever arm swing (10 nm). It is likely that the rapid stretching induces an unusual persistent conformational change of the bound myosin molecules. In this work, we introduce a free energy potential for the power-stroke model in which some of the bound myosin molecules become trapped in a deformed conformation when a rapid stretch is applied. These trapped myosin molecules cannot recover under normal thermal fluctuation unless their rods become relaxed or extremely stretched by subsequent sarcomeric movements. Through the beating-ventricle simulations, we show how this mechanism contributes to improved blood circulation.

Materials and Methods

Our strategy of coupling the different scales is summarized in Figure 1. The stretch rates were transferred from the macro- to micro-scale while the contractile forces were transferred back from the micro- to macro-scale. Finite element continuum mechanics were applied to the ventricle model. The half-sarcomere model of actomyosin complexes was imbedded into each tetrahedral element of the finite element ventricle model along the fiber direction. The molecular variables that represent the deformation of bound myosin molecules were computed by the Langevin dynamics. The shortening rate $-\dot{\lambda }$ along the fiber direction in the ventricle model was transferred to the sarcomeric shortening velocity ż by scaling with the unloaded half-sarcomere length SL₀/2. The sarcomeric shortening velocity ż was applied in the actomyosin model to slide the myosin thick filament. The contractile force of the half-sarcomere model was given by the sum of the tensile forces of the bound myosin rods. The contractile force in the half-sarcomere model was transferred to the macroscopic contractile tension along the fiber direction. In our coupling approach, the computational time step size ΔT of the sarcomeric dynamics and the ventricle continuum dynamics is given by an integer multiple of the time step size Δt of the actomyosin Langevin dynamics (ΔT = nΔt) to reduce the computational and communication overheads. As will be discussed in section Multiple Time Step (MTS) Method, such a multiple time-step strategy can be applied without suffering numerical instabilities by also transferring the stiffness given by the bound myosin rods. Readers who are not interested in the numerical schemes may skip sections Multiple Time Step (MTS) Method and Coupling With the Finite Element Ventricle Model.

FIGURE 1

Figure 1. Coupling strategy for three scales. In the actomyosin system, x and ξ are variables representing the deformation of the bound myosin molecule. In particular, ξ is the strain of the myosin rods, and W_rod(ξ) is its strain energy. These variables were updated by the time step Δt ~ 0.25 ns, while the variables in the half-sarcomere and the ventricle models were updated by the time step ΔT = nΔt ~ 1 μs. The shortening of the half-sarcomere model is represented by the variable z, which affected the Langevin dynamics of the bound myosin molecules through the constraint condition: Δξ = Δx − Δz, while the sarcomeric contractile force ^{T, ΔT}F on the time interval [T, T + ΔT] was given by the sum of tensile forces ^T+kΔtF_i,j of the bound myosin rods averaged over the time interval (k = 1, ⋯ , n). The half-sarcomere model of actomyosin complexes was imbedded into each tetrahedral element of the finite element ventricle model in the reference configuration along the fiber direction f. The deformation at the time T of the ventricle is represented by the current position x = ^Tx(X) of the material point X, thus λ = ||∂^Tx/∂X · f|| is the stretch along the fiber orientation direction. This stretching was transferred to the shortening of the imbedded half-sarcomere model with the factor −SL₀/2, while the contractile tension T, ΔT T_f was given along the fiber direction by scaling the sarcomeric contractile force T, ΔTF by taking the cross-sectional area per thin filament (SA₀), and the volume ratio of the sarcomere within the ventricle wall (R_S) into account.

Langevin Dynamics of a Single Sarcomere

The parameters adopted for the molecular dynamics are summarized in Table 1. Here, the dynamic equations for a half-sarcomere model composed of n_F pairs of thick and thin filaments (Figure 2A) are introduced. In this half-sarcomere model, we assumed that the right ends of the thin filaments were connected to the Z-line, which was fixed in microscopic space. The shortening displacement of the left end of the thick filament from the unloaded position was denoted by z. On each thick filament, there were n_M myosin molecules, which underwent repeated attachment and detachment with the thin filament. The value of n_M = 38 was adopted from our previous work (Washio et al., 2016). During the attached phase, the lever arm (LA) of the myosin molecule rotated around the joint point “O” of the myosin head (MH) under a given free energy potential φ with additional random forces (Figure 2B). These rotations were either the power stroke or the reversal stroke, depending on the rotational direction. To represent the deflection of the LA, it was decomposed into two rigid components, LA1 and LA2, jointed at the point “P” (Figure 2B). As with the real structure of a myosin molecule, LA1 may contain a series of subdomains from the lower 50 kDa to the converter in the motor domain because some conformational changes of these parts were supposed to be accompanied by lever arm rotation. The displacement of the point “P” of the filament direction given by the rotation of LA1 from its pre-power stroke position around the joint point “O” was represented by y. Here, the conformation of the myosin molecule just after attachment was assumed to be the same as the pre-power stroke conformation. Similarly, the displacement about the joint point “Q” with the myosin rod was represented by x. Thus, θ = y − x was the deflection of the LA from the pre-power stroke conformation. The strain energy of the myosin rod was given by a function W_rod(ξ), where ξ was the strain (length change) in the filament direction from its unloaded natural length ξ₀. The rod strain energy was non-linear with the generated force, as with our previous work (Washio et al., 2016) for a rod with ξ < 0. For positive strain (ξ > 0), a constant stiffness with a spring constant 2.8 pN/nm was used (Figure 3A). Under these assumptions, the dynamics of the sarcomere was described by the following Langevin equations, where the suffixes i and j represent the indexes of the MHs and the thick filaments, respectively. Also, t is the time, and ^tδ_A,i,j is set to one if the MH was attached at time t to the thin filament, and zero otherwise.

$\{\begin{array}{l}{{\gamma }_X}^t{\dot{x}}_{i,j}+\frac{\partial \phi }{\partial x}{(}_{}^t{x}_{i,j}{,}_{}^t{y}_{i,j})+\frac{d{W}_{rod}}{d\xi }{(}_{}^t{\xi }_{i,j}){-}^t{R}_{X,i,j}=0\hfill \\ {{\gamma }_Y}^t{\dot{y}}_{i,j}+\frac{\partial \phi }{\partial y}{(}_{}^t{x}_{i,j}{,}_{}^t{y}_{i,j}){-}^t{R}_{Y,i,j}=0{\text{ , }}_{}^t{\delta }_{A,i,j}=1(1\le i\le n_M,1\le j\le n_F)\hfill \\ {}_{}^t{\xi }_{i,j}{-}^{t_{A,i,j}}{\xi }_{i,j}-{(}_{}^t{x}_{i,j}{-}^{t_{A,i,j}}{x}_{i,j}){+}^{t}z{-}^{t_{A,i,j}}z=0\hfill & (1)\end{array}$

${{\gamma }_D}^t{\dot{\xi }}_{i,j}+\frac{d{W}_{rod}}{d\xi }{(}_{}^t{\xi }_{i,j}){-}^t{R}_{D,i,j}=0,{{}^t\delta }_{A,i,j}=0\left(1\le i\le n_M,1\le j\le n_F\right)(2)$

$\begin{array}{ll}{{\gamma }_Z}^t\dot{z}+{K_Z}^{t}z-\frac{1}{n_F}{\displaystyle \sum _{j=1}^{n_F}{\displaystyle \sum _{i=1}^{n_M}t{\delta }_{A,i,j}\frac{d{W}_{rod}}{d\xi }}}\left({{}^t\xi }_{i,j}\right)=0& (3)\end{array}$

TABLE 1

Table 1. Parameters for the actomyosin dynamics.

FIGURE 2

Figure 2. The half sarcomere model (A) and the attached myosin molecular model composed of the myosin head (MH) and the lever arm (LA) and the rod (B). In the half-sarcomere model, the Z-line was fixed and the shortening distance of the left edge of the thick filaments is denoted by z. The configuration of the attached myosin molecules is represented by two variables x and y. The LA is decomposed into LA1 and LA2 to represent its deflection around the point “P.” The degree of the deflection is given by θ = y − x. The power stroke is given by the counter clockwise rotation of LA1 around the point “O.” The rod is a non-linear spring connecting the thick filament and the point “Q” of LA2. The strain of the rod is denoted by ξ.

FIGURE 3

Figure 3. Strain energy W_rod(ξ) of the myosin rod and the force given by its derivative (A). The details are described in Supplementary Material S1.2. The free energy landscape of the myosin molecule φ in the attached state with respect to the molecular variables y and θ (B). This free energy consists of φ_PS(θ, y) and W_LA(θ). φ_PS(θ, y) is the energy source of the power stroke (rotation of LA1). W_LA(θ) is the deformation energy of LA for its deflection. The pre-power stroke configuration corresponds to the local minimum at y = 0. The other two local minima at y = s₁ = 5.5 nm and at y = s₁ + s₂ = 11 nm correspond to the states after the first and second power strokes, respectively. Between the pre-and the post first power stroke states, the high energy barrier is assumed for the positive deflection (θ > 0) of LA.

Here, the probabilistic rules for transitions between the attached and detached states will be given below. At the time of attachment, the myosin molecule was assumed to be in the pre-power stroke state.

$\{\begin{array}{l}{{}^{t_{A,i,j}}x}_{i.j}=x_{\text{pre}}\equiv 0\hfill \\ {{}^{t_{A,i,j}}y}_{i.j}=y_{\text{pre}}\equiv 0\hfill & (4)\end{array}$

Here, t_A,i,j is the time at which the attachment occurred. The spring strain ^tξ_i,j was continuously updated at the transitions.

In Equations (1–3), γ_X, γ_Y, and γ_D were the damping coefficients, and ^tR_X,i,j, ^tR_Y,i,j, and ^tR_D,i,j were the random forces, which fulfilled the condition:

$\begin{array}{l}\{\begin{array}{l}\langle {}^t{R}_{\alpha ,i,j}\rangle =0\hfill \\ \langle {}^t{R}_{\alpha ,i,j}{,}^{t^{\prime }}{R}_{\beta ,k,l}\rangle ={\delta }_{\alpha \beta }{\delta }_{ik}{\delta }_{jl}{\sqrt{2{\gamma }_{\alpha }{k}_{B}T}}^{t-t^{\prime }}\delta \hfill \end{array},\\ \alpha ,\beta =X,Y,D,1\le i,k\le n_M,1\le j,l\le n_F& (5)\end{array}$

where Boltzmann's constant is k_B and the temperature is T. In this paper, the damping coefficient γ_D was set to 70 pN · ns/nm, following Howard (2001), while γ_X and γ_Y were set to 20 and 50 pN · ns/nm, respectively. Since the rotation of LA1 may involve structural changes in other parts in the MH, the drag coefficient for LA1 was larger than that for LA2.

In Equation (3), γ_Z was the drag coefficient per length change of a single thin filament of the sarcomere and K_Z was the spring constant for each thin filament of the sarcomere. Equation (3) follows from the fact that the sarcomeric contractile tension is just the sum of the tensile forces of the rods for all of the attached myosin molecules. The third line in Equation (1) indicates the constraint condition in the association state. This condition gives the rod strain ^tξ_i,j in relation to the conformational change of the myosin (^tx_i,j) and the sarcomeric movement (^tz).

Free Energy of a Myosin Molecule

We assume that the free energy of the myosin molecule φ in the attached state can be decomposed into the power stroke free energy φ_PS of LA1 and the deflection energy of the LA:

$\begin{array}{ll}\phi \left(x,y\right)={\phi }_{PS}\left(\theta ,y\right)+W_{\text{LA}}\left(\theta \right),\theta =y-x& (6)\end{array}$

For the deflection energy of the LA, a simple quadratic potential was assumed:

$\begin{array}{ll}{W}_{\text{LA}}\left(\theta \right)=\frac{1}{2}{K}_{\theta }{\theta }^2& (7)\end{array}$

Since there was no appropriate reference for setting the stiffness, a comparable stiffness (K_θ = 4 pN/nm) to that of the rod strain was adopted in our model. For the power-stroke free energy φ_PS, the three local minima at y = 0, s₁, and s₁ + s₂ for a fixed deflection θ = y − x are given as shown in Figure 2B, which is described by the following equations:

$\begin{array}{l}{\phi }_{PS}(\theta ,y)=\\ \{\begin{array}{ll}{E}_{Pre}+\frac{1}{2}(E_{b1}(\theta )-E_{Pre})(1-2\pi \frac{y+s_1/4}{s_1})+\frac{1}{2}{k}_Y{(y+s_1/4)}^2,\hfill & y\le -\frac{s_1}{4}\hfill \\ E_{Pre}+\frac{1}{2}(E_{b1}(\theta )-E_{Pre})(1-\cos 2\pi \frac{y}{s_1}),\hfill & -\frac{s_1}{4}<y\le \frac{s_1}{2}\hfill \\ E_{PS1}+\frac{1}{2}(E_{b1}(\theta )-E_{PS1})(1-\cos 2\pi \frac{y-s_1}{s_1}),\hfill & \frac{s_1}{2}<y\le s_1\hfill \\ E_{PS1}+\frac{1}{2}(E_{b2}(\theta )-E_{PS1})(1-\text{cos2}\pi \frac{y-s_1}{s_2}),\hfill & s_1<y\le s_1+\frac{s_2}{2}\hfill \\ E_{PS2}+\frac{1}{2}(E_{b2}(\theta )-E_{PS2})(1-\cos 2\pi \frac{y-s_1-s_2}{s_2}),\hfill & s_1+\frac{s_2}{2}<y\le s_1+\frac{5s_2}{4}\hfill \\ E_{PS2}+\frac{1}{2}(E_{b2}(\theta )-E_{PS2})(1+2\pi \frac{y-s_1-5s_2/4}{s_2})+\frac{1}{2}{k}_Y{(y-s_1-\frac{5s_2}{4})}^2,\hfill & y>\frac{5s_2}{4}\hfill \end{array}.& (8)\end{array}$

Here, E_Pre, E_PS1and E_PS2 were the three local minimum energy values at y = 0, s₁, and s₁ + s₂, respectively. These local minima correspond to the configurations of the MH and LA1 in the pre-power stroke state, and the states after the first two power strokes. The power stroke step sizes, s₁ and s₂, and the energies E_Pre − E_PS1 and E_PS1 − E_PS2 consumed in the two strokes, are given values (Table 1) similar to those used in the Monte Carlo (MC) model in our previous work (Washio et al., 2016), in which the ATP hydrolysis energy was set to E_ATP = 22K_BT following Saupe et al. (1999) at a body temperature of T = 310 K.

In Equation (8), E_b1(θ) and E_b2(θ) are the energy barriers between the minima. The heights of the energy barriers were adjusted so that enhanced beating performance was realized in the coupled simulation for the ventricle model, which is introduced below. In our model, the first barrier was assumed to be a function of the LA deflection θ as:

$\begin{array}{ll}{E}_{b1}\left(\theta \right)=\{\begin{array}{l}{E}_{b01},\theta \le {\theta }_{trap}\hfill \\ E_{b01}+C_{trap}\frac{\theta -{\theta }_{trap}}{\Delta {\theta }_{trap}},{\theta }_{trap}<\theta \le {\theta }_{trap}+\Delta {\theta }_{trap}\hfill \\ E_{b01}+C_{trap},\theta >{\theta }_{trap}+\Delta {\theta }_{trap}\hfill \end{array}& (9)\end{array}$

This first energy barrier was introduced to reproduce the stretch-activation of the cardiac muscle (Stelzer et al., 2006). In their experiment, a small, rapid stretch of ~1% of the sample length was imposed to activated skinned myocardium. Then, a nearly 10% increase in the contractile tension persisted for a time on the order of seconds compared with that of the steady state before the stretch. This suggests the existence of a trapped conformation for the MH and LA in an attached state that can be generated by the rapid stretch. By experiencing a high barrier, as in Equation (9), a myosin molecule that exhibits a large deflection θ and a large strain ξ after the first power stroke can become trapped in that state if the MH is strongly attached, since these myosin molecules cannot make progress toward a larger forward stroke, which would requires a large increment in either the deflection energy of the LA [W_LA(θ)], or the strain energy of the rod [W_rod(ξ)]. Such large LA deflections and rod strains can be generated when the thick filament was pulled rapidly to the outside. In this work, the values C_trap = 200 pN/nm, θ_trap = −0.25 nm, and Δθ_trap = 5 nm were adopted (Figure 3B) so that the appropriate response to the stretch-activation is reproduced, as shown in the numerical simulation. The second energy barrier was assumed to be a constant:

$\begin{array}{ll}{E}_{b2}(\theta )\equiv E_{b02}.& (10)\end{array}$

Control Model of Attachment and Detachment

For the transition between the attached and detached states (Figure 4), an MC model similar to the one in our previous work (Washio et al., 2016) was used. In the half-sarcomere model, the MHs were arranged on the thick filament at regular intervals, and the thin filament was divided into segments called troponin/tropomyosin (T/T) units. The transitions between the states of a T/T unit were affected by the Ca²⁺ concentration, [Ca], and by the states of the MHs below the T/T unit. In this model, only the Ca-bound state increased the affinity of the MHs for the thin filament. There are two detached state of MHs - a nonbinding state N_XB, and a weakly binding state P_XB. The affinity was adjusted by modifying the factor K_np for the rate constant of the transition from N_XB to P_XB. The relationship between the MH_i,j location and the T/T unit was determined from the offset position of the MH_i,j $({{}^{t}z+{}^t\xi }_{i,j}-{{}^{t}x}_{i.j})$ from its unloaded position (Figure 2A). A cooperative mechanism with the nearest-neighbor MHs was added by introducing the factors γ^ng and γ^−ng (γ = 40), as in our previous work (Washio et al., 2016), in which the integer ng (= 0, 1 or 2) was the number of neighboring MHs in the weakly binding state P_XB or the attached state XB. The details of the transients of the T/T unit states and between N_XB and P_XB are described in Supplementary Material S1.1.

FIGURE 4

Figure 4. The state transition Monte Carlo model of the T/T unit and the myosin molecule. The MHs in either the N_XB or P_XB states are assumed to be detached. The rate constant factors K_np and K_pn between N_XB and P_XB are affected by the state of T/T unit above it. The detachment rate constant D_PXB between P_XB and XB are given as a function of y, so that the transition to P_XB is allowed only for the MHs in the pre-power stroke position. The detachment rate constant D_NXB to N_XB is given as a function of y and ξ, so that the detachment is allowed for the second post-power stroke state or the MHs connected to the extremely stained myosin rods. The time-step strategy for reducing the computational loads is shown by the arrows. The molecular variables x, y, and ξ of the bound myosin molecules are updated by the finest basic time step Δt (black arrow), while the rod strain ξ in the detached states is updated by its multiples n_DΔt. The state transitions of the MHs and the T/T units are calculated by the MC method with the time step n_DAΔt (red arrows).

Attachment was possible only from state P_XB with the rate constant A_Pre. Detachment from the attached state XB to the weakly bound state P_XB was allowed only from the pre-power stroke state, as follows:

$\begin{array}{ll}{D}_{PXB}\left(y\right)=\{\begin{array}{ll}{D}_{PXB,Pre},\hfill & y\le s_1/4\hfill \\ \left(1-\omega \right)D_{PXB,Pre},\hfill & y=\left(1/4+\omega /2\right)s_1:0<\omega \le 1\hfill \\ 0,\hfill & y>3s_1/4\hfill \end{array}& (11)\end{array}$

Here, the variable ω could take values between 0 and 1, and was introduced to interpolate the rate constant between the pre-power stroke state and the state after the first power stroke.

In this transition, no ATP molecules were consumed, whereas detachment to N_XB required one ATP molecule. This rate constant is given as a function of both the rod strain ξ and the power stroke displacement y:

$\begin{array}{l}{D}_{NXB}\left(\xi ,y\right)=\\ \{\begin{array}{ll}\max \left(0,D_{strain}\left(\xi ,y\right)\right),\hfill & y\le s_1+s_2/4\hfill \\ \max \left(\omega D_{NXB0},D_{strain}\left(\xi ,y\right)\right),\hfill & y=s_1+\left(1/4+\omega /2\right)\hfill \\ \hfill & s_2:0<\omega \le 1\hfill \\ \max \left(D_{NXB0},D_{strain}\left(\xi ,y\right)\right),\hfill & y>s_1+3s_2/4\hfill \end{array}& (12)\end{array}$

Similar to before, the variable ω could take values between 0 and 1, and interpolated the rate constant between the states after the first and second power strokes, while D_strain indicates the forced detachment due to the extreme strain of the myosin rod:

$\begin{array}{l}{D}_{strain}\left(\xi ,y\right)=\\ \{\begin{array}{ll}{c}_{min}\left(\exp \left(a_{min}{\left(\xi -d_{min}\right)}^2\right)-1\right),\hfill & \xi \le d_{min}\hfill \\ 0,\hfill & d_{min}<\xi \le d_{max}\left(y\right)\hfill \\ c_{max}\left(\exp \left(a_{max}{\left(\xi -d_{max}\left(y\right)\right)}^2\right)-1\right),\hfill & \xi >d_{max}\left(y\right)\hfill \end{array}& (13)\end{array}$

Here, the negative strain threshold d_min was a constant, and the positive strain threshold d_max depended on the stroke displacement y as

$\begin{array}{ll}{d}_{\max }\left(y\right)=\{\begin{array}{ll}{d}_{\text{max},\text{Pre}},\hfill & y\le s_1/4\hfill \\ \left(1-\omega \right)d_{\text{max},\text{Pre}}+\omega d_{\text{max},\text{PS1}},\hfill & y=\left(14+\omega /2\right)\hfill \\ \hfill & s_1:0<\omega \le 1\hfill \\ d_{\text{max},\text{PS1}},\hfill & s{3}_1/4<y\le s_1+s_2/4\hfill \\ \left(1-\omega \right)d_{\text{max},\text{PS1}}+\omega d_{\text{max},\text{PS2}},\hfill & y=s_1+\left(1/4+\omega /2\right)\hfill \\ \hfill & s_2:0<\omega \le 1\hfill \\ d_{\text{max},\text{PS2}},\hfill & y>s_1+3s_2/4\hfill \end{array}& (14)\end{array}$

In this study, the parameters d_max,Pre = 5 nm, d_max,PS1 = 9 nm, and d_max,PS2 = 9 nm were used. These values were adjusted so that the appropriate responses to stretch-activation were reproduced. These choices did not conflict with the fact that the binding affinity to the thin filament increased as the power stroke proceeded (Llinas et al., 2015).

Multiple Time Step (MTS) Method

First, we consider a multiple time step (MTS) approach for a single half-sarcomere model (Figure 2) in which different time step intervals Δt and ΔT were adopted, respectively, for updating the molecular variables x_i,j, y_i,j, ξ_i,j and the sarcomeric shortening displacement z, when solving Equations (1, 2) coupled with Equation (3). Below, this approach will be extended to coupling with a macroscopic finite element continuum model, in which a single sarcomere model was imbedded into each finite element.

The time step ΔT was assumed to be an integer multiple of the time step interval for the molecular variables Δt:

$\begin{array}{ll}\Delta T=n·\Delta t& (15)\end{array}$

Such approaches reduce the computational overhead of the shared-memory synchronization, as well as the data communication needed in distributed parallel systems, if a sufficiently large integer n can be applied. For our Langevin dynamics model, the microscale time step t was set at 0.25 ns. This choice was constrained by the relationships between the magnitudes of the drag coefficients γ_X, γ_Y with the curvature of the potential φ. For example, in the case of a simple Langevin equation:

$\begin{array}{ll}{\gamma }_{}^t\dot{u}+\frac{d\phi }{du}({}^{t}u)-{}^{t}R=0& (16)\end{array}$

with a given free energy potential φ, a variable u, and the random force that satisfies

$\begin{array}{ll}\{\begin{array}{l}\langle {}^{t}R\rangle =0\hfill \\ \langle {}^{t}R{}^{t^{\prime }}R\rangle ={\sqrt{2\gamma k_{B}T}}^{t-t\text{'}}\delta \hfill \end{array}& (17)\end{array}$

The stability of the explicit numerical integration scheme required that

$\begin{array}{ll}\Delta t\le \frac{\gamma }{K_{max}}& (18)\end{array}$

where K_max was the maximum magnitude of the curvature of φ (||d²φ/du²||) over the range of u. Even if an implicit time integration scheme was applied, Equation (18) must be satisfied for the maximum magnitude value of the negative curvature (d²φ/du² < 0). In our case, as shown in Figure 3B, negative curvatures were unavoidable on the ridge lines of the potential landscape. For example, if Δt = 0.25 ns was used when γ = 50 pN · ns/nm, the allowable maximal curvature from Equation (18) was K_max = γ/Δt = 200 pN/nm. This curvature value implies an energy change of $K_{max}\Delta u^2=100\text{ pN}·\text{nm}$ for a displacement Δu = 1 nm. Actually, values for the magnitude of the curvature were observed near the high energy barrier between the pre-power stroke state and the state after the first power stroke in our model (Figure 3B).

Another limitation on practicable time step size comes from considerations of fluctuations Δu during each time interval Δt. If we ignore the potential φ in Equation (16), the standard deviation of Δu given by a series of random forces in Equation (17) during time Δt is

$\begin{array}{ll}{}^{\Delta t}\sigma =\sqrt{\langle \Delta u^2\rangle }=\sqrt{2k_{B}T\Delta t/\gamma }& (19)\end{array}$

At body temperature, we have k_BT = 4.278 pN · nm. Thus, for the case of γ = 50 pN · ms/nm and Δt = 0.25 ns, we have ^Δtσ ~ 0.2 nm. These displacements are large enough to make a noticeable difference in the landscape of the potential φ.

Compared with the dynamics of the molecules, the sarcomeric movement in cardiac muscle is generally much slower, as shown by the following argument. The shortening velocity of the sarcomere model is related to the stretch rate $\dot{\lambda }$ of the cardiac muscle along the fiber direction by

$\begin{array}{ll}ż=-\frac{1}{2}S{L}_0\dot{\lambda }& (20)\end{array}$

Here, SL₀ = 2.1 μm is the unloaded sarcomere length. If we assume the maximal shortening velocity of the cardiac muscle ${(-\dot{\lambda })}_{max}=5ML/s$, where ML is the muscle length (Edman et al., 1974), the maximal shortening velocity of a half-sarcomere is ${ż}_{max}=5.25\mu \text{m}/s=5.25\times 10^{-6}\text{nm}/ns$. However, the previous consideration regarding the fluctuations during the time interval Δt = 0.25 ns gives the average magnitude of the molecular velocity to be ^Δtσ/t ≈ 0.8 nm/ns. This comparison between the sarcomeric and molecular velocities suggests the possibility of applying a multi-valued time step approach, in which tens of thousands of fine time steps of size Δt are calculated when integrating the molecular variables x_i,j, y_i,j, ξ_i,j for each large one step interval ΔT used for integrating the sarcomeric variable z.

During the time integration process using the small time step Δt over the time interval [T : T + ΔT], the LA conformation variables x_i,j, y_i,j are updated explicitly, and then the rod strains ξ_i,j are temporarily updated so that the constraint in Equation (1) is fulfilled, using the most recently calculated shortening velocity ż from time T. The temporarily updated variables are denoted with bars over them, such as ${\overline{\xi }}_{i,j}$ and $\overline{z}$. When the process switches to an implicit computation of the sarcomeric shortening displacement z and its time derivative ż at time T + ΔT for use in Equation (3), the tensile forces exerted by the attached MHs during the time interval [T : T + ΔT] are computed by using the corrected rod strain ξ_i,j, for which the shortening velocity ż over the time interval [T : T + ΔT] is replaced with ^T+ΔTż. By doing so, the stiffness due to the strained rods of the attached MHs is involved in the implicit time integration of Equation (3). This implicit strategy allows us to apply a time interval ΔT which is four orders of magnitude larger than Δt.

The molecular variable time integrations can be performed using the temporal sarcomeric shortening displacement $\overline{z}$ on the time interval [T : T + ΔT] given by

$\begin{array}{ll}{}^{T+\Delta T}\overline{z}={}^{T}z+k\Delta t_{}^T\dot{z},k=1,\cdots ,n& (21)\end{array}$

The LA conformation variables for the attached MHs at time t + Δt are explicitly updated from those at time t, so that the following equations are satisfied:

$\begin{array}{ll}\{\begin{array}{l}\begin{array}{l}{\gamma }_X\frac{{{}^{t+t}x}_{i.j}-{{}^{t}x}_{i.j}}{\Delta t}+\frac{\partial \phi }{\partial x}\left({{}^{t}x}_{i.j},{{}^{t}y}_{i.j}\right)+\frac{d{W}_{rod}}{d\xi }\left({{}^t\overline{\xi }}_{i,j}\right)\\ -{{}^{t}R}_{X,i,j}=0,{{}^t\delta }_{A,i,j}=1\end{array}\hfill \\ {\gamma }_Y\frac{{{}^{t+t}y}_{i.j}-{{}^{t}y}_{i.j}}{\Delta t}+\frac{\partial \phi }{\partial y}\left({{}^{t}x}_{i.j},{{}^{t}y}_{i.j}\right)-{{}^{t}R}_{Y,i,j}=0\hfill \end{array}& (22)\end{array}$

Then, the temporal rod strains $\left\{{\overline{\xi }}_{i,j}\right\}$ at time t + Δ are updated according to

$\{\begin{array}{l}{\gamma }_X\frac{{{}^{t+t}\overline{\xi }}_{i.j}-{{}^t\overline{\xi }}_{i.j}}{t}+\frac{d{W}_{rod}}{d\xi }\left({{}^t\overline{\xi }}_{i.j}\right)-{{}^{t}R}_{D,i,j}=0,{{}^t\delta }_{A,i,j}=0\hfill \\ \begin{array}{l}{{}^{t+t}\overline{\xi }}_{i.j}-{{}^{t_{A,i,j}}\overline{\xi }}_{i.j}-\left({{}^{t+\Delta t}x}_{i.j}-{{}^{t_{A,i,j}}x}_{i.j}\right){{}^t\delta }_{A,i,j}=1\\ +{}^{t+\Delta t}\overline{z}-{}^{t_{A,i,j}}\overline{z}=0,\end{array}\hfill & (23)\end{array}$

After performing the above time integrations for k = 1, ⋯ , n over the interval [T : T + ΔT], the true sarcomeric shortening displacement z is implicitly computed by solving the following equations:

$\begin{array}{ll}\{\begin{array}{l}{{\gamma }_Z}^{T+\Delta T}\dot{z}+{K_Z}^{T+\Delta T}z{-}_{}^{T,\Delta T}F=0\hfill \\ {}^{T+\Delta T}z={}^{T}z+\Delta T{}^{T+\Delta T}\dot{z}\hfill \end{array}& (24)\end{array}$

In Equation (24), the mean total tensile force ^T,ΔTF over the time interval [T : T + ΔT] is found by applying the true rod strains {ξ_i,j} over the time interval [T : T + ΔT] according to

$\begin{array}{ll}\begin{array}{l}{}^{T,\Delta T}F={}^{T,\Delta T}\overline{F}+\frac{1}{n·n_F}{\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^{n_F}{\displaystyle \sum _{i=1}^{n_M}{{}^{T+k\Delta t}\delta }_{A,i,j}}}}\frac{d^2{W}_{rod}}{d{\xi }^2}\\ \left({{}^{T+k\Delta t}\overline{\xi }}_{i,j}\right)\left({{}^{T+k\Delta t}\xi }_{i,j}-{{}^{T+k\Delta t}\overline{\xi }}_{i,j}\right)\end{array}& (25)\end{array}$

where the temporary total tensile force is evaluated using

$\begin{array}{ll}{}^{T,\Delta T}\overline{F}=\frac{1}{n·n_F}{\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^{n_F}{\displaystyle \sum _{i=1}^{n_M}{{}^{T+kt}\delta }_{A,i,j}}}}\frac{d{W}_{rod}}{d\xi }\left({{}^{T+k\Delta t}\overline{\xi }}_{i,j}\right)& (26)\end{array}$

from the temporary rod strain values $\left\{{{}^{T+k\Delta t}\overline{\xi }}_{i,j}\right\}$. Note that Equation (25) is a linear approximation of the tensile force for the true rod strains about the temporary rod strains, for which the differences are given by

$\begin{array}{ll}{{}^{T+k\Delta t}\xi }_{i,j}-{{}^{T+k\Delta t}\overline{\xi }}_{i,j}=-\left(k-k_{A,i,j}\right)\Delta t\left({}^{T+\Delta T}\dot{z}-{}^T\dot{z}\right),k=1,\dots ,n& (27)\end{array}$

where k_A,i,j is the most recent microscale step index for k for which MH_i,j is attached. This number is initialized to zero before starting the small time steps with k = 1. By substituting Equation (27) into Equation (25), the mean total tensile force can be rewritten as

$\begin{array}{ll}{}^{T,\Delta T}F={}^{T,\Delta T}\tilde{F}-\Delta T^{T,\Delta T}{K}_F{}^{T+\Delta T}\dot{z}& (28)\end{array}$

with total mean stiffness

$\begin{array}{ll}\begin{array}{l}{{}^{T,\Delta T}K}_F=\frac{1}{n^2·n_F}{\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^{n_F}{\displaystyle \sum _{i=1}^{n_M}{{}^{T+k\Delta t}\delta }_{A,i,j}}}}\left(k-k_{A,i,j}\right)\\ \times \frac{d^2{W}_{rod}}{d{\xi }^2}\left({{}^{T+k\Delta t}\overline{\xi }}_{i,j}\right)\end{array}& (29)\end{array}$

and extrapolated mean total tensile force using ^Tż

$\begin{array}{ll}{}^{T,\Delta T}\tilde{F}={}^{T,\Delta T}\overline{F}+\Delta T^{T,\Delta T}{K}_F{}^T\dot{z}& (30)\end{array}$

By substituting Equations (28–30) into Equation (24), the implicit scheme is established as follows:

$\begin{array}{ll}{({\gamma }_Z+K_Z\Delta T+\Delta T^{T,\Delta T}{K}_F)}^{T+\Delta T}\dot{z}=-(K_Z{}^{T}z-{}^{T,\Delta T}\tilde{F})& (31)\end{array}$

To see the necessity of the above implicit coupling scheme, consider the instability of the usual explicit scheme here. If an explicit scheme for the total mean tensile force is used

$\begin{array}{ll}\{\begin{array}{l}{{\gamma }_Z}^{T+\Delta T}\dot{z}+{K_Z}^{T+\Delta T}z{-}_{}^{T,\Delta T}\overline{F}=0\hfill \\ {}^{T+\Delta T}z={}^{T}z+\Delta T{}^{T+\Delta T}\dot{z}\hfill \end{array}& (32)\end{array}$

instead of Equation (24), the time step size ΔT is limited by the total mean stiffness by

$\begin{array}{ll}\Delta T<\frac{{\gamma }_Z+\Delta T{K}_Z}{{{}^{T,\Delta T}K}_F}& (33)\end{array}$

As an illustration, in the case of ${\gamma }_Z=10^4\text{ pN}·\text{ns}/\text{nm}$, as assumed in our previous work (Washio et al., 2017), ^T,ΔTK_F = 28 pN/nm, 20 attached MHs, the stiffness of each rod set to 2.8 pN/nm, and K_Z ≈ 0, the constraint in Equation (33) would be ΔT < 360 ns. However, the proposed algorithm is stable for any time step size, as far as the linear approximation in Equation (28) is concerned.

In coupling with the macroscopic finite element model, a half-sarcomere model is assigned to each element, for which Equation (20) is applied based on the relationship between the stretching along the fiber orientation f and the deformation gradient tensor:

$\begin{array}{ll}{}^T\lambda =\Vert \frac{{\partial }_{}^T\text{x}}{\partial \text{X}}\text{f}\Vert & (34)\end{array}$

Here, ^Tx = ^Tx(X) is the current position at time T of the material point X in the unloaded condition. Specifically, the following equation, obtained from Equation (34), is substituted into Equation (20).

$\begin{array}{ll}{}^T\dot{\lambda }=\frac{1}{{}^T\lambda }\left(\frac{{\partial }_{}^T\dot{\text{x}}}{\partial \text{X}}\text{f}\right)·\left(\frac{{\partial }_{}^T\text{x}}{\partial \text{X}}\text{f}\right)& (35)\end{array}$

From Equations (20, 28), the mean total tensile force of each thin filament is given by

$\begin{array}{ll}{}^{T,\Delta T}F={}^{T,\Delta T}\tilde{F}-\Delta T\frac{{{}^{T,\Delta T}K}_F}{2}S{L}_0{}^{T+\Delta T}\dot{\lambda }& (36)\end{array}$

Here, ^Tż in Equation (30) is also replaced with −SL₀^T$\dot{\lambda }/2$ to determine ${}^{T,\Delta T}\tilde{F}$. Thus, the total active tension per unit area in the unloaded configuration, the nominal stress, is given by

$\begin{array}{ll}{{}^{T,\Delta T}T}_f=2\frac{R_s}{S{A}_0}{}^{T,\Delta T}F=2\frac{R_s}{S{A}_0}\left({}^{T,\Delta T}\tilde{F}-\Delta T\frac{{{}^{T,\Delta T}K}_F}{2}S{L}_0{}^{T+\Delta T}\dot{\lambda }\right)& (37)\end{array}$

Here, SA₀ is the cross-sectional area of a single thin filament and R_S denotes the volume ratio of the sarcomere. The factor of two in Equation (37) comes from the fact that ${}^{T,\Delta T}F$ is the total tensile force given by the MHs surrounding one of the double spirals along the thin filament.

Although a small time step on the order of 1 ns must be used for the time integration of the molecular variables, a larger time step can be applied to the MC state-transition phase. Thus, it is reasonable to apply a much larger time step size, as long as it is an integer multiple of Δt, to the computation of the MC state-transitions. Furthermore, even for the time integration of the molecular variables, a coarser time step than the one used for the attached MHs can be applied to the detached MHs, since the magnitudes of the curvatures are different for the potentials W_rod and φ (Figure 3).

Coupling With the Finite Element Ventricle Model

In the beating-ventricle simulation, the Ca²⁺ transient is given for each element of the ventricle model (Figure 5). By referencing the Ca²⁺ transients, together with the stretching λ and the stretching rate $\dot{\lambda }$ along the fiber direction, the molecular variables were integrated using the small time step Δt, while the macroscopic displacements of the continuum were computed using the large time step ΔT. As derived in the Supplementary Material S3, the active stress on the continuum at time T + ΔT is represented by the first Piola–Kirchhoff stress tensor:

$\begin{array}{ll}{\Pi }_{\text{act}}=\frac{{{}^{T,\Delta T}T}_{\text{f}}}{T+\Delta T_{\lambda }}\text{f}\otimes \text{f}·{\left(\frac{{\partial }_{}^{T+\Delta T}\text{x}}{\partial \text{X}}\right)}^T& (38)\end{array}$

In the definition of the tension ${{}^{T,\Delta T}T}_f$ in Equation (37), the stiffness due to the attached MHs is implicitly included by the use of ${}^{T+\Delta T}\dot{\lambda }$ for ${}^{T+\Delta T}\dot{z}$ in Equation (28). See also the explanation of the stiffness in the Supplementary Material S3. Thus, the proposed scheme is stable for any size of time step.

FIGURE 5

Figure 5. Electrical analog circuits connected to the valve interfaces. The systemic circulation model represents blood flow from the left to the right ventricle through the aortic (AO) and the tricuspid (TR) valves, while the pulmonary circulation model represents blood flow from the right to the left ventricle through the pulmonary (PU) and the mitral (MI) valves. The boxed inset shows the Ca²⁺ transient profile with the delay time T_D.

The governing equation in the macroscale to be solved can be represented by

$\begin{array}{l}{\displaystyle {\int }_{\Omega }^{}}\delta \dot{u}·\rho \ddot{u}d\Omega +{\displaystyle {\int }_{\Omega }^{}}\delta \dot{\text{Z}}:{\left(\Pi +2pJ{F}^{-1}\right)}^{T}d\Omega \\ =P_L{\displaystyle {\int }_{{\Gamma }_L}^{}}\delta \dot{u}·nd{\Gamma }_L+P_R{\displaystyle {\int }_{{\Gamma }_L}^{}}\delta \dot{u}·nd{\Gamma }_R& (39)\end{array}$

$\begin{array}{ll}{\displaystyle {\int }_{\Omega }^{}}\delta p\left(2(J-1)-\frac{p}{\kappa }\right)d\Omega =0& (40)\end{array}$

Here, u = ^Tu(X)=^Tx(X)−X is the displacement of the material at point X ∈ Ω at time T, ρ is the density of the heart muscle, F = ∂x/∂X is the deformation gradient tensor, Z = ∂u/∂X is the displacement gradient tensor, J = det F is the Jacobian, p is the hydrostatic pressure, κ is the bulk modulus, and P_L and P_R are the blood pressures in the left and right ventricles, respectively. Ω is the muscle domain in the reference configuration, while Γ_L and Γ_R are the blood–muscle interfaces of the left and right ventricles, respectively, in the current configuration at time T, and n is the normal unit vector directed from the cavity to the muscle at these surfaces (Figure 5). The Dirichlet boundary condition ^Tu(X)=0 is imposed on the boundary nodes around the valve rings. The first Piola–Kirchhoff stress tensor Π consists of the active, passive, and viscous stresses:

$\begin{array}{ll}\Pi ={\Pi }_{act}+{\Pi }_{pas}+{\Pi }_{vis}& (41)\end{array}$

where Π_act is given by Equation (38), and the others are, respectively, the passive and viscous stresses, as described in our previous work (Washio et al., 2016). The details of these two stress tensors are given in the Supplementary Material S4.

The ventricle blood pressures P_L and P_R were determined through their interactions with the circulatory system of the body. These were modeled as electrical analog circuits, using the same parameters described in our previous work (Washio et al., 2016). The details of the circuit model that includes the atrial model are given in the Supplementary Material S5. In particular, the flow rates at the inlets and the outlets were associated with the rates of volume change in the cavity according to:

$\begin{array}{ll}\{\begin{array}{l}{\displaystyle {\int }_{{\Gamma }_L}\dot{\text{u}}\cdot \text{n}d{\Gamma }_L=F_{MI}-F_{AO}}\hfill \\ {\displaystyle {\int }_{{\Gamma }_R}\dot{\text{u}}\cdot \text{n}d{\Gamma }_R=F_{TR}-F_{PA}}\hfill \end{array}& (42)\end{array}$

Here, F_MI, F_AO, F_TR, and F_PA were the flow rates, respectively, through the mitral, aortic, tricuspid, and pulmonary valves (Figure 5). These flow rates were determined by Ohm's law while taking the rectification of the valve into account.

$\begin{array}{ll}F=H(\overline{F})\overline{F}& (43)\end{array}$

Here, $\overline{F}$ was the flow rate in the case of no rectification, and H was the relaxed Heaviside function:

$\begin{array}{ll}H\left(\overline{F}\right)=\{\begin{array}{ll}0,\hfill & \overline{F}<0\hfill \\ {\left(\overline{F}/{\overline{F}}_0\right)}^2\left(3-2\overline{F}/{\overline{F}}_0\right),\hfill & 0\le \overline{F}\le {\overline{F}}_0\hfill \\ 1,\hfill & \overline{F}>{\overline{F}}_0\hfill \end{array}& (44)\end{array}$

In our simulation, the value ${\overline{F}}_0=5\text{ mL}/\text{s}$ was used.

The macroscopic variables, including the acceleration $\ddot{u}$, velocity $\dot{u}$, and displacement u at time T + ΔT were found using Newton–Raphson iteration until the equilibrium condition was satisfied with the Newmark-beta time integration scheme (Supplementary Material S6). During the iterations, the active stress in Equation (38) was redefined with Equation (37), in which the microscopic computational results ${}^{T,\Delta T}\tilde{F}$ and ${{}^{T,\Delta T}K}_F$ were reused. Thus, switching between computations at the two scales only happened once for each macroscopic time step.

Results

Computer System

To perform the simulations, a distributed parallel system was used. Each node consisted of two Intel^® Xeon^® E5-2670 (20 MB Cache, 2.6 GHz) processors, and each processor was composed of 8 cores. In the single sarcomere simulations, only one node was used for shared memory OpenMP parallelization. In the beating-ventricle simulations, the elements of the ventricle wall were equally distributed to the nodes, while the macroscopic computations were performed only at the master node. In the microscopic computations, the time integrations of the molecular variables were parallelized using OpenMP by dividing the filaments equally among the 16 cores.

Validation of the MTS Scheme via Single Sarcomere Oscillation

The accuracy and computational efficiency of the MTS scheme were validated by numerical experiments with a single half-sarcomere model, in which 48 thin filaments were connected to a common Z-line (Figure 2A). In our previous work (Washio et al., 2017), we showed that the spontaneous oscillatory behavior of the sarcomere (Ishiwata et al., 2011) can be explained by the power stroke principle after applying a simple ordinary differential equation model. In this case, the collective reversal power strokes induced quick sarcomeric lengthening. Here, we show that this could also be reproduced by the Langevin dynamics model, regardless of the choice of macroscale time step size in the MTS scheme. In this numerical experiment, the spring constant K_Z was set to 1 pN/nm per thin filament, and the viscosity coefficient γ_Z was set to 10⁴ pN · ns/nm per thin filament. During the simulations, the Ca²⁺concentration was kept at the constant value of 1 μM.

In Figure 6, the shortening displacements obtained by using a conventional single-scale integration scheme (Δt = ΔT = 0.25 ns) and the MTS scheme (Δt = 0.25 ns, ΔT = 5, 000 ns) are compared for both the no-trap and trap models. In the no-trap model, the dependence of the first energy barrier height E_b1(θ) on the LA deflection θ in Equation (9) was eliminated, and the baseline of the energy barrier E_b01 was higher when compared with the one in the trap model (Table 1), so that a similar maximal tensile force is obtained in both models. Next, the state-transitions were computed with Δt = 0.25 ns. In these numerical experiments, the simulations started from an initial state in which all of the MHs were in N_XB, and an identical series of random forces and pseudorandom numbers for the MC state-transitions were applied to all the simulations. In case of the no-trap model (Figure 6A), similar amplitudes and periods were obtained for the shortening displacements, although there were deviations in the timing of the sharp declines. In case of the trap model (Figure 6B), the large dips in the displacements disappeared. Instead, rapid small vibrations appeared. In this case, similar initial rises, periods, and amplitudes of vibrations were obtained for the both time step sizes of ΔT.

FIGURE 6

Figure 6. The shortening displacements obtained by the standard scheme (blue: Δt = ΔT = 0.25 ns) and the MTS scheme (red: Δt = 0.25 ns, ΔT = 5, 000 ns) for the spontaneous oscillations of the single half-sarcomere model with n_M = 38 and n_F = 48. (A) The comparison for the no-trap model. (B) The comparison for the trap model.

As depicted in Figure 4, the attached MHs in the XB state were classified according to their power stroke displacement y, as follows:

$\begin{array}{ll}\{\begin{array}{l}\text{Pre}=\{{\text{MH}}_{i,j}\in \text{XB}:y_{i,j}<s_1/2\}\hfill \\ {\text{PS}}_1=\{{\text{MH}}_{i,j}\in \text{XB}:s_1/2\le y_{i,j}<s_1+s_2/2\}\\ {\text{PS}}_2=\{{\text{MH}}_{i,j}\in \text{XB}:y_{i,j}\ge s_1+s_2/2\}\hfill \end{array}& (45)\end{array}$

These states can be regarded as the pre-power stroke, the state after the first power stroke, and the state after the second power stroke, respectively. As suggested by our previous work (Washio et al., 2017), a large pulsed flux of the reversal power strokes from PS₂ to Pre over PS₁ generated the sharp decline in z for the no-trap model (Figure 7A). In the trap model, this reversal flux was trapped at PS₁, so that the decline in z was stopped at small changes, leading to Δz > −10 nm (Figure 7B), which corresponds to the stroke size of the LA.

FIGURE 7

Figure 7. The temporary change of the state ratios classified to the three power stroke stages (Pre, PS₁, and PS₂) obtained by the standard scheme (ΔT = 0.25 ns) for the spontaneous oscillations of the single half-sarcomere with the no-trap model (A) and the trap model (B). In case of the no-trap model (A), a large pulsed flux of the reversal power strokes from PS₂ to Pre through PS₁ generates the sharp decline of z around T = 370 ms. In case of the trap model (B), the flux of the reversal power strokes is trapped at PS₁, so that the decline of z is stopped within small changes Δz > −10 nm.

To test the stability of the MTS scheme, simulations using the explicit scheme given by Equation (32) were performed with a much smaller time step of ΔT = 500 ns (Figure 8). Although the explicit scheme also yielded good results at first, the computational results became totally invalidated when the active stiffness ${{}^{T,\Delta T}K}_F$ exceeded the threshold indicated by Equation (33), as estimated previously. Furthermore, oscillatory behavior could not be reproduced with the explicit scheme. This result suggests the drawback of explicitly using the active tensions, which occurs when solving a system of ordinary differential equations with a finer time step in coupled simulations. As shown in Figure 8A, the calculated force using the explicit scheme did not diverge, although the oscillatory behavior was completely lost. Thus, it is difficult to judge the accuracy of numerical results by examining only one case. As shown here, we must compare the results of different macroscale time step sizes ΔT to confirm the accuracy of the coupling scheme.

FIGURE 8

Figure 8. The shortening displacements obtained by the standard scheme (blue: ΔT = 0.25 ns) and the explicit scheme (red: ΔT = 500 ns) for the spontaneous oscillations of the single half-sarcomere model with n_M = 38 and n_F = 48. (A) The comparison for the no-trap model. (B) The comparison for the trap model.

The above simulations were executed on one node consisting of 16 cores using shared memory in OpenMP parallelization. Thus, in the parallelization, three filaments were assigned to each core. The averaged elapsed times for the 1-ms time integration were 125 and 97 s, with the standard integration scheme (Δt = ΔT = 0.25 ns) and the MTS scheme (Δt = 0.25 ns, ΔT = 5, 000 ns), respectively. The difference in the elapsed times came from the machine synchronization overhead, and the differences in the computational loads for the various filaments. With the MTS scheme that lumps 20,000 steps, the differences in computational loads between the filaments during each small time step were tremendously diminished. For a single-sarcomere simulation, using a much smaller time step size for ΔT was sufficient to attain good parallel efficiency because the overhead associated with updating z was negligible. However, a large step size was necessary when the sarcomere model was coupled with the macroscopic ventricle model because the communication overhead between the large number of nodes became greatly increased, along with the computation time for updating the macroscopic variables.

Validation of Basic Sarcomere Properties

The basic properties of the actomyosin trap model, which includes the SL and [Ca] dependences of the contractile force, the isometric twitch, the responses for the isotonic contraction, and the quick shortening of the half-sarcomere, along with the details of these numerical experiments, are presented in the Supplementary Material S2. The results of these numerical experiments confirm the validity of our half-sarcomere model. Here, the force-velocity curve obtained at a constant Ca²⁺ concentration ([Ca] = 1 μM) is examined in context with the behavior of the bound myosin molecules during the isotonic contractions at the various shortening velocities (Figure 9). As the shortening velocity increased, the state ratio of PS₂ increased (Figure 9B), because the joint point P was pushed forward (y increased) more strongly by the deflection potential W_LA(θ) in Equation (7) with the larger negative deflection θ = y − x (Figure 9D). Note that the negative averaged rod strain ξ at PS₂ for a shortening velocity larger than 1 μm/s (Figure 9C) does not imply a negative contractile force, because dW_rod/dξ(ξ) ≫ −dW_rod/dξ(−ξ) for any positive strain ξ>0, except for ξ~0 as shown in Figure 3A.

FIGURE 9

Figure 9. The force-velocity relation for [Ca] = 1 μM with the shortening velocity of the half-sarcomere in the horizontal axis (A). The state ratios of the pre-power stroke (Pre) and the post-second power stroke state (PS₂) are plotted with respect to the shortening velocities (B). Similarly, the averaged rod strains ξ in the states Pre and PS₂ (C), the averaged molecular variables x and y of the molecular deformation (F) in the states Pre and PS₂ are also plotted, respectively, in (D,E).

Stretch-Activation by Trapped Myosins

To see the effectiveness of the trapping mechanism in the state after the first power stroke PS₁ created by the energy barrier in Equation (9), together with the zero detachment rates for PS₁ in Equations (11, 12), a stretch-activation test was performed for the single half-sarcomere model consisting of 48 filament pairs (Figure 10). Here, a 1% stretch was applied over the 1-ms time interval starting at T = 150 ms, at which time the contractile force had sufficiently matured. In the simulation, the time step sizes were set at Δt = 0.25 ns and ΔT = 25 ns. The state-transitions were also computed using Δt = 0.25 ns. During the simulations, the Ca²⁺ concentration ([Ca]) was kept at the constant value of 10 μM.

FIGURE 10

Figure 10. Numerical results of the stretch-activation simulation. A shift of Δz = −9 nm was applied over the time span of 1 ms at T = 150 ms. (A) The time courses of the contraction force per one thin filament (black) and the applied shortening displacement z (blue) are shown. The dot-dash line indicates the baseline of the force before the stretch. (B) The state ratios for the three power stroke stages of the attached MHs. (C) The averaged molecular variables x and y for the MHs at PS₁.

A roughly 15% increase in the contractile force lasted at least 2 s after the quick stretch (Figure 10A). This long-lasting increase in the force compared with the pre-stretch steady state was apparently due to the lasting increase in the population of PS₁ (Figure 10B: orange line). The persistent increase of the averaged LA deflection, θ = y − x, for MHs in PS₁ (Figure 10C) indicates that it was generated by the MHs trapped by the higher free energy barrier E_b1(θ) defined by Equation (9). Compared with the experimental results given in Stelzer et al. (2006), our numerical result misses “Phase 2,” in which the force drops one time to the steady state level before the stretch. However, the magnitude of the force incrementation after that agrees with the experimental facts.

Beating-Ventricle Simulations

Beating-ventricle simulations were performed using a finite element ventricle model consisting of 7,600 tetrahedral elements. In each element, a sarcomere model consisting of 8 filament pairs was imbedded along the appropriate fiber orientation f. The distribution of the fiber orientations (Figure 1) was found by an optimization algorithm (Washio et al., 2016) based on the impulses given by the active tension, which was computed using the MC crossbridge model instead of the Langevin model to reduce the heavy computational loads. Portions of the helical fiber structure are depicted in Figure 1. As confirmed in our previous work (Washio et al., 2016), this algorithm constructed a fiber distribution that was quite similar to the one obtained by diffusion tensor magnetic resonance imaging (DTMRI) measurements. The heart rate was set to 60 beats per minute, and the Ca²⁺ transient (Figure 5) generated by the mid-myocardial cell model proposed by ten Tusscher and Panfilov (2006) was applied. The transmural delays of the Ca²⁺ transient determined by the distances from the endocardial surfaces of the left and right ventricles under a transmural conduction velocity of 52 cm/s, as measured by Taggart et al. (2000), was adopted. The deformation of each element was linked to the sarcomeric shortening displacement using Equations (34, 35). In the simulations, the optimized time step algorithm represented in Figure 4 was applied. Essentially, the values Δt = 0.25 ns and T = 5, 000 ns were used, so that n = 20, 000. However, the state-transitions were computed every 2.5 ns (n_DA = 10), and the time integration for the detached MHs were performed every 1.25 ns (n_D = 5).

In the crossbridge model, the trap and the no-trap models using the various power-stroke free energy potential functions φ_PS were used, as with the simulations of the single sarcomere oscillation (Table 1). By comparing it with the no-trap model in Figure 11A, the trap mechanism can be seen as contributing to maintaining the high pressure in the last half of the systolic phase. As a result, the blood volume ejected from the left ventricle in the trap model increased to 77 from 68 mL, while the ATP energy consumption of the left ventricular wall decreased to 5.9 from 6.4 J (Figure 11E). This implies that the trap mechanism serves to increase the blood ejection, while also decreasing the energy consumption. Note that the ATP consumption rates were computed by counting the detachments of MHs in PS₂ to those in N_XB, which was controlled by the rate constant D_NXB defined in Equation (12).

FIGURE 11

Figure 11. Numerical results of the beating-ventricle simulation using the FEM ventricle model. In each element, the sarcomere model consisting of 8 filament pairs was imbedded. (A) The time courses of the left ventricular pressure (solid lines) and volume (broken lines) with the no-trap MH model (red) and the trap model (black). (B–D) The time courses of the population ratio of attached MHs in the left ventricular wall classified to the pre-power stroke state (B: Pre), the first post-power stroke state (C: PS₁), and the second post-power stroke state (D: PS₂). (E) The time courses of the cumulative ATP energy consumption in the left ventricular wall.

As shown in Figures 11B–D, two increases in the population of MHs in state PS₁ can be seen; one at the beginning of the systolic phase, and one at the final half. These increases correspond to reversals in the left ventricular pressures of the trap and the no-trap models, as shown in Figure 11A. In the systolic phase, the cardiac myocytes supported their contractile tension along the shared fiber bundle, in which the active stress in Equation (38) provided the great majority of the total stress in Equation (41). Therefore, from the mechanical equilibrium condition along a fiber bundle, the active tensions must be almost equal. If there was a delay in the provision of the active tension, or a relaxation during the intermediate systolic phase at one point of the fiber bundle, this portion quickly became lengthened, and the sarcomeres in the remaining parts shortened until reaching a mechanical equilibrium. Since this transition accompanied decreases in the active tension of the sarcomeres, stopping the process as early as possible was desirable. The trap mechanism could achieve this goal, as shown in Figure 12, in which the distributions of the population of MHs in states PS₁ and PS₂ at the end of the systolic phase (T = 0.25 s) were compared. As shown in Figure 12B for the trap model, the higher populations in the PS₁ state were seen in the regions where the populations in state PS₂ were lower than in the other regions. This indicates that the decrease in the population of MHs at PS₂ was sufficiently compensated for by the trapped MHs in state PS₁. However, although the population in PS₂ for the no-trap model was similar to one of the trap model, the active tension was nearly half that of the trap model for the entire region (Figure 12A). In particular, the active tensions with the no-trap model were much smaller than those with the trap model, even in the regions with large PS₂ populations. This indicates the importance of maintaining the active tension along the fiber bundle. The distributions of the active tension values and the state populations over the entire cycle are shown in Supplementary Video 1.

FIGURE 12

Figure 12. The active tension (left), and the population of MHs in the first post-stroke state PS₁ (center) and the second post-stroke state PS₂ (right) at 0.25 s for the no-trap model (A) and the trap model (B) in the three cross sections perpendicular to the longitudinal axis.

The importance of the trap for synchronizing contraction and relaxation over the entire ventricle is further confirmed by Figure 13, in which the behaviors of the sarcomere model with the no-trap and the trap model imbedded with identical elements at the apical septal segment are compared. With the no-trap model (Figure 13A), there was a prominent decline in the sarcomere shortening displacement z that accompanied the large drops in the active tension around T = 0.18 s. This drop in the active tension was caused by shifts in the population of MHs from PS₂ to the pre-power stroke state Pre, as indicated in Figure 13C. As shown previously in the simulations of sarcomere oscillation, each sarcomere had the ability to undergo quick lengthening after a certain duration of contraction. However, the slow decline of LVP in the no-trap model (Figure 11A) at the end of the systolic phase indicates that this characteristic was not necessarily exploited for the quick relaxation of the whole ventricle before the next diastolic phase because the timing of the relaxation changed depending on the Ca²⁺ transients and the sarcomeric movements. Furthermore, a relaxation prior to a sufficient drop in the Ca²⁺-concentration was followed by the next contraction, as shown in Figure 13C, around T = 0.2 s. This contraction of the sarcomere did not efficiently contribute to increasing the ejected blood volume, as indicated by LVV in Figure 11A. However, the blood ejection lasted until T = 0.3 s in the trap model. Thus, maintaining the active tension with the trapped MHs in PS₁, which corresponded to a rise in the population of PS₁ during the time interval [0.23, 0.3] (Figure 13D), substantially contributed to the ejected blood volume.

FIGURE 13

Figure 13. The behavior of the sarcomere in the systolic phase for the no-trap and the trap models imbedded in an identical element at the apical septal segment. (A,B) The active tension (red) and the sarcomeric shortening Z (blue). (C,D) The population ratio of attached MHs.

Figure 14 compares the distributions of the attached MHs, which are imbedded in 33 elements at the apical septal segment, in the (y, θ) coordinate at T = 0.1, 0.2, and 0.3 s of the no-trap and trap models. Although the distributions at the beginning of the systolic phase (T = 0.1 s) were nearly the same for both models, differences were found in regions of higher deflection θ in the Pre and PS₁ states at the peak of the systolic phase (T = 0.2 s), and in PS₁ and PS₂ at the end of the systolic phase (T = 0.3 s). Note that the large deflection (θ > 0) of the LA created high strain (ξ > 0) in the rod due to the equilibrium condition for the variable x in Equation (1). However, these MHs in the Pre state of the no-trap model disappeared quickly due to their large rate of detachment into state P_XB in Equation (11) and Table 1 ($D_{PXB,Pre}=3,000{\text{s}}^{-1}$), so that they did not contribute to maintaining the active tension. However, the MHs in state PS₁ were trapped there so long as these myocytes were strongly pulled by the surrounding activated myocytes.

FIGURE 14

Figure 14. The distribution of the attached MHs at the initial time (T = 0.1 s), the peak (T = 0.2 s), and the end (T = 0.3 s) of the systolic phase for the trap model (A) and the no-trap model (B) for the (y, θ) coordinate. The contours represent the landscape of the free energy φ_PS(θ, y) + W_LA(θ). The plots are for the attached MHs of the sarcomere models imbedded into the 33 elements in the apical septal segment.

Finally, the computational load and the parallel efficiency were examined. For the microscale computations, the elements of the finite element model were equally distributed to the available cores. But, for the macroscale finite element computations, only one node consisting of 16 cores was used, and the remaining nodes were in the waiting state, since the finite element model was relatively small (7,700 elements). Thus, the parallel efficiency came from the proportion of the macroscale computational time, compared with the total computation time. With the original setup (n = 20, 000, n_DA = 10, n_D = 5, Δt = 0.25 ns), the parallel computation with 1,920 cores required 105 h per heartbeat. Within this total elapsed time, 16% was occupied by the macroscale computations. Thus, good parallel efficiency was achieved. Further evaluations of the parallel efficiency are given in the Supplementary Material S7.

Discussion

Accuracy, Stability, and Efficiency of the MTS Scheme

The MTS scheme coupled the integration of the molecular variables that use the small time step Δt with the integration of the sarcomere shortening variable z that used the coarse time step ΔT, which is a large integer multiple of Δt. Since sarcomere shortening is linked to the shortening of the continuum along the fiber orientation by Equation (35), the same coupling scheme can be applied to the coupling with the finite element model. The key point of the proposed MTS scheme is that the active tension at time T + ΔT is implicitly determined by combining the stretch rate of the continuum along the fiber orientation at T + ΔT, as given in Equation (37), in which the stiffness of the attached myosin rods during the time interval [T : T + ΔT] given by Equation (29) is used. By applying this implicit scheme, an appropriate time step interval ΔT can be chosen for the macroscale computation to diminish the synchronization and communication overhead in the distributed memory parallel system. The accuracy of the MTS scheme, in which the time step ratio was set to 0.25 ns: 5 μs, was validated using a simulation of the spontaneous oscillation of a single sarcomere, and by comparing the numerical results with those computed using equal time intervals.

Required Computational Power for the Coupled Simulation

For the beating-ventricle simulation of ventricle model consisting of 7,600 elements, 105 h were required for each beat using 1,920 cores and a 0.25-ns time step integration in the molecular computations, and a 5-μs time step integration for the macroscopic finite element computation. Within this computation, 84% of the total time was consumed by the microscopic molecular computation. In this simulation, 4 elements were assigned to each core, in which the sarcomere model consisted of 8 filament pairs imbedded in each element. Therefore, the CPU time per filament pair was ~2.8 h. This is the fastest case, not counting the macroscale computational case in which one core was assigned to each filament. Even for the rather coarse mesh model consisting of 7,600 elements, this fastest computation required 60,800 (= 7, 600 × 8) cores. This shows that our application still required huge computational power.

Potential of the Coupled Approach

In this paper, an effective utilization of the coupled approach to explore the macroscopic effects of a molecular mechanism was shown. Regarding the molecular mechanism, the power-stroke free energy potential was constructed so as to reproduce the stretch-activation for the single-sarcomere model. In this model, the energy barrier between the pre-power stroke state and the state after the first power stroke was made higher for large positive lever arm deflections, which meant that large loads were imposed on the myosin rods and heads. If the pre-power stroke state and the state after the first power stroke correspond to, respectively, the so-called “P_i-release state” and “ADP state,” the forward and reversal power stroke transitions accompany the release and the rebinding of inorganic phosphate (P_i), respectively (Llinas et al., 2015). Thus, if the larger load on the MH closes the channel in which P_i travels during the transitions, the height of the free energy barrier could increase. In the proposed numerical model, this hypothesis was reflected by the landscape of the free energy φ_PS(θ, y), as mentioned above. The coupled approach revealed that the proposed mechanism for the myosin molecule contributed to maintain the high systolic blood pressure for the appropriate period by synchronizing relaxations along the fiber bundles. Stelzer et al. (2006) discussed the possibility of stretch-activation reinforcing regions where stronger contractile tensions were required during the entire systolic phase, while our numerical results suggest that its function is to reinforce the regions that start relaxation earlier than other regions. Of course, this is still just a hypothesis linking the stretch-activation to the performance of the beating heart. However, this function of stretch-activation function at the end of the systolic phase has gone unnoticed until now.

Limitations

In the coupling approach, a single half-sarcomere model was directly imbedded into each element of the macroscopic ventricular mesh. This means that the periodically repeated pattern of single sarcomere movement was imposed along the filament direction within each element. Thus, the synchronization of the sarcomeres within each element can be assumed. In reality, relaxations of sarcomeres within the same myofibril are not necessarily synchronized. Thus, even though each of the sarcomeres was stretched quickly during relaxation, as shown in the spontaneous oscillation simulation, the stretch speed of the entire cardiac cell may be slowed due to time lags. One way to account for such an effect in the simulation model is to imbed a myofibril model, in which an adequate number of sarcomeres are connected in series, into each element. Obviously, such an approach requires even greater computational resources.

New Insights of Cardiac Muscle Relaxation in a Beating Heart

Using the numerical experiments on the single-sarcomere model, spontaneous oscillatory behavior was recovered via the Langevin dynamics model with a simple power-stroke free energy, as in Equation (8) with a constant energy barrier [E_b1(θ) ≡ E_b01]. The prominent characteristic of this oscillation is the quick lengthening induced by collective reversal strokes (Figure 7A). At first glance, it appears that this mechanism operated by quickly relaxing the muscle against the slow decline of the Ca²⁺ concentration (Inset in Figure 5). However the timing of the lengthening events differ from those in the ventricle wall due to the various feedback signals from the local muscle movements, resulting in the slow decline of the LVP (Figure 11A). Using the numerical experiments on the ventricular model, we see that the trap mechanism contributes to the synchronization of muscle relaxation by halting sarcomeric lengthening if it occurs earlier than in the neighboring muscle. We also see that the same trap mechanism causes the stretch-activation phenomenon at the tissue level.

Author Contributions

TW and TH: designed the project; TW and RK: designed and conceived the numerical model; TW and J-IO: constructed the simulation code and the input data; TW and SS: wrote the paper with input from TH.

Funding

This work was supported in part by the Ministry of Education, Culture, Sports, Science and Technology of Japan (MEXT) as Priority Issue on Post-K computer (Integrated Computational Life Science to Support Personalized and Preventive Medicine) (Project ID: hp170233). RK's work was supported in part also by the Research Complex Promotion Program.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors thank Louis R. Nemzer, Ph.D., from Edanz Group (www.edanzediting.com/ac) for editing a draft of this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.00333/full#supplementary-material

References

Cahill, T. J., Ashrafian, H., and Watkins, H. (2013). Genetic cardiomyopathies causing heart failure. Circ. Res. 113, 660–675. doi: 10.1161/CIRCRESAHA.113.300282

PubMed Abstract | CrossRef Full Text | Google Scholar

Clayton, R. H., Bernus, O., Cherry, E. M., Dierckx, H., Fenton, F. H., Mirabella, L., et al. (2011). Models of cardiac tissue electrophysiology: progress, challenges and open questions. Prog. Biophys. Mol. Biol. 104, 22–48. doi: 10.1016/j.pbiomolbio.2010.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Edman, K. A., Mattiazzi, A., and Nilsson, E. (1974). The influence of temperature on the force-velocity relationship in rabbit papillary muscle. Acta Physiol. Scand. 90, 750–756. doi: 10.1111/j.1748-1716.1974.tb05643.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Gardiner, C. W. (2004). Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, 3rd Edn, Springer Series in Synergetics, Vol. 13. Berlin: Springer-Verlag.

Grandi, E., Pasqualini, F. S., and Bers, D. M. (2010). A novel computational model of the human ventricular action potential and Ca transient. J. Mol. Cell. Cardiol. 48, 112–121. doi: 10.1016/j.yjmcc.2009.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Gurev, V., Lee, T., Arevalo, H., and Trayanova, N. A. (2011). Models of cardiac electromechanics based on individual hearts imaging data: image-based electromechanical models of the heart. Biomech. Model. Mechanobiol. 10, 295–306. doi: 10.1007/s10237-010-0235-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Howard, J. (2001). Mechanics of Motor Proteins and the Cytoskeleton. Sunderland, MA: Sinauer Associates.

Google Scholar

Huxley, A. F., and Simmons, R. M. (1971). Proposed mechanism of force generation in striated muscle. Nature 233, 533–538. doi: 10.1038/233533a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishiwata, S., Shimamoto, Y., and Fukuda, N. (2011). Contractile system of muscle as an auto-oscillator. Prog. Biophys. Mol. Biol. 105, 187–198. doi: 10.1016/j.pbiomolbio.2010.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerckhoffs, R. C., Neal, M. L., Gu, Q., Bassingthwaighte, J. B., Omens, J. H., and McCulloch, A. D. (2007). Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation. Ann. Biomed. Eng. 35, 1–18. doi: 10.1007/s10439-006-9212-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Llinas, P., Isabet, T., Song, L., Ropars, V., Zong, B., Benisty, H., et al. (2015). How actin initiates the motor activity of Myosin. Dev. Cell. 33, 401–412. doi: 10.1016/j.devcel.2015.03.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, C. H., and Rudy, Y. (1994). A dynamic model of the cardiac ventricular action potential - simulatons of ionic currents and concentration changes. Circ. Res. 74, 1071–1097. doi: 10.1161/01.RES.74.6.1071

PubMed Abstract | CrossRef Full Text | Google Scholar

Marcucci, L., Washio, T., and Yanagida, T. (2016). Including thermal fluctuations in actomyosin stable states increases the predicted force per motor and macroscopic efficiency in muscle modelling. PLoS Comput. Biol. 12:e1005083. doi: 10.1371/journal.pcbi.1005083

PubMed Abstract | CrossRef Full Text | Google Scholar

Saupe, K. W., Eberli, F. R., Ingwall, J. S., and Apstein, C. S. (1999). Hypoperfusion-induced contractile failure does not require changes in cardiac energetics. Am. J. Physiol. 276, H1715–H1723.

Google Scholar

Stelzer, J. E., Larsson, L., Fitzsimons, D. P., and Moss, R. L. (2006). Activation dependence of stretch activation in mouse skinned myocardium: implications for ventricular function. J. Gen. Physiol. 127, 95–107. doi: 10.1085/jgp.200509432

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugiura, S., Washio, T., Hatano, A., Okada, J., Watanabe, H., and Hisada, T. (2012). Multi-scale simulations of cardiac electrophysiology and mechanics using the university of Tokyo heart simulator. Prog. Biophys. Mol. Biol. 110, 380–389. doi: 10.1016/j.pbiomolbio.2012.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Taggart, P., Sutton, P. M., Opthof, T., Coronel, R., Trimlett, R., Pugsley, W., et al. (2000). Inhomogeneous transmural conduction during early ischaemia in patients with coronary artery disease. J. Mol. Cell. Cardiol. 32, 621–630. doi: 10.1006/jmcc.2000.1105

PubMed Abstract | CrossRef Full Text | Google Scholar

ten Tusscher, K. H., Noble, D., Noble, P. J., and Panfilov, A. V. (2004). A model for human ventricular tissue. Am. J. Physiol. 286, H1573–H1589. doi: 10.1152/ajpheart.00794.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

ten Tusscher, K. H., and Panfilov, A. V. (2006). Alternans and spiral breakup in a human ventricular tissue model. Am. J. Physiol. Heart Circ. Physiol. 291, H1088–H1100. doi: 10.1152/ajpheart.00109.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Washio, T., Hisada, T., Shintani, S. A., and Higuchi, H. (2017). Analysis of spontaneous oscillations for a three-state power-stroke model, physical review. Phys. Rev. E 95:022411. doi: 10.1103/PhysRevE.95.022411

CrossRef Full Text | Google Scholar

Washio, T., Yoneda, K., Okada, J., Kariya, T., Sugiura, S., and Hisada, T. (2016). Ventricular fiber optimization utilizing the branching structure. Int. J. Numer. Method Biomed. Eng. 32:e02753. doi: 10.1002/cnm.2753

PubMed Abstract | CrossRef Full Text | Google Scholar