Three parameters, s trap , h trap , and h escape , characterize the trapping. We adjusted these parameters so that the experimental stretch activation results reported by Stelzer et al. (2006) were reproduced. The adjusted parameter values are s trap = 1.3 nm, h trap = 50 1/s, and h escape = 5,000 1/s. From Eq. 5, we find that Δ G trap − Δ G 1 = 4.6 k B T . According to Eq. 6, the rate constants of trapping r trap and escaping r escape are equal at a distortion of x = 4.8 nm in the XB_Trap state and x + s trap = 6.1 nm in the XB_PostR1 state, assuming a spring constant of k xb = 2.8 pN/nm (Kaya and Higuchi, 2010) and a physiological body temperature of T = 310 K.

In our model, we assume that attachment, which represents the transition from the P_XB state to the XB_PreR state (Figure 1C), is allowed only in the single overlap region of the thin and thick filaments. The myosin head (MH) (# i ) is situated in the single overlapping region only if the following condition is fulfilled: max ( L A − S L 2 , S L 2 − L A ) ≤ L B 2 + 2 ( i − 0.5 ) n M ( L M − L B ) ≤ S L 2   . (7)

Here, the middle term is the distance from the center of the sarcomere .   L M , L B , and L A represent the lengths of the thick filament, bare zone, and thin filament, respectively (Figure 1B). S L is the sarcomere length. The parameters for the sarcomere geometry were determined from cardiac sarcomeres (Rodriguez et al., 1993; Lodish et al., 2000; Rice et al., 2008; Kolb et al., 2016).

Two states (Ca-off and Ca-on) are assumed by each T/T unit (Figure 1D). The transitions between the states of the T/T unit are determined by the Ca²⁺ concentration, [ Ca 2 + ] , and the two parameters K ¯ on and K ¯ off . Those two parameters are defined as follows: K ¯ on = { K on ′             i f   t h e r e   i s   a n   M H   b i n d i n g   b e l o w , K on   o t h e r w i s e . K ¯ off = { K off ′                 i f   t h e r e   i s   a n   M H   b i n d i n g   b e l o w , K off       o t h e r w i s e .

The transitions between the N_XB and P_XB states (Figure 1C) are affected by the status of the T/T unit above, via modifications of k np and k pn , as well as by the states of the neighboring MHs through the integer n g . The value of n g   ( = 0 ,   1 ,   o r   2 ) represents the number of neighboring MHs in the P_XB state or the four attached states. The corresponding T/T unit index τ for the i -th MH is given by: τ = ⌊ 0.5 L B + ( i − 0.5 ) S M − ( S L / 2 − L A ) S T ⌋ . (8)

Here, ⌊ ⌋ indicates the floor function, which rounds down after the decimal point. The parameter S M = 0.5 ( L M − L B ) / N M represents the spacing of the MHs, and S T = L A / N T is the spacing of the T/T units. The corresponding T/T unit exists only if 1 ≤ τ ≤ N T . Based on this correspondence, the factors k np and k pn of the rate constants are given by: k n p = { δ O V K np 1           i f   t h e   T / T   u n i t   a b o v e   i s   i n   t h e   C a − o n   s t a t e , δ O V K np 0           o t h e r w i s e . (9) k p n = { K p n 1             i f   t h e   T / T   u n i t   a b o v e   i s   i n   t h e   C a − o n   s t a t e , K pn 0           o t h e r w i s e . (10)

Here, δ O V = 1 if the MH is located in the single overlapping region with the thin filament; otherwise, δ O V = 0 . This, along with γ n g or γ − n g ( γ = 80 ) , represents the nearest-neighbor cooperativity of the MHs, following Rice et al. (Rice et al., 2003), which plays an important role in the force-pCa relationship, as shown in Figures 5A,B. From the cooperativity with γ = 80 , the cardiac muscle in our ventricle model can relax in the diastole when [Ca²⁺] is lower than 0.1 μM (Figure 4C), whereas the diastole deteriorates slightly with γ = 40 (Figures 5C,D). The times to reach equilibria differ substantially and are dependent on the Ca²⁺ concentration and parameter γ . Thus, cooperativity also affects the maximal rising rate of the left ventricular pressure (dP/dt_max), which is an important contractility index used in the medical community (Figures 5C,D).

We assume that one thin filament in the three-dimensional arrangement corresponds to two thin filaments in our half-sarcomere model. This case arises because we assumed that cooperative behavior exists along the tropomyosin molecules wrapped around the thin filament in a double-spiral fashion, and one spiral is modeled in our half-sarcomere. The constants K np 0 , K np 1 , K pn 0 , and   K pn 1 are determined from Q , K basic , and μ , as follows: K n p 0 = F K ( S L ) Q K b a s i c μ ,   K n p 1 = F K ( S L ) Q K b a s i c ,    K p n 0 = K p n 1 = K b a s i c γ 2 . (11)

Here, μ > 1 controls the degree of cross-bridge inhibition for the T/T units in states other than Ca-on, and Q controls the ratio of binding states for the MHs. To reproduce the SL dependence in the active contraction tension, the following function F K ( S L ) is multiplied to define K np 0 and K np 1 . F K ( S L ) = { 1 ,                                                                             S L ≥ S L Q , 1 − α Q ( S L Q − S L ) ,         S L < S L Q . (12)

The values α Q = 0.25 [1/ μ m] and S L Q = 2.1  μ m are used in this study.

The rate constants of attachment r a   and detachment r d   between the P_XB state and XB_PreR state are based on the assumed free energies, Δ G w b   and Δ G 0 , of the P_XB state and XB_PreR state, respectively: r a = c pre exp ( − Δ G 0 − Δ G w b k B T ) ,       r d = c pre . (13)

The initial rod distortion at attachment is given stochastically based on a Boltzmann distribution determined from the rod strain energy, W r o d ( x ) (Washio et al., 2016). The detachment rate constant for transitions from the XB_PostR2 state to the N_XB state is r NXB .

Using the rate constants, we also considered forced detachments from the XB_PreR, XB_PostR1, and XB_PostR2 states to the N_XB state caused by extreme strain on the myosin rod: d X B , i ( x ) = { 0 ,                                                                                                       x < x X B , i , c X B , i ( exp ( a X B , i ( x − x X B , i ) ) − 1 ) , x ≥ x X B , i ,               i = 0 ,   1 ,   2. (14)

Similarly, the rate constant of detachment from the XB_Trap state is given by: d trap ( x ) = { 0 ,                                                                                                 x < x X B , t , c X B , t ( exp ( a X B , t ( x − x X B , t ) ) − 1 ) ,           x ≥ x X B , t . (15)

We paid close attention to the adjustments of parameters in Eqs 14 and 15 in reproducing the stretch activation. These parameters affect the degree of the increased delayed force. In this study, we used the following parameters to reproduce the characteristics of cardiac muscles: c X B , 0 = 100 1/s, c X B , 1 = c X B , 2 = c X B , t = 20 1/s, a X B , 0 = 1 1/nm, a X B , 1 = a X B , 2 = a X B , t = 0.1 1/nm, x X B , 0 = 3 nm, x X B , 1 = 8.5 nm, x X B , 2 = 13 nm, x X B , t = 8.5 nm.

In our model, the rate constants of the power and reverse strokes are determined from the Kramers escape theory (Kramers 1940), in which rate constants are defined by the Boltzmann factor associated with the height of the energy barrier from the origin. r ^ f , i ( x ) = h i e x p ( Δ G i − 1 + W r o d ( x ) − W r o d ( x + s i / 2 ) k B T ) , (16) r ^ b , i ( x + s i ) = h i e x p ( Δ G i + W r o d ( x + s i ) − W r o d ( x + s i / 2 ) k B T ) . (17)

Here, W r o d ( x + s i / 2 ) is introduced because we assume that the barrier is given at the middle of the power stroke distances for the two states. In addition, we adopt upper limits to those transition rate constants because the power stroke transitions are thought to accompany a release of Pi and ADP from the nucleotide-binding pocket in the MH during the first and second strokes, respectively. Thus, the reverse stroke transitions must accompany a rebinding of these molecules. We assumed that the rates of chemical reactions for the release and rebinding are given by r ¯ f , i and r ¯ b , i , respectively. With these upper limits, the temporary rate constants are modified so that the Boltzmann equilibrium in Eq. 1 is preserved:   r f , i ( x ) = { r ¯ f , i ,                                           x ≤ x ¯ f , i r ^ f , i ( x ) ,       x ¯ f , i < x ≤   x ¯ b , i r ^ f , i ( x ) r ¯ b , i r ^ b , i ( x + s i ) ,                 x > x ¯ b , i , (18) r b , i ( x + s i ) = { r ^ b , i ( x + s i ) r ¯ f , i r ^ f , i ( x ) ,                                   x ≤ x ¯ f , i r ^ b , i ( x + s i ) ,     x ¯ f , i < x ≤   x ¯ b , i r ¯ b , i ,                                                                                     x > x ¯ b , i . (19)

Here, x ¯ f , i and x ¯ b , i are the distortions at which the temporary rates reach the upper limits ( r ^ f , i ( x ¯ f , i ) = r ¯ f , i ,   r ^ b , i ( x ¯ b , i + s i ) = r ¯ b , i   ).

The elastic force of a myosin rod is nonlinear with respect to the strain (Kaya and Higuchi, 2010). We assume that the myosin rod behaves as a linear spring for positive stretches, whereas nonlinear behavior is introduced for negative stretches because of the slack induced along the myosin rod. The strain energy W r o d is found by integrating the force F rod from x = 0 , defined by: F r o d ( x ) = { b x b k x b ( x + ξ 1 ) − F 1 ,               x < − ξ 1 , k x b a x b ( exp ( a x b x ) − 1 ) ,               − ξ 1 ≤ x < 0 , k x b x ,               0 ≤ x . (20)

Here, a x b   and F 1 are determined from the other parameters, so that the function F r o d and its first derivative are continuous at ξ = 0 and − ξ 1 : { a x b = − ( ln ⁡ b x b ) ξ 1 ,                                         F 1 = k x b ( 1 − exp ( − a x b ξ 1 ) ) a x b . (21)

We constructed a half-sarcomere model (Figure 5A) with n F one-dimensional filament pairs in which the active tension generated by the bound myosin molecules is given as: T a c t = 2 R S S A 0 ⋅ n F ∑ β = 1 n F ∑ α = 1 n M δ A , α , β d W r o d ( x α , β ) d x   . (22)

Here, δ A , α , β = 1 if the MH is in an attached state; otherwise, δ A , α , β = 0 . The parameter S A 0 is the cross-sectional area per thin filament in an unloaded half-sarcomere (Sato et al., 2013). The factor of 2 comes from the fact that our one-dimensional model corresponds to one of the double spirals of actin monomers along the thin filament. The parameter R S is the sarcomere volume ratio under the unloaded condition. R S = 0.5 is used, which indicates that 50% of the total volume is occupied by the sarcomere in the cardiac muscle.

For the feedback between sarcomere dynamics and actomyosin dynamics, the relation between the time transients of the molecular and sarcomere variables can be expressed by:   t x = t   A x − (   t z −   t   A z ) . (23)

Here, t A is the most recent time at which the myosin molecule was attached, x   t A is the initial distortion at the attachment, and z   t is the half-sarcomere shortening length from the unloaded condition:   t z = − S L 0 2 (   t λ − 1 ) , (24) where S L 0 is the unloaded SL. Here, the stretch λ   t is the parameter given by macroscale models. The initial distortion is given from the Boltzmann distribution, assuming a fluctuation in the MH position under the potential W r o d (Washio et al., 2016).

The half-sarcomere shortening ( z ˙   t > 0 ) reduces the rod distortion (Eq. 23), which facilitates the power stroke (Eq. 16) and the deterioration of the reverse stroke (Eq. 17). Therefore, both the active tension and energy consumption are affected by the half-sarcomere shortening velocity. In Figure 6A, changes in the half-sarcomere length with various isotonic tensions computed using the MC cross-bridge model are depicted. In these isotonic contractions, the shortening velocities were computed over the common range of the half-sarcomere length (1.09–1.12 μm) to derive the tension-shortening velocity relationship (Figure 6B). By counting the number of detachments from the XB_PostR2 state (Figure 1C: r N X B ), the energy consumption rates were evaluated assuming the detachment requires the energy G A T P = 76.5   pN ⋅ nm , which is used in power stroke transitions (Figure 6C), as represented by Eq. 1. In the MC cross-bridge model, we assumed that 30% of G A T P is consumed during the first power stroke (6 nm), and the remaining G A T P is consumed during the second power stroke (4 nm) (See Table 1). The efficiencies were computed by dividing the work rates by the energy consumption rates (Figure 6D).

The energy loss in the cross-bridge cycle is given during the power stroke transitions and detachment from the XB_PostR2 state. The difference between Δ G i − Δ G i + 1 and W r o d ( x + s i ) − W r o d ( x ) is lost in the former, whereas the rod strain energy, W r o d ( x ) , at the detachment is lost in the latter. To observe the influence of the stiffness parameter k x b in Eq. 20 on the dynamics in isotonic shortening, the cross-bridge performance was evaluated with the stiffness parameter k x b = 1.4 pN/nm, which is half of our parameter k x b = 2.8 pN/nm found by Kaya and Higuchi, 2010. In this case, the phase changes of the shortening velocity were not reproduced (Figure 6A), and the changes in the energy consumption rate for the tension is gentler than our model because the rate constants of power stoke transitions are less sensitive to the tension than the case of k x b = 2.8 pN/nm. As a result, the peak efficiency is approximately half for our case (Figure 6D). The strain energy after the second power stroke ( 0.5 ( s 1 + s s ) 2 k x b ) with k x b = 1.4 pN/nm is 70 pN ⋅ nm assuming isometric contraction and zero distortion at XB_PreR. This energy is almost equal to G A T P . However, rod distortion decreases by half-sarcomere shortening until reaching XB_PostR2 (Eq. 23) during isotonic contraction. Therefore, 75% of the energy loss is accounted for even at the optimal shortening velocity.

In this study, coupling of the half-sarcomere model (Figure 4A) and the myofibril model (Figure 5B) or the ventricle model (Figure 4C) was achieved by applying the multiple-step active stiffness integration scheme in the exchange of the active tension T a c t and stretch λ between the two scales. The macroscale model is driven by the active stress given by the active tension in the half-sarcomere model, whereas the stretch in the fiber orientation in the macroscale model provides feedback to the length change of the half-sarcomere model.

In the macroscale model, a much larger time step Δ T = n Δ t compared with the MC time step Δ t is used to save computational time. For a given stretch   T + Δ T λ at T + Δ T , the macroscopic active tension T act , [ T , T + Δ T ] over the time interval [ T , T + Δ T ] is determined by taking the time average of the active tension in the half-sarcomere model given for each MC step within an interval of Δ t . T act , [ T , T + Δ T ] ( λ   T + Δ T ) = 1 n ∑ k = 1 n   T act ( T + k n Δ T , λ T + k n (   T + Δ T λ − λ T ) ) . (25)

The active tension T act in the half-sarcomere model is given as a function of time and stretch for each MC step. The stretch for each MC step is determined by interpolating the stretches at T and T + Δ T to evaluate an appropriate stiffness of the binding myosin molecules in macroscopic Newton iterations. Note that the state transitions in the MC computations are calculated only once before the Newton iterations, assuming a stretch of λ   T + Δ T = λ T + k Δ t λ T for each MC step. However, the rod distortion λ   T + Δ T is re-evaluated from Eqs 23 and 24 by replacing λ ˙   T with ( λ   T + Δ T − T λ ) / Δ T using the updated stretch λ   T + K Δ T in the macroscopic Newton iterations. In this study, Δ t = 2.5 μs is used for the MC step, while Δ T = 0.1 ms and Δ T = 1.25 ms are used, respectively, for the myofibril and ventricle models.

For simplicity, we introduce the Newton iteration for the one half-sarcomere model under the isotonic condition, as in Figure 6A. γ L λ ˙ + d φ L d λ ( λ ) + T a c t − T i s o = 0 , (26) where γ L = 0.01  kPa ⋅ s is the sarcomere friction, φ L is the deformation energy for stretching in the longitudinal direction (Washio et al., 2019), and T i s o is the isotonic tension imposed on the cardiac muscle. Eq. 26 is discretized as: γ L   T + Δ T λ − T λ Δ T + d φ L d λ (   T + Δ T λ ) + T act , [ T , T + Δ T ] (   T + Δ T λ ) − T i s o = 0. (27)

After executing n -times MC steps in [ T , T + Δ T ] , where the half-sarcomere shortening length z   T + K Δ T at the kth step is given by assuming λ   T + K Δ T = T λ + k Δ t T λ ˙ with   T λ ˙ = ( T λ − T − Δ T λ ) / Δ T , the initial solution of the Newton iteration is set as λ   T + K Δ T = λ T , and the following linearized equation of Eq. 27 at   T + Δ T λ = T + Δ T λ ( i ) is iteratively solved until the magnitude of the residual R ( i ) is reduced sufficiently. ( γ L Δ T + d 2 φ L d λ 2 ( T + Δ T λ ( i ) ) + K act , [ T , T + Δ T ] ( T + Δ T λ ( i ) ) ) ( T + Δ T λ ( i + 1 ) − T + Δ T λ ( i ) ) = R ( i ) , (28) where the residual is given by R ( i ) = − ( γ L   T + Δ T λ ( i ) − λ   T Δ T + d φ L d λ (   T + Δ T λ ( i ) ) + T act , [ T , T + Δ T ] (   T + Δ T λ ( i ) ) − T i s o ) . (29)

Here, K act , [ T , T + Δ T ] (   T + Δ T λ ) is the total stiffness of the binding myosin rods given as: K act , [ T , T + Δ T ] (   T + Δ T λ ) = d T act , [ T , T + Δ T ] d λ (   T + Δ T λ ) . (30)

See our previous work (Yoneda et al., 2021) for the actual computation of the right-hand side in Eq. 30. Under normal situations, the stiffness coefficient K act , [ T , T + Δ T ] increases greatly to be much larger than γ L / Δ T + d 2 φ L / d λ 2 during contraction. Therefore, stiffness caused by the binding myosin rods must be correctly incorporated in the total stiffness, as in Eq. 28; otherwise, the time integration scheme generates inaccurate solutions without using a small time step Δ T in a similar order of magnitude to the MC time step Δ t (Yoneda et al., 2021).

The above implicit time integration scheme can be naturally extended to more general cases where the half-sarcomere models are imbedded in different elements that interact with each other at the element interfaces. In the myofibril model (Figure 4B), two degrees of freedom, the stretches λ i and μ i , respectively, in the longitudinal and transverse directions, are given to each half-sarcomere element (Washio et al., 2019). For each half-sarcomere, the passive deformation energy per unit volume is given by ψ ( λ , μ ) = φ L ( λ ) + φ T ( μ ) + 1 2 k L T R L T ( λ , μ ) 2 . (31)

Here, the function φ T is the deformation energy for the transverse stretch, and the last term is a weak penalty term associated with the inverse SL-lattice space (LS) relationship, which constrains the half-sarcomere deformation to make | R L T ( λ , μ ) | small. In our model, a simple linear relation: R L T ( λ , μ ) = λ − 1 + 2 β R ( μ − 1 ) , (32) with β R = 2 is applied (Washio et al., 2019). Concerning the elasticity with differences in LS between adjacent half-sarcomeres, the following deformation energy per unit volume is further applied. φ T A , i ( μ i , μ i + 1 ) = { k M M ( μ i − μ i + 1 ) 2 ,   i = 1,3 , ⋯ k M Z ( μ i − μ i + 1 ) 2 , i = 2,4 , ⋯ (33) where the half-sarcomeres are separated by the M-band for the top expression, and are separated by the Z-disc for the lower expression.

Within each half-sarcomere, the following tensions act at the left and right edges: T S R , i = γ L λ ˙ i + d φ L d λ ( λ i ) + k L T R L T ( λ i , μ i ) + T a c t , i   ,   i = 1 , ⋯ , n s . (34)

The longitudinal mechanical equilibrium condition at the element boundaries: T S R , i − T S R , i + 1 = 0 , (35) and the transversal mechanical equilibrium condition at each element: γ T μ ˙ i + d φ T d μ ( μ i ) + 2 β R k L T R L T ( λ i , μ i ) + ∂ φ T A , i − 1 ∂ μ i ( μ i , μ i − 1 ) + ∂ φ T A , i ∂ μ i ( μ i , μ i + 1 ) = 0 (36) are simultaneously solved with the boundary condition, where one end of the myofibril is fixed and the other end is connected to a fixed linear spring (Figure 4B). The parameter γ T is the sarcomere friction in the transverse direction, and φ T is the deformation energy for stretching in the transverse direction (Washio et al., 2019). In our myofibril model, the rapid sarcomere lengthening in SPOCs is reproduced by the avalanche of reverse strokes at the peak contractile phase in which the population of MHs with large mechanical loads is high. From the SL-LS relationship (Eq. 32) and LS alignment with the adjacent sarcomeres (Eq. 33), the rapid lengthening of one sarcomere transversally compresses the neighboring half-sarcomere in the peak contraction phase. Then, from the SL-LS relationship, the transverse compression enlarges the longitudinal stretch, resulting in the load increments of the attached MHs. As a result, an avalanche of reverse strokes is triggered. In this way, the rapid lengthening wave travels in the myofibril model.

In the finite element ventricle model, the half-sarcomere model is imbedded in each tetrahedral element along the fiber orientation f represented by a unit vector in the reference configuration (Figure 4C). Thus, the half-sarcomere models pull each other at the element interfaces in the fiber orientation. As depicted in Figure 4C, the angle of the fibers relative to the equatorial plane varies depending on the depth of the wall, so that the direction of force developed in the wall covers a wide range. In the finite element analysis, the current position of each material point   X ∈ Ω in the reference (unloaded) configuration is represented by x = x ( X ) . Therefore, the stretch in the fiber orientation is given by λ = ∥ Ff ∥ , (37) where  F = ∂ x / ∂ X is the deformation gradient tensor. The active stress S a c t (as the second Piola-Kirchhoff stress tensor) is formulated based on considering the virtual work done by the active tension T a c t in the fiber orientation per unit volume in the reference space as follows: δ w a c t = − T a c t δ λ = − T a c t λ 1 2 δ ( f ⊗ f : F T F ) = − T a c t λ f ⊗ f : δE , (38) where δE = 1 2 ( δ F T F + F T δF ) is the infinitesimal increment of the Green–Lagrange strain tensor: E = 1 2 ( F T F − I ) . (39)

As a result, with the active stress represented by S a c t = T a c t λ f ⊗ f , (40) the virtual work by the active stress per unit volume is given by δ w a c t = − δE : S act . (41)

The equation of motion for the displacement, u ( X ) = x ( X ) − X , is given as: ∫ Ω δu ⋅ ρ u ¨   d Ω + ∫ Ω δE : S  d Ω − P L ∫ Γ L δu ⋅ n  d Γ L − P R ∫ Γ R δu ⋅ n  d Γ R = 0. (42)

Here, P L = L V P and P R = R V P are the blood pressures 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 ( 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: S ≡ S act + S pas + S vis , (43) where S act is given by Eq. 40, and S pas = ∂ W p a s ∂ E 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 systems of the body and the lung. See details of the deformation potential W p a s , the viscous stress S vis , and the coupling with the circulatory system in our previous work (Yoneda et al., 2021).

Assuming a periodical solution over a heartbeat, the time integration on a cycle with the substitutions of u ˙ and E ˙ , respectively, into δu and δE in Eq. 42 gives the relationship between the cardiac output C o u t and the energy production inside the cardiac muscle as follows: C o u t ≡ − ∮ P L V ˙ L d T − ∮ P R V ˙ R d T = − ∮ ∫ Ω E ˙ : ( S act + S vis )   d Ω d T . (44)

Here, the inertia and the passive energy terms disappear for the sake of the periodicity assumption. If the viscous energy loss is negligible, the following work done by the active stress is almost equal to the cardiac output. W a c t ≡ − ∮ ∫ Ω λ ˙   T a c t d Ω d T = − ∮ ∫ Ω E ˙ : S act   d Ω d T ∼ C o u t (45)

The above relationship between the blood dynamics and muscle dynamics of the left ventricle (the left ventricular free wall and the septum) in our ventricle model is depicted in Figures 7A,B, where the two cases of the myosin rod stiffness with k x b = 2.8 pN/nm (black lines) and k x b = 1.4 pN/nm (red lines) are compared. The solid lines and broken lines represent, respectively, the blood dynamics and the muscle dynamics. Roughly, in the mid-systole ([0.15 s, 0.2 s]), the average active tension reaches 30–40 kPa, with a systolic blood pressure of 12–15 kPa and an average half-sarcomere shortening velocity of 1 μm/s (see also Supplementary Video 2). The relationship between the active tension and the systolic pressure agrees with the simple estimation given by the Laplace law assuming the dimension of the left ventricle in the early-systole (40 mm-diameter, 10 mm-wall thickness, the fiber helix angle twisting from −60° (epicardium) to +60° (endocardium). For example, to support blood pressure of 15 kPa, the active tension T a c t must be roughly equal to 15  kPa  × ( 20   mm / 10   mm ) × 2 π 3 3   ∼ 36   kPa . Here, the factor 2 π 3 3 comes from considering the fiber helix angles. The complex distributions of the active tension and the half-sarcomere shortening velocity over the ventricular wall (see Supplementary Video 2) make it difficult to find the correspondence to the tension-shortening velocity relationship under isotonic conditions (Figure 6B). Although the cardiac outputs and the muscle works of the two cases (the solid lines in Figures 7A,B, and the broken lines in Figure 7C) are not so different, the difference in energy consumption (the solid lines in Figure 7C) is remarkable, as indicated also in Figure 6D. The comparison of the time transients of the population ratio of the binding states (Figure 7D) indicates a shorter average dwell time in the XB_PostR1 state and a larger population in the XB_PostR2 state for the case of k x b = 1.4 pN/nm than those for the case of k x b = 2.8 pN/nm, resulting in the remarkable energy loss. Therefore, in this study, we analyzed the effects of the trapping mechanism when k x b = 2.8 pN/nm.

To assess the effectiveness of the trapping mechanism, a stretch-activation test was performed for a single half-sarcomere model consisting of 32,768 filament pairs. Here, stretch lengths of 5, 6, 7, or 8 nm were applied over a 1-ms time interval after the active tension had sufficiently matured. During the simulations, the Ca²⁺ concentration ([Ca²⁺]) was held at a constant value of 0.3   μM (Figure 8A) or 0.4   μM (Figure 8B). A greater increase in tension (the maximum in Phase 3) was observed as the stretching increased up to 7-nm stretch, corresponding to roughly +13% for a 6-nm stretch and +15% for a 7-nm stretch under [Ca²⁺] = 0.3 μm. Here, 7 nm is 0.7% of the half-sarcomere length. These increases lasted for a few seconds after the rapid stretching. In agreement with the experimental results given by Stelzer et al. (2006), our numerical results reproduced “phase 2 (Figure 1A)” in which the force decreased to the steady-state level before stretching. The increase in tension observed by Stelzer et al. (2006) was roughly 16% for a 2% stretch and 8% for a 1% stretch. The tension increases stopped at an 8-nm stretch in the simulation result. In the numerical model, the stretch is directly linked to the shift of the thin filament relative to the thick filament (Figure 3A), while the stretch of the sarcomere may be relaxed by intermediate substances such as Z-lines or intercalated disks in the experimental conditions, as analyzed in the Discussion Section 4.3. Note that we assumed a cross-bridge model under physiological conditions identical to those of a living heart, which differ from the experimental conditions in many aspects, such as temperature.

Focusing on the result for a 7-nm stretch (Figures 9A–F), compared with the pre-stretch steady state, the lasting increase in tension apparently arose from a lasting increase in the population of the XB_Trap state (Figures 9C,D). The effect of the XB_Trap state is more emphasized when the individual contributions in the active tension of the attached states are plotted (Figures 9E,F). In the half-sarcomere model, we assumed that one thin filament is surrounded by 76 (=   38   × 2 ) myosin molecules. In the steady state before stretching, filament forces of roughly 20 and 22 pN are generated by myosin molecules in the XB_PostR1 (2.9%) and XB_PostR2 (1.7%) states, respectively, when [Ca²⁺] = 0.3 μm. Thus, the average force per molecule is 9.1 and 17.0 pN, respectively, in the XB_PostR1 and XB_PostR2 states. Because we assumed a force constant of k x b = 2.8 pN/nm for the myosin rod distortion x , the above forces are generated by distortions of 3.3 and 6.1 nm, respectively. After the stretch, a filament force of 12 pN is produced by the myosin molecules in the XB_Trap state (∼0.9%). Thus, the average force is roughly 17.5 pN per trapped myosin molecule. This average force is slightly larger than experimentally observed maximal forces (∼15 pN) (Hwang et al., 2021). As analyzed in Section 2.4 and Section 2.5, our setting of the force constant ( k x b = 2.8 pN/nm) seems reasonable with respect to the efficiency in both the isotonic contraction and heart pumping. Furthermore, it may be difficult to measure the maximal force produced by a single myosin molecule under the sarcomeric condition. In particular, the force measured in the filament direction depends on the angle between the myosin rod (S2) and the thick filament, and the maximal force under a similar condition to the stretch activation environment has not been reported.

As observed by Stelzer et al. (2006), the MC cross-bridge model also reproduced the different behaviors in Phase 3 (Figure 1A) in the different activation levels. In our model, this difference is made by the differences in the time courses of the population ratio of the attached states (Figures 9C,D) that are produced by the cooperativity effects, as shown in Figure 5A.

The effects of the trap mechanism on SPOCs of a single myofibril model consisting of 40 half-sarcomeres were investigated (Figures 10A,B). Here, 2048 filament pairs were imbedded in each half-sarcomere model. SPOCs were produced for all Ca²⁺ concentrations ([Ca²⁺]) in the no-trap model. In contrast, SPOCs were reproduced only for the intermediate concentration ([Ca²⁺] = 0.4 μM ) between the states of relaxation and contraction in the trap model, as observed by Fabiato and Fabiato (Fabiato and Fabiato, 1978). Therefore, the trap mechanism may contribute to the [Ca²⁺] dependence of the SPOCs by inhibiting the reverse stroke for certain amounts of myosin molecules in the XB_PostR1 state. For the low Ca²⁺ concentration ([Ca²⁺] = 0.3 μM ), the lengthening was prevented by myosin molecules trapped in the XB_Trap state despite the oscillating XB_PostR2 concentration in the trap model, as indicated by a small increase in the XB_Trap state for each reduction in the XB_PostR2 state (Figure 11A). For the intermediate Ca²⁺ concentration ([Ca²⁺] = 0.4 μM ), the SPOCs are similar for both cases (Figure 11B). For the high Ca²⁺ concentration ([Ca²⁺] = 0.6  μM ), the sarcomeres were shortened to the minimum length, which is almost equal to the thick filament length of LM ∼ 1.6  μM , as shown in Figure 1B. At this length, the population of the attached MHs is small because of the small single overlap region of the two filaments. Thus, lengthening was prevented in the same way as that observed for the low Ca²⁺-concentration case in the trap model (Figure 11C).

Beating-ventricle simulations were performed using a finite element ventricle model with the same setup as our previous work (Yoneda et al., 2021). In each element, a sarcomere model consisting of 16 filament pairs was imbedded along the appropriate fiber orientation f . The distribution of the fiber orientations was found by an optimization algorithm based on the isovolumetric active tension (Washio et al., 2020) according to the impulses given by the active tension. Portions of the helical fiber structure are depicted in Figure 4C. The heart rate was set to 60 beats per minute, and the Ca²⁺ transient (Figure 4C) generated by the mid-myocardial cell model proposed by ten Tusscher panflov, (2006) was applied. Transmural delays of the Ca²⁺ transient determined by the distances from the endocardial surfaces of the left and right ventricles under a transmural conduction velocity of 52 cm/s, as measured by Taggart et al. (2000), were adopted. By comparing the trap and no-trap models in Figures 12A,B, one can see that the trap mechanism contributes to maintaining both the high pressure in the last half of the systolic phase and the rapid pressure decrease at the end of the systolic phase (Figure 12A). As a result, the blood volume ejected from the left ventricle in the trap model increased from 71 to 78 mL , while the ATP energy consumption of the left ventricular wall was almost equivalent (Figure 12C). This trend implies that the trap mechanism increases blood ejection without increasing energy consumption. It should be noted that the ATP consumption rates were computed by counting the detachments of MHs in the XB_PostR2 state transferred to those in the N_XB state, which was controlled by the rate constant r N X B and the forced detachments defined by Eqs 14 and 15. As shown in Figures 12A,B slight increase in the population of MHs in the XB_Trap state can be seen at the end of the systolic phase (Figures 12D,E). This increase corresponds to a prolonged ejection time in the trap model, as shown in Figure 12A.

In the systolic phase, the cardiac myocytes within the shared fiber bundle support each other by pulling their neighbors via contractile forces, with the active tension in Eq. 22 playing the greatest role. Therefore, based on the mechanical equilibrium condition in the fiber direction, the active tensions must be almost equal. Consequently, if there is a loss of active tension at one point of the fiber bundle prior to the remaining parts in the early-systole or in the end-systolic phase, this portion would be lengthened quickly, and the sarcomeres in the remaining parts would be shortened until a mechanical equilibrium is reached. Because this transition accompanies decreases in the active tension owing to the loss of distortion in the myosin rods, stopping this process as early as possible is desirable to maintain blood pressure. The trapping mechanism can achieve this goal, as shown in Figures 13A–C, which depict the distributions of the active tension generated by the individual states XB_PostR2 and XB_Trap in the end-systolic phase. The degenerated forces of the XB_PostR2 state in the no-trap model (Figure 13A) are compensated by the force generated by the XB_Trap state in the trap model (Figure 13B). Furthermore, the forces of the XB_PostR2 state in the trap model around the degenerated regions in the no-trap model (Figure 13C) are well maintained, owing to the contributions of the reinforced regions via the trap mechanism. As shown in the frequencies of transitions (Figure 14), the reverse power stroke contributes to diastolic relaxation to the same extent as the detachment from the XB_PostR2 state (Figure 14A). Although the difference between the trapping and escaping frequencies averaged over the left ventricular wall (Figure 14B) is relatively small, the two peaks of the difference in Figure 14C agree with the differences in LVP between the trap model and the no-trap model in Figure 12A. Figure 14C also indicates that a certain extent of the trapped MHs is forced to detach because of the extreme distortion (Eq. 15).

In the stretch activation tests, the step length dependence was reproduced (Figure 8), similar to the experimental results of Stelzer (Stelzer et al., 2006). However, the minimum values for the rapid force decay (phase 2) in the numerical results were lower than the original isometric force before stretching and higher than the isometric force in the experimental results of Stelzer. Conversely, minimum values smaller than the original isometric force were reported for wild-type myocardium by Mamidi et al. (Mamidi et al., 2018). Thus, the level of the minimum force in phase 2 compared with the original isometric force seems to depend on the experimental conditions. In our model, the decay was determined from decreases in the XB_PostR1 and XB_PostR2 states owing to the reverse strokes from XB_PostR2 to XB_PostR1 and XB_PostR1 to XB_PreR. Meanwhile, the loss of force was compensated by the increase in population and the distortions of the trapped myosin molecules in the XB_Trap state (Figures 9B,C). Here, the rates of reverse strokes are limited by the upper bounds r ¯ b , i , as in Eq. 19, and the population in the XB_Trap state is determined by the parameter values associated with the trap mechanism ( s t r a p , h trap , and h escape ) in Eqs 3 and 4 and the forced detachment ( x X B , t , c X B , t , and a X B , t ) in Eq. 15. In our trap model, these parameters were carefully chosen, focusing not only on the stretch activation phenomena, but also on the SPOCs and performance in the beating ventricle model.

In our model, the trap mechanism was added to the XB_PostR1 state, which is the state after the first power stroke. Alternatively, it might be conceivable to add a similar trap mechanism in the XB_PostR2 state (the state after the second power stroke) if a force-dependent detachment rate constant (Greenberg et al., 2014) is adopted in the numerical model. Hwang et al. (2021) reported a much gentler increase in the reverse stroke rate constant of the first stroke compared with that of the second stroke, with a force increase of 8–14 pN. Their single-molecule experimental results support the adequacy of adding the trap mechanism to the XB_PostR1 state.

The numerical SPOC results indicate that the trap mechanism affects the calcium activation sensitivity for the relaxation dynamics in an advantageous manner. The trap mechanism prevents sarcomere lengthening at high or low levels of calcium activation (Figure 10). The sarcomere lengthening maximizes the single overlap region between the two filaments, resulting in a facilitation of re-activation. Thus, the prevention of lengthening at low calcium levels (Figure 11A) stops weak re-activation, leading to a smooth relaxation of muscle in the diastolic phase (Figure 12A). Conversely, the prevention of lengthening at high calcium levels (Figure 11C) stops sarcomere shortening in the neighboring cardiomyocytes, causing a retention of high active tension in the fiber bundle (Figures 13B,C).

The above-mentioned features for heartbeats are similar to those found in our previous work (Washio et al., 2018) with the Langevin dynamics model. However, the new MC model produced better results in the tension response for the step length changes and SPOCs owing to the simpler adjustment of related parameters. In a previous study, MC simulations were performed for a three-dimensional half-sarcomere model consisting of three myosin filaments and 13 thin filaments implemented with 360 myosin molecules and 1,170 binding sites to examine the impact of filament compliance on Ca²⁺ activation (Chase et al., 2004). However, these simulations were conducted only under steady-state conditions and were not coupled with a macroscopic finite element model (FEM). Moreover, in the beating ventricle simulation, the computation cost per beat was 1.5 h when 320 cores of a conventional parallel computer system were used for the MC model (Yoneda et al., 2021), while the cost was 105 h when 1920 cores of the same computer system were used (Washio et al., 2018). The former simulation was performed for a bi-ventricular FEM consisting of 45,000 tetrahedral elements with 16 filament pairs, while the latter was performed for a smaller FEM consisting of 7,900 elements with eight filament pairs. The improvements of the MC model, such as those demonstrated in the current study, will lead to better predictions in clinical applications (Kariya et al., 2020).

Smith et al. argued that multiple working strokes are required for a cross-bridge cycle to satisfy energetic constraints and demonstrated that a cross-bridge cycling model with three strokes can reproduce various experimental findings observed for frog muscle, including absolute values of active tension, stiffness, and ATPase rate; the phase-2 tension response to a length release; and the transient tension rise during ramp stretching (Smith et al., 2008a; Smith et al., 2008b). This line of reasoning led us to adopt a model with two strokes with lengths of 6.0 and 4.0 nm to realize the 10-nm stroke suggested by the crystal structure of myosin molecules. In the three-stroke model of Smith et al., the first two strokes, which each have a length of 5.0 nm, occur around the Pi release, and a third small stroke was implemented to account for the strain-dependent ADP release rate. Although we did not explicitly define the relationship between strokes and the nucleotide released in our model, these points should be examined by comparing our results against various experimental findings in future work.

The average force of the trapped MH was larger than the maximal force ever observed experimentally. In our model, this discrepancy comes from our stiffness parameter k x b = 2.8 pN/nm for nearly a 10 nm stretch, assuming linear elasticity. The force function, F r o d , in Eq. 20 may be improved by reducing the stiffness for substantial distortions. Further investigation and modification of the detachment rate function in Eq. 15 represent future tasks.

In our model, we didn’t account for the elasticity of the components, such as the thin and thick filaments and the Z-band in the sarcomere. Namely, we assumed that all of the sarcomere components (except for the myosin rods) are rigid. Thus, the macroscopic stretch change Δ λ is directly reflected in the distortion increase Δ x = S L 0 Δ λ / 2 . However, in the actual setting, the distortion increase might be somewhat relaxed because of the elasticities. For example, if we assume a linear elasticity of the Z-line with a spring constant of k Z and N x b attached myosin molecules per thin filament, we have: { N x b k x b Δ x = k Z Δ Z , Δ x + Δ Z = S L 0 Δ λ 2 . (46)

By eliminating the Z-line distortion increment Δ Z in Eq. 46, we obtain: ( 1 + N x b k x b k Z ) Δ x = S L 0 Δ λ 2 . (47)

Thus, the distortion increment Δ x decreases as the number of attached myosin molecules increases. Here, the magnitude of the distortion increment Δ x corresponds to the strength of trapping for myosin molecules in the XB_Trap state. Because N x b is smaller for a lower activation, Eq. 47 may give a straightforward explanation of the experimental findings reported by Stelzer (Stelzer et al., 2006), who reported that the stretch activation is most pronounced at low levels of Ca²⁺ activation. Apart from these considerations, the effect of filament compliance on the realignment of cross-bridges reported by Chase et al. (Chase et al., 2004) is an important issue that should be included in future modeling studies.

Razumova et al. modeled and compared three possible mechanisms of cooperativity: 1) interactions between adjacent T/T RUs, 2) interactions between adjacent cycling cross-bridges, and 3) a facilitation of the transition to the on-state of a RU by the adjacent attached cross-bridge, which all control the open and closed states of the thin filament (Razumova et al., 2000). Their results clearly demonstrated distinct roles of these interactions in the maximal force and cooperativity; however, their formulations are conceptual and thus do not represent specific molecular interactions. McKillop and Geeves proposed a three-state (blocked, closed, and open) model of thin filament activation, which is compatible with X-ray diffraction data (McKillop and Geeves, 1993). Smith et al. further attempted to establish the relation between RU–RU interactions and physical entities (Smith et al., 2003) by modeling the flexible-chain-like structure of the tropomyosin molecule with a continuous-flexible-chain model. In this regard, the RU–RU model used in this study is also empirical and lacks a relation to a physical entity. Moreover, only part of the above-mentioned mechanism of cooperativity was considered.

In the finite element ventricular model, a single half-sarcomere model was imbedded in each tetrahedral element, and the half-sarcomere length changed according to the stretch of the element in the fiber direction. This is nothing but assuming that the movements of sarcomeres contained in each element are perfectly synchronized. In reality, there may be time lags in the length changes between the neighboring sarcomeres particularly in the relaxation phase. In our future work, the issue will be studied by using the homogenization method (Washio et al., 2013) whereby the bundle of myofibril models is imbedded in each element.