A Multiple Step Active Stiffness Integration Scheme to Couple a Stochastic Cross-Bridge Model and Continuum Mechanics for Uses in Both Basic Research and Clinical Applications of Heart Simulation

In a multiscale simulation of a beating heart, the very large difference in the time scales between rapid stochastic conformational changes of contractile proteins and deterministic macroscopic outcomes, such as the ventricular pressure and volume, have hampered the implementation of an efficient coupling algorithm for the two scales. Furthermore, the consideration of dynamic changes of muscle stiffness caused by the cross-bridge activity of motor proteins have not been well established in continuum mechanics. To overcome these issues, we propose a multiple time step scheme called the multiple step active stiffness integration scheme (MusAsi) for the coupling of Monte Carlo (MC) multiple steps and an implicit finite element (FE) time integration step. The method focuses on the active tension stiffness matrix, where the active tension derivatives concerning the current displacements in the FE model are correctly integrated into the total stiffness matrix to avoid instability. A sensitivity analysis of the number of samples used in the MC model and the combination of time step sizes confirmed the accuracy and robustness of MusAsi, and we concluded that the combination of a 1.25 ms FE time step and 0.005 ms MC multiple steps using a few hundred motor proteins in each finite element was appropriate in the tradeoff between accuracy and computational time. Furthermore, for a biventricular FE model consisting of 45,000 tetrahedral elements, one heartbeat could be computed within 1.5 h using 320 cores of a conventional parallel computer system. These results support the practicality of MusAsi for uses in both the basic research of the relationship between molecular mechanisms and cardiac outputs, and clinical applications of perioperative prediction.

Introduction

Demands for the prediction of outcomes from various types of operations are emerging in clinical problems of heart disease. At the present time, we can reconstruct a precise three-dimensional (3D) model from computed tomography or magnetic resonance imaging data for an individual patient, and use it for perioperative simulations to help doctors to choose the best among various possible operations. However, there are still some difficulties in the modeling of excitation-contraction coupling, even if we can precisely predict the excitation propagation from patient electrocardiogram (ECG) data (Okada et al., 2013). Such difficulties are because the macroscopic shortening in the fiber orientation is fed back to the extremely large stochastic combination consisting of various states of contractile proteins (Figure 1), and their stochastic responses depend on individual scenarios that include their neighbors (Figure 2). Although previous efforts to construct numerical models using a type of mean field approximation have successfully reproduced specific tissue-level phenomena (Hunter et al., 1998; Niederer et al., 2006; Negroni and Lascano, 2008; Rice et al., 2008; Guérin et al., 2011; Chapelle et al., 2012; Washio et al., 2012; Syomin and Tsaturyan, 2017; Regazzoni et al., 2018; Caruel et al., 2019), these models have not yet been fully exploited in real-life heart simulations. The uses of ordinary differential equation (ODE) models that adopt the phenomenological approximations of the force-pCa relationship and the force-velocity relationship have become mainstream instead (Smith et al., 2004; Kerckhoffs et al., 2007; Gurev et al., 2011; Shavik et al., 2017; Dabiri et al., 2019; Azzolin et al., 2020; Regazzoni et al., 2020). However, these approaches appear to have difficulties, particularly in reproducing the realistic relaxation phase that is important to ease the influx of blood from the atria to the ventricles.

FIGURE 1

Figure 1. (A) Tetrahedral element in (B) the ventricle wall of the FE model. The cylinder in panel (A) indicates the fiber orientation f that is shown in panel (B) using the integral curves colored according to the longitudinal component of f (blue: –1, red: +1). (C) The half-sarcomere model was assigned for each element. The active tension T_act is given by summing the forces of the binding myosin molecules composed of the myosin head (ellipse), and the lever arm (bar) and rod (spring). In the half-sarcomere model, the C-line and thick filaments (green) are fixed. The Z-line and thin filament (gray) slide with the half-sarcomere (HS) shortening velocity $-S{L}_0\stackrel{.}{\mathrm{\lambda }}/2$, where the stretch rate $\stackrel{.}{\mathrm{\lambda }}$ along the fiber orientation f is obtained from the FE model. The active tension T_act is given by summing the forces generated by the binding myosin molecules, and it is used to define the macroscopic active stress tensor S_act that drives the heartbeat in the FE model.

FIGURE 2

Figure 2. State transition Monte Carlo model of (A) the myosin molecule and (B) the T/T unit. The myosin molecules in either the N_XB or P_XB states are assumed to be detached. The rate constants in the T/T unit state transition are affected by the states of the myosin heads (MHs) below it. The rate constant factors k_np and k_pn between N_XB and P_XB are affected by the state of the T/T unit above it. ng is an integer that takes the values 0, 1, or 2, according to the number of neighboring MHs attached. γ = 40 was adopted to model the co-operativity of the MHs. The forward transitions from XB_PreR to XB_PostRs via XB_PostR1 are called ‘power strokes,’ whereas the back transitions are called ‘reverse strokes.’ The MHs connected to the extremely strained myosin rods are detached to N_ATP. (C) The typical transient of [Ca²⁺] applied to the T/T unit state transition in panel (B). (D) The strain energy W of the myosin rod and the non-linear force dW/dx.

Two major problems exist in directly coupling a stochastic molecular model and a living heart model: (i) the time scale difference between the two models; and (ii) the treatment of active stiffness in the heart model associated with the stochastic cross-bridge activity in the molecular model. Regarding the first problem, cross-bridge activity is typically modeled either using an ODE model or a Monte Carlo (MC) model that requires a small time step of microsecond order, whereas millisecond order is appropriate for the heart model discretized by the finite element method (FEM) in terms of the computational load and communication overhead. Regarding the second problem, an implicit method is typically applied for the FE model because of the strong anisotropic and volumetric stiffness of living tissue. When a muscle is excited, active stiffness associated with cross-bridge activity is generated in the fiber orientation (Figure 1B). This active stiffness is much greater than passive stiffness, with the exception of the volumetric stiffness of the incompressibility. Thus, if the prescribed active tensions computed in the cross-bridge model are explicitly applied to the FE model, active stiffness is not taken into account in the total stiffness matrix of the Newton iteration, which causes some problems of either convergence or accuracy. Such a problem was analyzed by Regazzoni and Quarteroni (2021), and the instability was fixed by introducing an appropriate active stiffness. However, their study was limited to some phenomenological ODE models that cannot reproduce a spontaneous oscillation (SPOC) (Ishiwata et al., 2010; Kagemoto et al., 2015) in low Ca²⁺ concentrations, as Regazzoni et al. (2021) remarked. By contrast, we successfully reproduced a SPOC (Washio et al., 2019; Shintani et al., 2020) using our MC cross-bridge model (Figure 2), and we addressed the similarity between the rapid lengthening of sarcomeres in the SPOC and the quick relaxation of cardiac muscle at early diastole in the cardiac cycle. Therefore, in this study, we focus on the stability in the direct coupling of the MC and FE models. The similarities of the cross-bridge dynamics in the relaxation phases between the SPOC and the biventricular FE simulations demonstrate the usefulness of our scheme both in areas of basic research and clinical applications.

Materials and Methods

Coupling of the MC Model and the FEM Model

In this study, we apply a multiple time step approach in which about 100 time steps of the MC model are performed within a single time step of the FE model to reduce the computational load and communication overhead (Figure 3). In our approach, to update the variables in both the MC and FE models from the FE time step at T to the next time step at T + ΔT, first, the stretch λ_T and stretch rate ${\stackrel{.}{\mathrm{\lambda }}}_T$ in the fiber orientation of the FE model are used as the initial half-sarcomere length (HSL) λ_T⋅SL₀/2 at T and its shorting velocity $-{\stackrel{.}{\mathrm{\lambda }}}_T\cdot S{L}_0/2$ in the time interval [T,TΔT] of the half-sarcomere model (Figure 3A), where SL₀ is the sarcomere length under the unloaded condition. This filament sliding information is used to calculate the myosin rod strains that are referenced to compute the state transitions of binding myosin molecules. Then, the active tension T_{act,[T,T+ΔT]} and associated stiffness ∂⁡T_{act,[T,T + ΔT]}/∂⁡λ_{T + ΔT} are computed iteratively in Newton iterations by summing the individual contributions of binding myosin molecules. The myosin rod strains are recomputed using the interpolated stretches between λ_T and λ_{T + ΔT} in the time interval [T,T + ΔT], whereas the myosin molecule states already computed in the MC steps are fixed (Figure 3B) for convergence in the Newton iterations. Thus, the stretch λ_{T + ΔT} is implicitly integrated into the active tension, which results in the appropriate evaluation of active stiffness in the FE model and the stability of the Newton iterations. Hereafter, we call our approach the “Multiple step Active stiffness integration” (MusAsi) scheme.

FIGURE 3

