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.

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 . 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.
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 crossbridge 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 crossbridge 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 , 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 2+ concentrations, as  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.

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λ T in the fiber orientation of the FE model are used as the initial half-sarcomere length (HSL) λ T ·SL 0 /2 at T and its shorting velocity −λ T ·SL 0 /2 in the time interval [T, T T] of the half-sarcomere model (Figure 3A), where SL 0 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.
In our previous works (Washio et al., , 2016, the active tension T act,[T,T+ T] was implicitly computed by assuming the shorting velocity −λ T+ T ·SL 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 0 /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 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 −SL 0λ /2, where the stretch rateλ 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. 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 crossbridge 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 2+ -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 2+ 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 2+ concentration remains ( Figure 2C). Contraction force is generated by the power stroke transitions in which the strain of the myosin rod increases by s 1 and s 2 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: 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 .
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 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λ T . A simple approach to perform the MC steps with the time step size t = T/n is to define the stretch λ at the k-th step as In this case, the rod strainx 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: x ij,T+k t = x A,ij,T+k t +s ij,T+k t + SL 0 2 λ T+k t − λ A,ij,T+k t (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 Frontiers in Physiology | www.frontiersin.org 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 λ 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 B T . The variable s ij = 0, s 1 or s 1 + s 2 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 0 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 straiñ x ij,T+k 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: where SA 0 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 t ≡x ij,T+k 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 halfsarcomere 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 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: λ A,ij,T+k t λ T+k T , δ A,ij,T+κ t = 1 and δ A,ij,T+(k−1) t = 0 λ A,ij,T+(k−1) t , otherwise.
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 where the derivative with respect to λ T+ T is given by and starting from Therefore, the stiffness ∂T act,[T,T+ T] /∂λ T+ T is always non-negative, provided the potential W is convex downward (d 2 W/dx 2 ≥ 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 From the relationship between the stretch λ and the Green- Thus, the infinitesimal increment of stretch is represented by Therefore, if we define the second Kirchhoff active stress tensor as the infinitesimal work is represented by δW act = S act : δ E.
The derivative of the infinitesimal work is given by 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 with the normal fiber orientation vector in the current coordinate a = Ff /λ. Therefore, δ 2 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 < 0 may 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 ∼ 1 ms, in our experience. When some impaired crossbridge 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 λ T+k t ≡ λ T + k tλ T+ T was used instead of Eq. 10, the first term in Eq. 21 was replaced with ∂T act /∂λ δλ∂λ. 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): where J = detF is the Jacobian, p is the hydrostatic pressure in the ventricular walls, C = F T F 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: 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: where the vector U contains all variables as follows: 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: U Compute the residual and the stiffness matrix: Solve the linear system: Because the initial guesses U T+ T 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 where the derivatives of stretch λ with respect to nodal displacements u e of the element e are given by Eqs 19, 22. In 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 derivativeu 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.
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.

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 . 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.25 ms, 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 2+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 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).

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 SV: stroke volume, EDV: end-diastolic volume, EDP: end-diastolic pressure, P max : maximal pressure, W out : work at output, ATP: ATP consumption, time: computational time for one cycle. These values were taken from the third cycle after the initial process of blood filling. ATP consumption in the septum was included. T = 1.25 ms and t = 5 µs were applied in (A). N F = 16 and T 0 = 1.25 ms were applied in (B).
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. 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 0 = 5 µs, and the number of MC steps n performed within a single FE step was determined by where " " represents the floor function that rounds down after the decimal point. Therefore 250, 125, and 63 MC steps were performed with T = T 0 , T 0 /2, and T 0 /4 ( T 0 = 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 x ij,T+k t in Eq. 7 calculated in the MC steps from the stretch λ and stretch rateλ 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λ 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 0 =1.25 ms was reproduced with quite a small time step T = T 0 /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 0 /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.

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).
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 2+ -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).

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): 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 1 = 20 [1/s] and g 2 = 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 nonbinding 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 2+ 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(Washio et al., , 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.
In Figure 13, the distribution of binding states on the rod strain space x ∈ [−10 nm, 10 nm] 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.

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 2+ 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 2+ 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 with the adopted parameters: sarcomere volume ratio R S = 0.5, cross-sectional area SA 0 = 693 nm 2 (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 = 2 pN/nm, and HSL SL 0 /2 = 0.95 µm. Because we applied the viscosity coefficient µ S = 36.66 Pa·s, the stability condition for the time step size in the explicit approach is roughly estimated as 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 B T 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, midmyocardial, 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.