Figure 3. Difference between the sliding distances of the thin filament in (A) the MC state transition steps and (B) the FE active tension integration steps. An FE time interval [T,T + ΔT] is divided into n MC multiple steps (k = 1,⋯,n). The vertical blue line is the position of the Z-line under the unloaded condition. The sliding distances in [T,T + ΔT] are extrapolated from the sliding velocity $-{\stackrel{.}{\mathrm{\lambda }}}_T\cdot S{L}_0/2$ at T in the MC steps, whereas they are interpolated from the stretches λ_T and λ_T+ΔT at T and T + ΔT in the active tension integration steps, respectively. The gray arrows around the time axis indicate the flow of the computational process. The MC state transition steps are performed only once from T to T + ΔT, whereas the FE active tension integration steps are iterated until the convergence of Newton iterations, where λ_T+ΔT is computed from the updated displacement u_{T + ΔT}. Note that u_{T + ΔT} is initialized with u_T at the beginning of the Newton iterations in this study.

In our previous works (Washio et al., 2013, 2016), the active tension T_{act,[T,T+ΔT]} was implicitly computed by assuming the shorting velocity $-{\stackrel{.}{\mathrm{\lambda }}}_{T+\mathrm{\Delta }T}\cdot S{L}_0/2$ of the half-sarcomere model in the time interval [T,T + ΔT], whereas in the MusAsi scheme, it is computed by assuming −(λ_{T + ΔT}−λ_T)⋅SL₀/2ΔT. Although both approaches produce almost the same result, the previous approach is based on the velocities of the continuum, which result in the inconsistent stiffness of the active tension. By contrast, active stiffness in the MusAsi scheme is proportional to the total stiffness of the binding myosin molecules in the half-sarcomere model, thus the interpretation is consistent with our intuition. In the following, we introduce the MC model applied in this study and the details of the MusAsi scheme.

MC Cross-Bridge Model

We briefly present an overview of our MC cross-bridge model (Washio et al., 2016) used in this study. We provide the details in Supplementary Material 1. A myosin molecule in our cross-bridge model has three non-binding states (N_XB, P_XB, and N_ATP) and three strong binding states (XB_PreR, XB_PostR1, and XB_PostR2) (Figure 2A). Ca²⁺-sensitivity is reproduced using the state transitions in the troponin/tropomyosin (T/T) units on the thin filament (Figure 2B). The coefficients k_np and k_pn in the rate constants between the non-binding state N_XB and the weakly binding state P_XB are changed according to the state of the T/T unit above the myosin molecule. Co-operativity in the nearest neighbor interactions is incorporated with the factors γ^ng and γ−^ng to reproduce the force-pCa²⁺ relationship (Rice et al., 2003), where γ= 40 is used, and ng = 0, 1, or 2 is the number of neighboring myosin molecules either in the weakly binding (P_XB) or strong binding (XB_PreR, XB_PostR1, and XB_PostR2) states. We assume that one real thin filament in the 3D arrangement corresponds to two thin filaments in our half-sarcomere model. This is because we assume that co-operative behavior exists along the tropomyosin and tropomyosin molecules wrapped around the thin filament in a double spiral manner, and only one of the spirals is considered in our half-sarcomere model. As shown in the “Results” section, this co-operative mechanism contributes to almost completely removing the population of binding states in the diastolic phase in which nearly 10% of the peak Ca²⁺ concentration remains (Figure 2C).

Contraction force is generated by the power stroke transitions in which the strain of the myosin rod increases by s₁ and s₂ in the first and second strokes, respectively. In our model, the rate constants of the power and reverse strokes are given by functions of the rod strain x (the displacement from the unloaded state) so that the Boltzmann equilibrium condition is fulfilled:

$$\frac{h_{f,i}\left(x\right)}{h_{b,i}\left(x+s_i\right)}=\text{exp}\left(\frac{-W\left(x\right)-E_{i-1}+W\left(x+s_i\right)+E_i}{k_{B}T}\right),(1)$$

where E_i–1 and E_i are the free energies before and after the power stroke under the unloaded condition, respectively. W is the strain energy of the myosin rod (Figure 2D). With the power stroke, the free energy decrease E_i−1−E_i is transferred to the increase of strain energy W(x + s_i)−W(x). The strain energy is used for the external work via the half-sarcomere shortening that corresponds to muscle shortening in the fiber direction. In this study, the individual rate constants for the power and reverse strokes are determined by

$$\{\begin{array}{c}\hfill h_{f,i}\left(x\right)=h_i\text{exp}\left(-{\displaystyle \frac{E_{i-1}-W\left(x+s_i\right)-E_i}{k_{B}T}}\right)\hspace{1em}\hspace{1em}(2)\hfill \\ \hfill h_{b,i}(x+s_i)=h_i\text{exp}(-{\displaystyle \frac{W\left(x\right)}{k_{B}T}}).\hspace{1em}\hspace{1em}\hspace{0.5em}\hspace{0.25em}(3)\hfill \end{array}.$$

To achieve stable MC steps, if either h_f,i or h_b,i exceeds the maximum rate r_max = 100,000 [1/s], it is replaced by r_max and the other parameter is modified so that Eq. 1 is fulfilled.

MusAsi Scheme

For the MC model, the time step size Δt must be chosen so that it is sufficiently smaller than the reciprocals of the rate constants. In our case, the choice Δt~5μs is appropriate. By contrast, the FE time step size ΔT∼1ms is sufficient to catch the time transients of macroscopic variables, such as the ventricular cavity volume and pressure. Because the linear solution, which requires many communications among processes, must be performed in each Newton iteration, a FE time step size ΔT that is much larger than Δt∼5μs is desirable. This leads to the use of an approach in which multiple MC steps are performed in a single FE step.

The feedback to the state transitions in the MC model from the dynamics of the FE model is given by the stretch in the fiber orientation f. In the FE model, the stretch and stretch rate are given by

$$\mathrm{\lambda }=||\mathbf{Ff}||(4)$$

$$\stackrel{.}{\mathrm{\lambda }}=\frac{1}{\mathrm{\lambda }}\left(\text{Ḟ}\mathbf{\text{f}}\right)\cdot \left(\mathbf{Ff}\right)(5)$$

where F = I + ∂⁡u/∂⁡X is the deformation gradient tensor defined for the displacements u = u(T,X) from the unloaded configuration. The macroscopic information of the stretch is provided at the two ends of the time interval at T and T + ΔT. The information at T + ΔT is determined when the macroscopic displacement u = u(T + ΔT,X) is determined, which can be computed only if the active stress during [T,T + ΔT] is provided. Thus, the state transitions of myosin molecules in MC model must be computed before the FE step from T to T + ΔT. Hereafter, we denote the time using a subscript, if necessary, for example, λ_T and ${\stackrel{.}{\mathrm{\lambda }}}_T$. A simple approach to perform the MC steps with the time step size Δt = ΔT/n is to define the stretch $\tilde{\mathrm{\lambda }}$ at the k-th step as

$${\tilde{\mathrm{\lambda }}}_{T+k\mathrm{\Delta }t}\equiv {\mathrm{\lambda }}_T+k\mathrm{\Delta }t{\stackrel{.}{\mathrm{\lambda }}}_T.(6)$$

In this case, the rod strain ${\tilde{x}}_{ij}$ (displacement from the unloaded position) of the (i,j)-th myosin molecule is given by the following if it is in the binding state:

$${\tilde{x}}_{ij,T+k\mathrm{\Delta }t}=x_{A,ij,T+k\mathrm{\Delta }t}+s_{ij,T+k\mathrm{\Delta }t}+\frac{S{L}_0}{2}\left({\tilde{\mathrm{\lambda }}}_{T+k\mathrm{\Delta }t}-{\tilde{\mathrm{\lambda }}}_{A,ij,T+k\mathrm{\Delta }t}\right)(7)$$

where the integers i = 1,⋯,N_M and j = 1,⋯,N_F denote the index of a myosin molecule in a filament and the index of a thin filament in the half-sarcomere model imbedded in a single finite element in the FE model (Figure 1), respectively. The variables x_A,ij and ${\tilde{\mathrm{\lambda }}}_{A,ij}$ are the initial strain and the stretch at the most recent attachment (the transition from P_XB to XB_PreR), respectively. The initial strain x_A,ij is yielded probabilistically from the Boltzmann distribution exp(−W(x)/k_BT). The variable s_ij = 0,s₁ or s₁ + s₂ is the power stroke distance. The third term on the right-hand side is the sliding distance between the thin filament and the thick filament after the attachment. The constant SL₀ is the sarcomere length under the unloaded condition. In the MC computation, a state transition of a binding myosin molecule ij at k-th step is computed based on the rate constant determined by the rod strain ${\tilde{x}}_{ij,T+k\mathrm{\Delta }t}$.

Once the state transitions in the MC model in the time interval [T,T + ΔT] are computed, the active tension T_{act,[T,T+ΔT]} of the FE model is calculated by summing all the forces produced by individual myosin molecules so that the impulses of both scales are the same:

$$T_{act,[T,T+\mathrm{\Delta }T]}=\frac{2\cdot R_S}{S{A}_0\cdot n_F\cdot n}\sum _{j=1}^{N_F}\sum _{i=1}^{N_M}\sum _{k=1}^n{\mathrm{\delta }}_{A,ij,T+\mathrm{\kappa }\Delta t}\frac{dW}{dx}\left(x_{ij,T+k\mathrm{\Delta }t}\right),(8)$$

where SA₀ and R_S are the cross-sectional area per thin filament and the sarcomere volume ratio under the unloaded condition, respectively. The numerator is multiplied by a factor 2 because we consider a half the myosin molecules (N_M = 38), which are accessible along a single spiral in the thin filament for the purpose of co-operative attach-detach control along the tropomyosin. The variables {δ_{A,ij,T + κΔt} = 0: non-binding, =1: binding} are obtained from the MC model. Regarding the forces generated by the individual myosin molecules, if we define the rod strain as $x_{ij,T+k\mathrm{\Delta }t}\equiv {\tilde{x}}_{ij,T+k\mathrm{\Delta }t}$, the active tension T_{act,[T,T+ΔT]} is determined regardless of the stretch λ_{T + ΔT} at T + ΔT. Therefore, the active tension stiffness dT_{act,[T,T + ΔT]}/dλ_{T + ΔT} associated with the binding myosin molecules in the half-sarcomere model is not taken into account in the total stiffness matrix used in the FE Newton iterations. This approach, which considers the active tension explicitly, results in the instability of the numerical solution, as seen in the “Results” section. Therefore, in the MusAsi scheme, the rod strains {x_{ij,T + kΔt}} used to determine the active tension are given by

$$x_{ij,T+k\mathrm{\Delta }t}\equiv x_{A,ij,T+k\mathrm{\Delta }t}+s_{ij,T+k\mathrm{\Delta }t}+\frac{S{L}_0}{2}\mathrm{}\left({\mathrm{\lambda }}_{T+k\mathrm{\Delta }t}-{\mathrm{\lambda }}_{A,ij,T+k\mathrm{\Delta }t}\right),(9)$$

where the variables {x_{A,ij,T + kΔt}s_{ij,T + kΔt}} produced in the MC steps k = 1,⋯,n are used, whereas the current stretch λ_{T + kΔt} and the stretch λ_{A,ij,T + kΔt} at the most recent attachment are sequentially redefined from the stretch λ_{T + ΔT} at the time step T + ΔT as follows:

$${\mathrm{\lambda }}_{T+k\mathrm{\Delta }t}\equiv {\mathrm{\lambda }}_T+\frac{k}{n}\left({\mathrm{\lambda }}_{T+\mathrm{\Delta }T}-{\mathrm{\lambda }}_T\right)(10)$$

$${\mathrm{\lambda }}_{A,ij,T+k\mathrm{\Delta }t}\{\begin{array}{c}\hfill {\mathrm{\lambda }}_{T+k\mathrm{\Delta }T},{\mathrm{\delta }}_{A,ij,T+\mathrm{\kappa }\mathrm{\Delta }t}=1and{\mathrm{\delta }}_{A,ij,T+\left(k-1\right)\mathrm{\Delta }t}=0\hfill \\ \hfill {\mathrm{\lambda }}_{A,ij,T+\left(k-1\right)\mathrm{\Delta }t},otherwise.\hfill \end{array}(11)$$

By including λ_{T + ΔT} in the definition of rod strain {x_{ij,T + kΔt}}, the stiffness associated with cross-bridge activity at T + ΔT is properly represented as

$$\frac{\partial T_{act,[T,T+\text{Δ}T]}}{\partial {\text{λ}}_{T+\text{Δ}T}}=\frac{2}{S{A}_0\cdot N_F\cdot n}{\displaystyle \sum }_{j=1}^{N_F}{\displaystyle \sum }_{i=1}^{N_M}{\displaystyle \sum }_{k=1}^n{\text{ δ}}_{A,ij,T+\text{κΔ}t}\frac{d^{2}W}{d{x}^2}(x_{ij,T+k\text{Δ}t})\frac{\partial x_{ij,T+k\text{Δ}t}}{\partial {\text{λ}}_{T+\text{Δ}T}}(12)$$

where the derivative with respect to λ_{T + ΔT} is given by

$$\frac{\partial x_{ij,T+k\mathrm{\Delta }t}}{\partial {\mathrm{\lambda }}_{T+\mathrm{\Delta }T}}=\frac{S{L}_0}{2}\left(\frac{\partial {\mathrm{\lambda }}_{T+k\mathrm{\Delta }t}}{\partial {\mathrm{\lambda }}_{T+\mathrm{\Delta }T}}-\frac{\partial {\mathrm{\lambda }}_{A,ij,T+k\mathrm{\Delta }t}}{\partial {\mathrm{\lambda }}_{T+\mathrm{\Delta }T}}\right)(13)$$

with

$$\frac{\partial {\mathrm{\lambda }}_{T+k\mathrm{\Delta }t}}{\partial {\mathrm{\lambda }}_{T+\mathrm{\Delta }T}}=\frac{k}{n}(14)$$

and

$$\frac{\partial {\mathrm{\lambda }}_{A,ij,T+k\mathrm{\Delta }t}}{\partial {\mathrm{\lambda }}_{T+\mathrm{\Delta }T}}=\{\begin{array}{c}\hfill \frac{k}{n},{\mathrm{\delta }}_{A,ij,T+\mathrm{\kappa }\mathrm{\Delta }t}=1and{\mathrm{\delta }}_{A,ij,T+(k-1)\mathrm{\Delta }t}=0\hfill \\ \hfill \frac{\partial {\mathrm{\lambda }}_{A,ij,T+(k-1)\mathrm{\Delta }t}}{\partial {\mathrm{\lambda }}_{T+\mathrm{\Delta }T}},otherwise\hfill \end{array}(15)$$

starting from $\frac{\partial {\mathrm{\lambda }}_{A,ij,T+k\mathrm{\Delta }t}}{\partial {\mathrm{\lambda }}_{T+\mathrm{\Delta }T}}=0$ at k=0. From Eqs 14, 15, we obtain

$$\frac{\partial x_{ij,T+k\mathrm{\Delta }t}}{\partial {\mathrm{\lambda }}_{T+\mathrm{\Delta }T}}\ge 0.(16)$$

Therefore, the stiffness ∂⁡T_{act,[T,T + ΔT]}/∂⁡λ_{T + ΔT} is always non-negative, provided the potential W is convex downward (d²W/dx²≥0).

Active Stress and Stiffness in the FE Model

The infinitesimal external work per unit volume required to make an infinitesimal increment of stretch δλ against the active tension T_act is given by

$$\mathrm{\delta }{W}_{ext}=T_{act}\mathrm{\delta }\mathrm{\lambda }.(17)$$

From the relationship between the stretch λ and the Green–Lagrange strain tensor E = (F^TF−I)/2, we obtain

$$\mathbf{\text{f}}\otimes \mathbf{\text{f}}:\mathbf{\text{E}}=\frac{1}{2}\left({\mathrm{\lambda }}^2-1\right).(18)$$

Thus, the infinitesimal increment of stretch is represented by

$$\mathrm{\delta }\mathrm{\lambda }=\frac{1}{\mathrm{\lambda }}\mathbf{\text{f}}\otimes \mathbf{\text{f}}:\mathrm{\delta }\mathbf{\text{E}}.(19)$$

Therefore, if we define the second Kirchhoff active stress tensor as

$${\mathbf{\text{S}}}_{act}=\frac{T_{act}}{\mathrm{\lambda }}\mathbf{\text{f}}\otimes \mathbf{\text{f}},(20)$$

the infinitesimal work is represented by δW_act = S_act:δE.

The derivative of the infinitesimal work is given by

$${\mathrm{\delta }}^2{W}_{act}=\frac{\partial T_{act}}{\partial \mathrm{\lambda }}\mathrm{\delta }{\mathrm{\lambda }}^2+T_{act}{\mathrm{\delta }}^2\mathrm{\lambda }.(21)$$

If we assume T_act≡T_{act,[T,T + ΔT]} and λ≡λ_{T + ΔT}, the first term on the right-hand side is non-negative from Eq. 16. The second term is also non-negative, provided the active tension T_act is non-negative because the Hessian of λ is represented as

$${\mathrm{\delta }}^2\mathrm{\lambda }=\frac{1}{\mathrm{\lambda }}(\mathbf{\text{f}}\otimes \mathbf{\text{f}}:{\mathrm{\delta }}^2\mathbf{\text{E}}-\frac{1}{{\mathrm{\lambda }}^2}{(\mathbf{\text{f}}\otimes \mathbf{\text{f}}:\mathrm{\delta }\mathbf{\text{E}})}^2)=\frac{1}{\mathrm{\lambda }}\left(\left(\mathrm{\delta }\mathbf{\text{Ff}}\right)\cdot \left(\mathrm{\delta }\mathbf{\text{Ff}}\right)-{\left(\mathbf{\text{a}}\cdot \left(\mathrm{\delta }\mathbf{\text{Ff}}\right)\right)}^2\right),(22)$$

with the normal fiber orientation vector in the current coordinate a = Ff/λ. Therefore, δ²W_act is non-negative, provided the active tension T_act is non-negative. This guarantees the stability of the MusAsi scheme during the contraction phase. Some cases of negative active tension T_act < 0may exist. In this case, the negative stiffness appears in the subspace of variations of deformation gradients δF that satisfy a⋅(δFf) = 0. However, because the negative tension is typically small compared with the positive term from the mass, viscosity, and passive stiffness, stability is guaranteed for an appropriately small time step size ΔT∼1ms, in our experience. When some impaired cross-bridge models were tested, we also observed that the deletion of the second term on the right-hand side of Eq. 21 from the stiffness matrix for negative active tensions (T_act < 0) further stabilized the convergence of Newton iterations in the MusAsi scheme. In our previous approach (Washio et al., 2016), because ${\mathrm{\lambda }}_{T+k\mathrm{\Delta }t}\equiv {\mathrm{\lambda }}_T+k\mathrm{\Delta }t{\stackrel{.}{\mathrm{\lambda }}}_{T+\mathrm{\Delta }T}$ was used instead of Eq. 10, the first term in Eq. 21 was replaced with $\left(\partial T_{act}/\partial \stackrel{.}{\mathrm{\lambda }}\right)\mathrm{\delta }\mathrm{\lambda }\partial \stackrel{.}{\mathrm{\lambda }}$. Therefore, the interpretation of this term as the stiffness was difficult, although we did not have any convergence difficulty.

Newton Iterations for the FE Model

In this study, the FE biventricular model was connected with the systemic and pulmonary circulation models, and the transfer of blood volume using these circulation models was described only by the volumetric changes of ventricular cavities. Thus, the combined system of equations for the FE model is given by the following six formulas for the biventricular FE (Eqs 23–26), systemic circulation (Eq. 27), and pulmonary circulation models (Eq. 28):

$$\{\begin{array}{l}{\displaystyle {\int }_{\Omega }^{}\text{δ}u\cdot \rho űd\Omega +}{\displaystyle {\int }_{\Omega }^{}\text{δ}E:(S-p{C}^{-1})d\Omega }\\ \begin{array}{ll}-P_L{\displaystyle {\int }_{{\Gamma }_{\text{L}}}^{}\text{δ}u\cdot nd{\Gamma }_L-P_R{\displaystyle {\int }_{{\Gamma }_R}^{}u\cdot nd{\Gamma }_R=0}}\hfill & (23)\hfill \\ {\displaystyle {\int }_{\Omega }^{}\text{δ}p}\left((J-1)+\frac{p}{k}\right)d\Omega =0\hfill & (24)\hfill \\ {\displaystyle {\int }_{{\Gamma }_L}^{}\dot{u}\cdot nd{\Gamma }_L-(F_{MI}-F_{AO})=0}\hfill & (25)\hfill \\ {\displaystyle {\int }_{{\Gamma }_R}^{}\dot{u}\cdot nd{\Gamma }_R-(F_{TR}-F_{PA})=0}\hfill & (26)\hfill \\ C_S(P_L,F_{AO},Q_S,F_{TR},P_R)=0\hfill & (27)\hfill \\ C_P(P_R,F_{PA},Q_P,F_{MI},P_L)=0,\hfill & (28)\hfill \end{array}\end{array}$$

where J = detF is the Jacobian, p is the hydrostatic pressure in the ventricular walls, C = F^TF is the right Cauchy-Green deformation tensor, κ is the bulk modulus, and P_L = LVP and P_R = RVP are blood pressure in the left and right ventricular cavities, respectively. is the muscle domain in the reference configuration, whereas Γ_L and Γ_R are the blood–muscle interfaces of the left and right ventricles, respectively, in the configuration at time T, and n is the normal unit vector directed from the cavity to the muscle at these surfaces. The Dirichlet boundary condition u_T(X)≡0 is imposed on the boundary nodes around the valve rings. The second Piola–Kirchhoff stress tensor S consists of the active, passive, and viscous stresses:

$$\mathbf{\text{S}}\equiv {\mathbf{\text{S}}}_{act}+{\mathbf{\text{S}}}_{pas}+{\mathbf{\text{S}}}_{vis}(29)$$

where S_act is given by Eq. 20, and S_pas and S_vis are the passive and viscous stresses, respectively. The ventricle blood pressures P_L and P_R are determined through their interactions with the circulatory system of the body. F_MI, F_AO, F_TR, and F_PA are the flow rates through the mitral, aortic, tricuspid, and pulmonary valves, respectively. Q_S and Q_P are the variables in the systemic and pulmonary circulatory systems, respectively. We provide the details and the parameters applied in this study in Supplementary Material 2.2.

The time integration of the combined system composed of Eqs 23–28 were performed with the Newmark-beta scheme:

$$\{\begin{array}{ll}{\dot{U}}_{T+\Delta T}={\dot{U}}_T+\Delta T({\gamma }^{\ddot{U}}T+\Delta T+(1-\gamma ){\ddot{U}}_T)\hfill & (30)\hfill \\ U_{T+\Delta T}=U_T+\Delta T{\dot{U}}_T+\Delta T^2(\text{β}{\ddot{U}}_{T+\Delta T}+(1/2-\text{β}){\ddot{U}}_T)\hfill & (31)\hfill \\ R(U_{T+\Delta T},{\dot{U}}_{T+\Delta T},{\ddot{U}}_{T+\Delta T})=0.\hfill & (32)\hfill \end{array}$$

where the vector U contains all variables as follows:

$$\mathbf{\text{U}}=\left(\begin{array}{c}\hfill \begin{array}{c}\hfill \mathbf{\text{u}}\hfill \\ \hfill p\hfill \\ \hfill P_L\hfill \end{array}\hfill \\ \hfill P_R\hfill \\ \hfill *\hfill \\ \hfill *\hfill \\ \hfill *\hfill \\ \hfill *\hfill \\ \hfill *\hfill \\ \hfill *\hfill \end{array}\right),\stackrel{.}{\mathbf{\text{U}}}=\left(\begin{array}{c}\hfill \begin{array}{c}\hfill \stackrel{.}{\mathbf{\text{u}}}\hfill \\ \hfill *\hfill \\ \hfill *\hfill \end{array}\hfill \\ \hfill *\hfill \\ \hfill F_{MI}\hfill \\ \hfill F_{AO}\hfill \\ \hfill F_{PA}\hfill \\ \hfill F_{TR}\hfill \\ \hfill {\mathbf{\text{Q}}}_P\hfill \\ \hfill {\mathbf{\text{Q}}}_S\hfill \end{array}\right),\ddot{\mathbf{\text{U}}}=\left(\begin{array}{c}\hfill \begin{array}{c}\hfill \ddot{\mathbf{\text{u}}}\hfill \\ \hfill *\hfill \\ \hfill *\hfill \end{array}\hfill \\ \hfill *\hfill \\ \hfill *\hfill \\ \hfill *\hfill \\ \hfill *\hfill \\ \hfill *\hfill \\ \hfill {\stackrel{.}{\mathbf{\text{Q}}}}_P\hfill \\ \hfill {\stackrel{.}{\mathbf{\text{Q}}}}_S\hfill \end{array}\right),(33)$$

and the function R involves all the equilibrium and constraint conditions in Eqs 23–28. The missing components denoted by “^∗” in Eq. 33 do not appear in the function R. Although these components are calculated following the rules of time interpolation in Eqs 30, 31, they do not have any physical meaning.

Eqs 30–32 are simultaneously solved implicitly using Newton iterations as follows:

Set the initial guess: ${\mathbf{\text{U}}}_{T+\mathrm{\Delta }T}^{\left(0\right)}={\mathbf{\text{U}}}_T$,${\stackrel{.}{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(0\right)}={\stackrel{.}{\mathbf{\text{U}}}}_T$, ${\ddot{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(0\right)}={\ddot{\mathbf{\text{U}}}}_T$

Iterate k = 0,1,2,⋯ until ||R^(k)||≤ε

Compute the residual and the stiffness matrix:

$${\mathbf{\text{R}}}^{\left(k\right)}=\mathbf{\text{R}}({\mathbf{\text{U}}}_{T+\mathrm{\Delta }T}^{\left(k\right)},{\stackrel{.}{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(k\right)},{\ddot{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(k\right)})$$

$${\mathbf{\text{K}}}^{\left(k\right)}=\frac{\partial \mathbf{\text{R}}}{\partial \mathbf{\text{U}}},{\mathbf{\text{C}}}^{\left(k\right)}=\frac{\partial \mathbf{\text{R}}}{\partial \stackrel{.}{\mathbf{\text{U}}}},{\mathbf{\text{M}}}^{\left(k\right)}=\frac{\partial \mathbf{\text{R}}}{\partial \ddot{\mathbf{\text{U}}}}$$

Solve the linear system:

$$\left({\mathbf{\text{M}}}^{\left(k\right)}+\mathrm{\gamma }\mathrm{\Delta }T{\mathbf{\text{C}}}^{\left(k\right)}+\mathrm{\beta }\mathrm{\Delta }{T}^2{\mathbf{\text{K}}}^{\left(k\right)}\right)\mathrm{\Delta }{\ddot{\mathbf{\text{U}}}}^{\left(k\right)}=$$

$$\{\begin{array}{c}\hfill -{\mathbf{\text{R}}}^{\left(0\right)}-{\mathbf{\text{C}}}^{\left(0\right)}\mathrm{\Delta }T{\ddot{\mathbf{\text{U}}}}_T-{\mathbf{\text{K}}}^{\left(0\right)}\left(\mathrm{\Delta }T{\stackrel{.}{\mathbf{\text{U}}}}_T+\frac{1}{2}\mathrm{\Delta }{T}^2{\ddot{\mathbf{\text{U}}}}_T\right),k=0\hfill \\ \hfill -{\mathbf{\text{R}}}^{\left(k\right)},k\ge 1\hfill \end{array}$$

Update:

$${\ddot{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(k+1\right)}={\ddot{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(k\right)}+\mathrm{\Delta }{\ddot{\mathbf{\text{U}}}}^{\left(k\right)}$$

$${\stackrel{.}{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(k+1\right)}=\{\begin{array}{c}\hfill {\stackrel{.}{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(0\right)}+\mathrm{\Delta }T\left(\mathrm{\gamma }\mathrm{\Delta }{\ddot{\mathbf{\text{U}}}}^{\left(k\right)}+{\ddot{\mathbf{\text{U}}}}_T\right),k=0\hfill \\ \hfill {\stackrel{.}{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(k\right)}+\mathrm{\gamma }\mathrm{\Delta }T\mathrm{\Delta }{\ddot{\mathbf{\text{U}}}}^{\left(k\right)},k\ge 1\hfill \end{array}$$

$${\mathbf{\text{U}}}_{T+\mathrm{\Delta }T}^{\left(k+1\right)}\{\begin{array}{c}\hfill {\mathbf{\text{U}}}_{T+\mathrm{\Delta }T}^{\left(0\right)}\mathrm{\Delta }T{\stackrel{.}{\mathbf{\text{U}}}}_T+\mathrm{\Delta }{T}^2\left(\mathrm{\beta }\mathrm{\Delta }{\ddot{\mathbf{\text{U}}}}^{\left(0\right)}+1/2{\ddot{\mathbf{\text{U}}}}_T\right),k=0\hfill \\ \hfill {\mathbf{\text{U}}}_{T+\mathrm{\Delta }T}^{\left(k\right)}+\mathrm{\Delta }{T}^2\mathrm{\beta }\mathrm{\Delta }{\ddot{\mathbf{\text{U}}}}^{\left(k\right)},k\ge 1\hfill \end{array}$$

Because the initial guesses ${\mathbf{\text{U}}}_{T+\mathrm{\Delta }T}^{\left(0\right)}$, ${\stackrel{.}{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(0\right)}$, and ${\ddot{\mathbf{\text{U}}}}_{T+\mathrm{\Delta }T}^{\left(0\right)}$ do not fulfill the interpolation rules in Eqs 30, 31, whereas the solutions after that do, different right-hand sides and update rules are applied. The key issue in the Newton iteration is the treatment of active stress S_act in the computation of stiffness matrix K = ∂⁡R/∂⁡U. The active force vector and the stiffness matrix associated with the active stress tensor S_act on an element e are given by

$${\mathbf{\text{F}}}_{act,e}={\int }_e{T}_{act}\frac{\partial \mathrm{\lambda }}{\partial {\mathbf{\text{u}}}_e}\mathit{d}\Omega (34)$$

$${\mathbf{\text{K}}}_{act,e}={\int }_e\frac{\partial T_{act}}{\partial \mathrm{\lambda }}{\left(\frac{\partial \mathrm{\lambda }}{\partial {\mathbf{\text{u}}}_e}\right)}^T\left(\frac{\partial \mathrm{\lambda }}{\partial {\mathbf{\text{u}}}_e}\right)\mathit{d}\Omega +{\int }_e{T}_{act}\frac{{\partial }^2\mathrm{\lambda }}{\partial \mathbf{\text{u}}_e{}^2}\mathit{d}\Omega (35)$$

where the derivatives of stretch λ with respect to nodal displacements u_e of the element e are given by Eqs 19, 22. In the MusAsi scheme, the active tension T_act and its derivative ∂⁡T_act/∂⁡λ must be computed from the current value of λ_{T + ΔT} as defined by Eqs 8, 12, respectively, in every Newton iteration step. The processing flows and the data transfers between the macro and micro processes are shown in Figure 4.

FIGURE 4

Figure 4. (A) Micro and (B) macro processes, and (C) communications between them to perform the MusAsi scheme. Before the Newton iterations, the displacement u_T,e and its time derivative ${\stackrel{.}{\mathbf{\text{u}}}}_{T,e}$ on individual elements {e} are sent from the macro process to the micro process to perform the MC steps in the micro process. In the Newton iterations, the current displacement u_{T + ΔT,e} is sent to the micro process to compute the active force vector F_act,e and the active stiffness matrix K_act,e on each element e.

Computer Resource and the FE Model

The tested biventricular FE model consisted of 45,000 tetrahedral elements, where the MINI(5/4c) element (Brezzi and Fortin, 1991) was adopted to avoid instability caused by the nearly incompressible condition in Eq. 24. Although the higher-order interpolation of MINI elements was applied to the displacement u to evaluate the integration associated with the passive stress tensor, standard linear interpolation ignoring the central node was adopted for active stress. Thus, it was sufficient to assign one half-sarcomere model to each element. The fiber-sheet architecture was constructed by applying the optimization algorithm in our previous work (Washio et al., 2020). The computations were performed using 20 nodes (320 cores) of a parallel computer system (Intel Xeon E-2670 [2.6 GHz], 16 cores/node; Intel, Santa Clara, CA, United States). In the typical MusAsi scheme in which 16 filaments (N_F = 16) were assigned to each half-sarcomere model, the computational time was about 1.26 h per heartbeat with ΔT = 1.25ms, and Δt = 5μs. The MS steps and active tension integration steps (Figure 4A) in the micro process took 0.58 and 0.42 h, respectively. The remaining computational time was almost occupied with the linear solutions (Washio and Hisada, 2011; Kariya et al., 2020) in the macro process (Figure 4B). Because 38 myosin molecules were arranged in each filament, 27 million myosin molecules were used in total.

The heart rate was set to 60 beats per minute, and the Ca²⁺-transient generated by the midmyocardial cell model proposed by ten Tusscher and Panfilov (2006) was applied (Figure 2C). Transmural delays were used that were determined by the distances from the endocardial surfaces of the left and right ventricles under a transmural condition velocity of 52 cm/s, as measured by Taggart et al. (2000).

From the numerical results, the output work from the aortic valve was evaluated as

$$W_{out}={\int }_0^{T_C}{F}_{AO}\cdot P_L\mathit{d}T,(36)$$

where T_C = 1 s is the cardiac cycle period. ATP consumption was also calculated by counting the transients to N_XB from XB_PostR2 or N_ATP (Figure 2A).

Results

Accuracy in Overall Cardiac Outcomes

The influence of the number of filaments N_F imbedded in one element on cardiac outcomes, such as pressure and volume, are shown in Table 1A and Figure 5. In the simulations, ΔT = 1.25 ms and Δt = 5 μs were applied so that 250 MC steps were performed within a single FE step. Although there was a little difference in the overall time transients of pressure and volume, even for N_F = 4, the difference between the minimum and maximum of these variables from N_F = 64 was less than 1%. The difference in ATP consumption from N_F = 64 was slightly larger than that of pressure and volume. However, it was also less than 1% for N_F = 16. Thus, it seemed to be sufficient to use 16 filaments to obtain the overall cardiac outcomes.

TABLE 1

Table 1. Influence of the filament number N_F (A) and FE time step size ΔT (B) on overall pumping performance and computational time.

FIGURE 5

Figure 5. Time transients of (A) left ventricular pressure (LVP) and (B) volume (LVV) for N_F = 4 (thin red lines), 16 (medium thick blue lines), and 64 (thick gray broken lines). ΔT = 1.25 ms and Δt = 5 μs were applied. These almost equal time transients were taken from the third cycle after the initial process of blood filling.

The influence of the FE time step size ΔT on the overall cardiac outcomes was also examined, as shown in Table 1B. The baseline of the MC time step size was given by Δt₀ = 5 μs, and the number of MC steps n performed within a single FE step was determined by

$$n\equiv \lfloor \frac{\mathrm{\Delta }T-0.5\mathrm{\Delta }{t}_0}{\mathrm{\Delta }{t}_0}\rfloor +1,(37)$$

where “⌊⌋” represents the floor function that rounds down after the decimal point. Therefore 250, 125, and 63 MC steps were performed with ΔT = ΔT₀, ΔT₀/2, and ΔT₀/4 (ΔT₀ = 1.25 ms), respectively. As with the number of filaments N_F, the difference with ΔT was sufficiently small. As the computational time increased, ΔT decreased because the total Newton iterations and the communication between the MC and FE models increased. This result suggests that the choice ΔT = 1.25 ms was sufficiently good and preferable in terms of the computational cost.

To confirm the necessity of the implicit approach for active tension, an explicit approach with various sizes of the FE time step ΔT was tested (Figure 6). In the explicit approach, only the calculation of active tension was modified so that it was determined using the strains $\left\{{\tilde{x}}_{ij,T+k\mathrm{\Delta }t}\right\}$ in Eq. 7 calculated in the MC steps from the stretch λ and stretch rate $\stackrel{.}{\mathrm{\lambda }}$ at time T instead of {x_{ij,T + kΔt}} in Eq. 9 calculated from the stretch λ at time T + ΔT. In this explicit case, only the second term on the right-hand side in Eqs 21, 35 was used to construct the stiffness matrix because neither the stretch λ or stretch rate $\stackrel{.}{\mathrm{\lambda }}$ at time T + ΔT was involved in the definition of T_act at time T + ΔT. Although there was no breakdown of the Newton iterations with the explicit approach, incorrect results appeared at certain times during the contraction phase. The smaller the time step size ΔT, the later the time at which the wrong result appeared. Additionally, finally, almost the same result as the implicit approach with ΔT = ΔT₀ = 1.25 ms was reproduced with quite a small time step ΔT = ΔT₀/128∼10 μs. This result supports both the numerical accuracy and computational efficiency of the MusAsi scheme because the stable explicit approach with ΔT = ΔT₀/128 required communication between the micro and macro processes and the linear solution in the macro process every MC step (n = 2) and, thus, a single beat computation took 70 h, whereas the MusAsi took only 1.26 h with n = 250 without loss of accuracy.

FIGURE 6

Figure 6. Time transients of (A) left ventricular pressure (LVP) and (B) volume (LVV) obtained by the MusAsi with the FE time step size ΔT = ΔT₀ (thick gray broken lines) and the explicit active tension approach with the finer time steps ΔT = ΔT₀/8 (blue lines), ΔT₀/32 (green lines), and ΔT₀/128 (red lines).

Accuracy of Local Dynamics

In clinical applications, not only the overall cardiac outputs, but also the local mechanical load and energy consumption are important for predicting a remote prognosis. To confirm the accuracy of local dynamics, the influence of the filament number N_F on the distribution of active tension T_act and ATP consumption at the peak of the systolic phase were examined (Figure 7). Although the discontinuities of active tensions at element boundaries was slightly noticeable, even for N_F = 64 (Figure 7A), it became inconspicuous when the elementwise variables were averaged at the nodes for N_F≥16 (Figures 7B,C). The distributions of the active tension values and the ATP consumption values over the entire cycle are shown in Supplementary Video 1. A more detailed comparison of the time transients of active tensions at a single element further indicated that the choice N_F = 16 was sufficient for analyzing the local mechanical load (Figure 8).

FIGURE 7

Figure 7. Distributions of active tension and ATP consumption at T = 0.2 s in the contraction phase in the middle cross-section perpendicular to the long axis. (A) The active tensions calculated in the individual elements are shown. The lines indicate the segmentation to elements. The black lines are element boundaries. (B) The active tensions in elements were averaged on the nodes in the FE mesh. (C) The ATP consumption in the interval [0.0 s, 0.2 s] for the individual elements was averaged on the nodes in the FE mesh.

FIGURE 8

Figure 8. Time transients of local active tensions at the Endo, Mid, and Epi sites in the lateral left ventricular wall. (A) The active tensions T_act at a single element are plotted for N_F = 4 (blue lines), 16 (green lines), and 64 (red lines). (B) The averaged values with surrounding elements within 2 mm are plotted.

It is somewhat counter-intuitive that the highest ATP consumption is not happening in the regions of highest active tension production (Figures 7B,C). Since the ATP consumption in Figure 7C is the cumulative value over the time interval [0.0 s, 0.2 s], it is difficult to find temporal relationship with the active tension. In Supplementary Material 3, the active tension, the stretch rate, and the ATP consumption rate at T = 0.2 s are shown. Here, due to the force-velocity relationship (Figure 10B), the higher the shortening velocity (negative stretch rate) is, the smaller the active tension gets. Because the shortening of half-sarcomere shifts the rod strain distribution to the negative direction (Figure 1C) resulting in the facilitation of power stroke (h_f,1 and h_f,2 in Figure 13A), the higher the shortening velocity is, the higher the ATP consumption rate gets. Therefore, the lowest ATP consumption rate is happening in the region of highest active tension production when the stretch is relatively uniform over the region.

Sensitivity of the Nearest Neighbor Co-operative Parameter

To confirm the importance of the neighboring co-operative mechanism for the relaxation phase in the cardiac cycle, pumping performances were compared for different co-operative parameters γ=40, 20, and 10 (Figures 9A,B). Because the Ca²⁺-concentration did not disappear, even in the diastolic phase (Figure 2C), active tensions were not completely removed with the impaired co-operativity (Figure 9C). Thus, the insufficient drop of left ventricular pressure (LVP) blocked the filling of blood through the mitral valve. The time transients of [XB_PostR2] indicate that even a small binding population less than 1% was likely to hamper the extension of the ventricular cavity (Figure 9D).

FIGURE 9

Figure 9. Effects of the nearest neighbor co-operative parameter γ on the time transients of (A) left ventricular pressure (LVP), (B) volume (LVV), (C) the force-pCa relationship under the unloaded sarcomere length, and (D) time transient of state ratio [XB_PostR2] in the logarithmic scale. The black, blue, and red lines indicate γ = 40, 20, and 10, respectively.

Sensitivity of the Power and Reverse Stroke Rate Constants

To examine the capability of MusAsi to reflect the stochastic behavior of the power and reverse strokes on cardiac outcomes, we examined the following alternative of the load dependent power stroke model, which we adopted in our previous work to reproduce the SPOC of a single myofibril of rabbit iliopsoas muscle (Washio et al., 2019):

$$\{\begin{array}{c}\hfill h_{f,i}\left(x\right)=g_i\text{exp}\left(\frac{E_{i-1}+W\left(x\right)-E_i-W\left(x+s_i/2\right)}{k_{B}T}\right)\hspace{1em}\hspace{1em} (38)\hfill \\ \hfill h_{b,i}(x+s_i)=g_i\text{exp}\left(\frac{W\left(x+s_i\right)-W\left(x+s_i/2\right)}{k_{B}T}\right).\hspace{1em}\hspace{1em}\hspace{0.5em}(39)\hfill \end{array}$$

This model originated from the Kramers escape theory (Kramers, 1940; Scherer and Fischer, 2017), in which the rate constants were defined by the Boltzmann factor associated with the height of the energy barrier from the origin, whereas the previous definition in Eqs 2, 3 used the strain energy at the destination. In Eqs 38, 39, W(x + s_i/2) was introduced to represent the contribution of the strain energy at the energy barrier that was assumed to be located at the mid strain (Washio et al., 2017). The contribution of the free energy of the myosin head at the barrier was included in the constant g_i. In this study, g₁ = 20 [1/s] and g₂ = 0.1 [1/s] were adopted, as in our previous work (Washio et al., 2019). Furthermore, the rate constant (k_np in Figure 2A) from the non-binding state N_XB to the weak binding state P_XB was multiplied by the factor 1.1 so that the same maximal pressure (P_max: the maximum of LVP) was achieved by both models. Hereafter, we call the power stroke models of Eqs 2, 3 and of Eqs 38, 39) the destination strain energy (DSE) model and barrier strain energy (BSE) model, respectively.

Both models reproduced similar tendencies in the force-pCa relationship (Figure 10A) and the force-velocity relationship (Figure 10B) though the active tensions of the BSE model were slightly smaller than that of the DSE model for a large Ca²⁺ concentration (>0.5 μM) or a small half-sarcomere shorting velocity (<1 μm/s). The SPOCs on the single myofibril model consisting of 40 half-sarcomeres under the constant [Ca^{2 +} ]=0.3μM were also reproduced by both models (Figures 10C,D). However, their periods and amplitudes were different (Figures 11A–F). In particular, the remarkable increases of two reverse stroke rates, which we called the avalanche of reverse strokes in our previous works (Washio et al., 2017, 2019), were observed for the two reverse rates at the lengthening in the BSE model, whereas such an increase was slightly recognized only in the reverse stroke rate from XB_PostR1 and XB_PreR in the DSE model. Furthermore, the ratios of the reverse rates to the forward rates were higher for the DSE model than for the BSE model. These differences in the SPOCs of the two models were also recognized in the numerical results for the ventricle FE model (Figures 12A–F). The local contraction duration of the DSE model was longer than that of the BSE model (Figures 12C,D) as the difference in the SPOC periods between the two models (Figures 11A,B). The rise and drop of LVP for the BSE model were slower than for the DSE model. In particular, for the BSE model, there was a small rebound of the binding population once after the binding myosin molecules almost disappeared. Although the rebound population was small, the binding myosin molecules clearly hampered the drop of LVP. In the BSE model, though the quick lengthening of sarcomere accompanying the avalanche of reverse strokes was observed, it was not reflected in the rapid transition from the systolic phase to the diastolic phase because of the rebounds.

FIGURE 10

Figure 10. Comparisons of the DSE and BSE models in (A) the force-pCa relationship under the unloaded sarcomere length, (B) the force-velocity relationship at [Ca²⁺] = 0.7 μM. Both the (C) DSE and (D) BSE models reproduced the SPOC for 20 sarcomeres in a single myofibril under the constant Ca²⁺ concentration ([Ca²⁺] = 0.3 μM).

FIGURE 11

Figure 11. (A,B) Time transients of the sarcomere length (SL), (C,D) ratio of binding states, and (E,F) transition rates between the binding states for the half-sarcomere model imbedded at the center of the myofibril model consisting of 20 sarcomeres during the SPOC with [Ca²⁺] = 0.3 μM. Panels (A,C,E) and panels (B,D,F) are the results of the DSE and BSE models, respectively. (E,F) Transition rates were calculated by dividing the total number of transitions per unit time by the total number of myosin molecules.

FIGURE 12

Figure 12. (A,B) Time transients of the left ventricular pressure (LVP) and volume (LVV) of the FE ventricle model, (C,D) the ratio of binding states to the Ca²⁺ concentration, and (E,F) the transition rates between the binding states and the stretch given by averaging the half-sarcomere models imbedded at the endo lateral left ventricular wall (elements inside a sphere of 3 mm radius). Panels (A,C,E) and panels (B,D,F) are the results of the DSE and BSE models, respectively. In panels (A,B), the broken lines indicate the results of counterparts. The transition rates in panels (E,F) were calculated by dividing the total number of transitions per unit time by the total number of myosin molecules.

In Figure 13, the distribution of binding states on the rod strain space x ∈ [−10nm,10nm] on the elements at the endo lateral left ventricular wall are provided with the rate constants. The large magnitude of the reverse stroke rates (h_b,1,h_b,2) of the BSE model (Figure 13B) caused an avalanche of reverse strokes and local quick lengthening (Figure 12F). The small distribution of the rebound after lengthening was recognized (Figure 13D). Note that the MusAsi scheme allowed us to directly couple such distributions of rod strains with wall motion in the biventricular FE model.

FIGURE 13

Figure 13. (A,B) Rate constants in the logarithmic scale between the binding states, and (C,D) time transients of the distribution of binding states in the rod strain space in the half-sarcomere models imbedded at the endo lateral left ventricular wall (elements inside a sphere of 3 mm radius). Red and blue indicate high and low density, respectively. Panels (A,C) and panels (B,D) are for the DSE and BSE models, respectively.

Discussion

Practicality in Clinical Applications

Using 320 cores of a conventional parallel computer system, one cardiac cycle of the MusAsi scheme could be computed within 1.5 h for a sufficiently fine FE biventricle model consisting of 45,000 tetrahedral elements. Accuracy and robustness were also confirmed through sensitivity tests with the various parameters of sample numbers and time step sizes. In fact, we have already applied the approach to follow-up verifications of practical clinical problems (Kariya et al., 2020; Masuda et al., 2021). In these cases, pumping performance after operations was predicted not only using the standard indices, such as LVP, ejection fraction (EF), and stroke volume (SV), but also the energy consumption. Now, we are moving to the next stage of applying our simulator using the MusAsi scheme in prospective clinical trials in an ongoing project on congenital heart disease.

Relaxation Mechanism

How the rapid drop of LVP can be achieved for the intracellular Ca²⁺ transient with a slow attenuation (Figure 2C) may be still a controversial issue. Furthermore, the population of binding myosin molecules must almost vanish in the relaxation phase, whereas nearly 10% of the maximum is left in the Ca²⁺ concentration (Figures 12C,D). In our model, the latter problem was resolved by adopting the nearest neighbor co-operative mechanism in transitions between the non-binding state and weak binding state (Figures 2A, 9), whereas the former problem was resolved by the reverse stroke mechanism introduced by the load-dependent power stroke model (Figure 12). The MusAsi scheme enabled the stochastic cross-bridge mechanisms and the macroscopic dynamics to be directly coupled, and we confirmed that these molecular mechanisms work efficiently to achieve the physiological relaxation of cardiac muscle in the beating cycle. A comparison of relaxation in the two power stroke models indicated the usefulness of the MusAsi scheme as a basic research tool in fields that study the role of molecular-level observation in the heartbeat (Figures 10–13). In our previous work (Washio et al., 2018) in which we directly coupled the Langevin dynamics model and the FE ventricle model, we detected the same problem of the slowed LVP drop caused by the rebounds observed in the BSE model. This problem was resolved by introducing the trapping mechanism that inhibits the reverse strokes when the rod strains increased quickly over a certain threshold. The trapping mechanism may have similar effects to the DSE model in which the reverse rates are drastically reduced for large strains (Figure 13A). Furthermore, the experimental measurements made by Hwang et al. (2021) revealed a higher frequency of backward steps at lower loads of the cardiac myofilaments than those of fast skeletal myofilaments like the higher reverse rate h_b,1 of the DSE model than that of the BSE model (Figures 13A,B). Our numerical results suggest that such a characteristic of the reverse rate brings the benefit to the cardiac myofilaments for quick relaxation. Note that achieving the quick relaxation of muscle is crucial in simulations of congenital heart disease because heart rates are more than a hundred, in most cases.

Significance of Active Stiffness

A key factor for stability in the MusAsi scheme is the implicit treatment of the active stress tensor S_act in the standard FE framework using Newton iterations. Assuming a constant stiffness k_rod per binding myosin molecule and a binding ratio R_B in the half-sarcomere model, the macroscopic axial stiffness coefficient K_A in the fiber orientation for the active tension in Eq. 8 is estimated as

$$K_A=\frac{R_S}{S{A}_0}{R}_{B}2N_M{k}_{rod}\frac{S{L}_0}{2}=104.2R_B\text{MPa}(40)$$

with the adopted parameters: sarcomere volume ratio R_S = 0.5, cross-sectional area SA₀ = 693nm² (Sato et al., 2013), number of accessible myosin molecules 2N_M = 76 per thin filament, spring coefficient of the myosin rod with positive strains k_rod = 2pN/nm, and HSL SL₀/2=0.95μm. Because we applied the viscosity coefficient μ_S = 36.66Pa⋅s, the stability condition for the time step size in the explicit approach is roughly estimated as

$$\mathrm{\Delta }T\le \frac{{\mathrm{\mu }}_S}{K_A}\sim \frac{0.35}{R_B}\mathrm{\mu }\mathrm{s}.(41)$$

Note that the stiffness coefficients for passive stress are in the order of kPa for strains less than 0.2 (see Supplementary Material 2.1 for details of the passive material parameters). Therefore, the contribution of passive stiffness, which is dealt with implicitly, is negligible compared with active stiffness in Eq. 40, even if only a few percent of myosin molecules are in binding states (R_B∼0.02). The above estimations of the limitation of the time step size in the explicit approach for active tension are good fit for the instability depending on the time step size in Figure 6. The stiffness estimation in Eq. 40 also justifies the significant influence of a single binding myosin molecule contained in the half-sarcomere model (R_B∼1/38) regarding hampering the diastole as observed in Figure 9.

Future Directions

In the one-dimensional (1D) half-sarcomere model adopted in this study, the characteristics of the realistic 3D regular arrangements of myosin molecules on the thick filament and the binding sites on the double spirals on the thin filament (Hussan et al., 2006) were not taken into account. These geometrical parameters are likely to have been optimized in the process of evolution. Thus, they may have a significant influence on the rate constant of binding transitions from the P_XB state to the XB_PreR state and the initial rod strains that were provided probabilistically from the Boltzmann distribution exp(−W(x)/k_BT) of strain energy W in this study. The modeling of active stress was also limited only to the fiber orientations in this study. However, actin filaments are pulled not only in the longitudinal direction of the sarcomere but also in the lateral direction by myosin rods. These limitations should be removed by extending the MusAsi scheme from the 1D model to an appropriate 3D model in our future work.

In this study, we applied the MusAsi scheme to the simplified model to focus on the impact of the properties of contractile proteins on the macroscopic outcomes. In the heart model of our previous work (Kariya et al., 2020) in which a similar approach for coupling MC and FEM simulations (Washio et al., 2016) has been used, three species of ventricular myocytes, i.e., endocardial, mid-myocardial, and epicardial cells, were implemented. Therefore, we believe that the current scheme will also work in a model implemented with the realistic electrophysiology. In that work, we also assumed that the heart walls were surrounded by the pericardium, which was fixed in space by the planar springs, and we incorporated the impact of pericardial pressure that was generated based on volume conservation of pericardial liquid. It is expected that the negative pericardial pressure also facilitates the drop of LVP at the early diastole. The contributions of these more realistic boundary conditions will be evaluated in future studies that should also take the pre-stress of the myocardium into account.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

TH, SS, and MW designed the project. J-IO and KY prepared the input data for the computer simulations. TW and KY developed the multiple time step active stiffness scheme, ran the simulations, analyzed the simulation data, and wrote the manuscript with input from SS, TH, and MW. TW, J-IO, and KY developed the simulation code with input from SS, TH, and MW. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by MEXT under the “Program for Promoting Researches on the Supercomputer Fugaku” (hp200121) and used the computational resources of Supercomputer Fugaku provided by the RIKEN Center for Computational Science, and by AMED under Grant Number JP21he2102003.

Conflict of Interest

TW, J-IO, SS, and TH were employed by UT-Heart Inc. KY and MW were employed by Fujitsu Japan Ltd.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Maxine Garcia, from Edanz (https://jp.edanz.com/ac) edited 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.2021.712816/full#supplementary-material

References

Azzolin, L., Dedé, L., Gerbi, A., and Quarteroni, A. (2020). Effect of fibre orientation and bulk modulus on the electromechanical modelling of human ventricles. Math. Eng. 2, 614–638. doi: 10.3934/mine.2020028

CrossRef Full Text | Google Scholar

Brezzi, F., and Fortin, M. (1991). Mixed and Hybrid Finite Element Methods. Springer Series in Computational Mathematics, Vol. 15. Berlin: Springer, 200–273.

Google Scholar

Caruel, M., Moireau, P., and Chapelle, D. (2019). Stochastic modeling of chemical–mechanical coupling in striated muscles. Biomech. Model. Mechanobiol. 18, 563–587. doi: 10.1007/s10237-018-1102-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Chapelle, D., Le Tallec, P., Moireau, P., and Sorine, M. (2012). Energy-preserving muscle tissue model: formulation and compatible discretizations. Int. J. Multiscale Comput. Eng. 10, 189–211. doi: 10.1615/IntJMultCompEng.2011002360

CrossRef Full Text | Google Scholar

Dabiri, Y., Van der Velden, A., Sack, K. L., Choy, J. S., Kassab, G. S., and Guccione, J. M. (2019). Prediction of left ventricular mechanics using machine learning. Front. Phys. 7:117. doi: 10.3389/fphy.2019.00117

PubMed Abstract | CrossRef Full Text | Google Scholar

Guérin, T., Prost, J., and Joanny, J. F. (2011). Dynamical behavior of molecular motor assemblies in the rigid and crossbridge models. Eur. Phys. J. E Soft. Matter. 34, 60. doi: 10.1140/epje/i2011-11060-5

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunter, P., McCulloch, A., and Ter Keurs, H. (1998). Modelling the mechanical properties of cardiac muscle. Prog. Biophys. Mol. Biol. 69, 289–331. doi: 10.1016/S0079-6107(98)00013-3

CrossRef Full Text | Google Scholar

Hussan, J., de Tombe, P. P., and Rice, J. J. (2006). A spatially detailed myofilament model as a basis for large-scale biological simulations. IBM J. Res. Dev. 50, 583–600. doi: 10.1147/rd.506.0583

PubMed Abstract | CrossRef Full Text | Google Scholar

Hwang, Y., Washio, T., Hisada, T., Higuchi, H., and Kaya, M. (2021). A reverse stroke characterizes the force generation of cardiac myofilaments, leading to an understanding of heart function. Proc. Natl. Acad. Sci. U.S.A. 118:e2011659118. doi: 10.1073/pnas.2011659118

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishiwata, S., Shimamoto, Y., and Fukuda, N. (2010). 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

Kagemoto, T., Li, A., dos Remedios, C., and Ishiwata, S. (2015). Spontaneous oscillatory contraction (SPOC) in cardiomyocytes. Biophys. Rev. 7, 15–24. doi: 10.1007/s12551-015-0165-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kariya, T., Washio, T., Okada, J., Nakagawa, M., Watanabe, M., Kadooka, Y., et al. (2020). Personalized perioperative multi-scale, multi-physics heart simulation of double outlet right ventricle. Ann. Biomed. Eng. 48, 1740–1750. doi: 10.1007/s10439-020-02488-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerckhoffs, R. C. P., 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

Kramers, H. A. (1940). Brownian motion in a field of force and the diffusion model of chemical reactions. Physica 7, 284–304. doi: 10.1016/S0031-8914(40)90098-2

CrossRef Full Text | Google Scholar

Masuda, H., Miyagawa, S., Sugiura, S., Washio, T., Okada, J., Ueno, T., et al. (2021). An application of a patient-speciic cardiac simulator for the prediction of outcomes after mitral valve replacement: a pilot study. J. Artif. Organs doi: 10.1007/s10047-021-01248-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Negroni, J. A., and Lascano, E. C. (2008). Simulation of steady state and transient cardiac muscle response experiments with a Huxley-based contraction model. J. Mol. Cell Cardiol. 45, 300–312. doi: 10.1016/j.yjmcc.2008.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S. A., Hunter, P. J., and Smith, N. P. (2006). A quantitative analysis of cardiac myocyte relaxation: a simulation study. Biophys. J. 90, 1697–1722. doi: 10.1529/biophysj.105.069534

PubMed Abstract | CrossRef Full Text | Google Scholar

Okada, J., Sasaki, T., Washio, T., Yamashita, H., Kariya, T., Imai, Y., et al. (2013). Patient specific simulation of body surface ECG using the finite element method. Pacing. Clin. Electrophysiol. 36, 309–321. doi: 10.1111/pace.12057

PubMed Abstract | CrossRef Full Text | Google Scholar

Regazzoni, F., Dedé, L., and Quarteroni, A. (2018). Active contraction of cardiac cells: a reduced model for sarcomere dynamics with cooperative interactions. Biomech. Model. Mechanobiol. 17, 1663–1686. doi: 10.1007/s10237-018-1049-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Regazzoni, F., Dedé, L., and Quarteroni, A. (2020). Biophysically detailed mathematical models of multiscale cardiac active mechanics. PLoS Comput. Biol. 16:e1008294. doi: 10.1371/journal.pcbi.1008294

PubMed Abstract | CrossRef Full Text | Google Scholar

Regazzoni, F., Dedé, L., and Quarteroni, A. (2021). Active force generation in cardiac muscle cells: mathematical modeling and numerical simulation of the actin-myosin interaction. Vietnam J. Math. 49, 87–118. doi: 10.1007/s10013-020-00433-z

CrossRef Full Text | Google Scholar

Regazzoni, F., and Quarteroni, A. (2021). An oscillation-free fully staggered algorithm for velocity-dependent active models of cardiac mechanics. Comp. Methods Appl. Mech. Eng. 373:113506. doi: 10.1016/j.cma.2020.113506

CrossRef Full Text | Google Scholar

Rice, J. J., Stolovitzky, G., Tu, Y., and de Tombe, P. P. (2003). Ising model of cardiac thin filament activation with nearest-neighbor cooperative interactions. Biophys. J. 84(2 Pt 1), 897–909. doi: 10.1016/S0006-3495(03)74907-8

CrossRef Full Text | Google Scholar

Rice, J. J., Wang, F., Bers, D. M., and de Tombe, P. P. (2008). Approximate model of cooperative activation and crossbridge cycling in cardiac muscle using ordinary differential equations. Biophys. J. 95, 2368–2390. doi: 10.1529/biophysj.107.119487

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, K., Kuramoto, Y., Ohtaki, M., Shimamoto, Y., and Ishiwata, S. (2013). Locally and globally coupled oscillators in muscle. Phys. Rev. Lett. 111:108104. doi: 10.1103/PhysRevLett.111.108104

PubMed Abstract | CrossRef Full Text | Google Scholar

Scherer, P. O. J., and Fischer, S. F. (2017). Theoretical Molecular Biophysics. Berlin: Springer.

Google Scholar

Shavik, S. M., Wall, S. T., Sundnes, J., Burkhoff, D., and Lee, L. C. (2017). Organ-level validation of a cross-bridge cycling descriptor in a left ventricular finite element model: effects of ventricular loading on myocardial strains. Physiol. Rep. 5:e13392. doi: 10.14814/phy2.13392

PubMed Abstract | CrossRef Full Text | Google Scholar

Shintani, S. A., Washio, T., and Higuchi, H. (2020). Mechanism of contraction rhythm homeostasis for hyperthermal sarcomeric oscillations of neonatal cardiomyocytes. Sci. Rep. 10:20468. doi: 10.1038/s41598-020-77443-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, N. P., Nickerson, D. P., Crampin, E., and Hunter, P. (2004). Multiscale computational modelling of the heart. Acta Numerica 13, 371–431. doi: 10.1017/S0962492904000200

CrossRef Full Text | Google Scholar

Syomin, F. A., and Tsaturyan, A. K. (2017). A simple model of cardiac muscle for multiscale simulation: passive mechanics, crossbridge kinetics and calcium regulation. J. Theor. Biol. 420, 105–116. doi: 10.1016/j.jtbi.2017.02.021

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., 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., and Hisada, T. (2011). Convergence analysis of inexact LU-type preconditioners for indefinite problems arising in incompressible continuum analysis. Jpn. J. Indust. Appl. Math. 28, 89–117. doi: 10.1007/s13160-011-0024-2

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. Phys. Rev. Eng. 95:022411. doi: 10.1103/PhysRevE.95.022411

PubMed Abstract | CrossRef Full Text | Google Scholar

Washio, T., Okada, J., Sugiura, S., and Hisada, T. (2012). Approximation for cooperative interactions of a spatially-detailed cardiac sarcomere model. Cell Mol. Bioeng. 5, 113–126. doi: 10.1007/s12195-011-0219-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Washio, T., Okada, J., Takahashi, A., Yoneda, K., Kadooka, Y., Sugiura, S., et al. (2013). Multiscale heart simulation with cooperative stochastic cross-bridge dynamics and cellular structures. Multiscale Model. Simul. 11, 965–999. doi: 10.1137/120892866

CrossRef Full Text | Google Scholar

Washio, T., Shintani, S. A., Higuchi, H., and Hisada, T. (2019). Effect of myofibril passive elastic properties on the mechanical communication between motor proteins on adjacent sarcomeres. Sci. Rep. 9:9355. doi: 10.1038/s41598-019-45772-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Washio, T., Sugiura, S., Kanada, R., Okada, J., and Hisada, T. (2018). Coupling Langevin dynamics with continuum mechanics: exposing the role of sarcomere stretch activation mechanisms to cardiac function. Front. Physiol. 9:333. doi: 10.3389/fphys.2018.00333

PubMed Abstract | CrossRef Full Text | Google Scholar

Washio, T., Sugiura, S., Okada, J., and Hisada, T. (2020). Using systolic local mechanical load to predict fiber orientation in ventricles. Front. Physiol. 11:467. doi: 10.3389/fphys.2020.00467

PubMed Abstract | 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