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

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

(S1.5) The transitions between the N XB and P XB states are affected by the status of the T/T unit above it through modifications of np and pn , as well as by the state of the neighboring MHs through the integer , where (= 0, 1, or 2) is the number of neighboring MHs in a weakly or strongly bound state.
(S1.8) pn = � pn1 if the T/T unit above is in the Ca-on state, pn0 otherwise. (S1.9) Here, = 1 if the MH is located at the single overlap region with the thin filament, otherwise = 0. This, along with or − ( = 40), represents the nearest-neighbor cooperativity of the MHs, following Rice (2003), which plays an important role for the force-pCa relationship, as shown in S2. We assume that one thin filament in the three-dimensional arrangement corresponds to two thin filaments in our half-sarcomere model. This is because we are assuming that cooperative behavior exists along the tropomyosin and tropomyosin molecules wrapped around the thin filament in a double spiral fashion, and one of the spirals is considered in our half-sarcomere model. The constants np0 , np1 , pn0 , and pn1 are determined from , basic , and , as follows: np0 = basic , np1 = basic , pn0 = pn1 = basic 2 (S1.10) Here, > 1 controls the degree of cross-bridge inhibition for T/T units in states other than Ca-on, and controls the ratio of binding states of the MHs. The greater the value of , the larger the ratio of binding states for a given Ca 2+ concentration. To reproduce the sarcomere length ( ) dependence in the contraction force, is given as a function of by (S1.11)

S1.2 Nonlinear rod strain energy model
The rod strain energy was nonlinear with the generated force, following Kaya (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 region along the myosin rod (Kaya and Higuchi, 2010), as depicted in Figure S1.2. The strain energy is given by integrating the force from = 0 given by (S1.12) where and 1 are determined from the other parameters so that the function and its first derivative are continuous at = 0 and = − 1 : ) . (S1.13) The parameters adopted in our model are = 2.8 pN/nm, 1 =4.35 nm, and = 0.05 . The profiles of the resultant force and potential functions are depicted in Figure S1.2. Figure S1.2 Force and the potential energy for the myosin rod strain .

S2. Verification of the Basic Properties of the Half-Sarcomere Model
The basic properties of the half-sarcomere model were verified as follows. For the Langevin dynamics, the trap model was applied. All the simulations were executed with the MTS scheme (∆ = 0.25 ns, ∆ = 5000 ns:except for S2.4 where ∆ = 25ns was used). The contractile force per one thin filament, T,∆T , in Equation (25) of the main text was averaged over 48 thin filaments ( = 48). Note that the contractile force computed in the numerical experiments corresponds to the force for one of the double spirals along the thin filament. Therefore, the factor 2/ 0 ( 0 = 0.001µm 2 , the cross-sectional area of a single thin filament) should be multiplied if one evaluates the tension within the sarcomere. All the numerical simulations were performed with the same parameters as those in the main text, which attempted to reproduce the muscle contraction at normal body temperature, while most of existing experimental data were obtained at lower temperatures. Therefore, we avoid quantitative comparisons with the existing experimental results.

S2.1 SL and Ca 2+ Force Relationships
To obtain the force in each case, the simulations were started in the relaxed state (all MHs were in N XB ) at = 0 s with a given Ca 2+ concentration [Ca], and the contractile force was averaged over the time interval [0.75 s, 1 s]. In Figure S2.

S2.2 Isometric Twitch
A series of isometric twitches were simulated for the Ca 2+ transient generated by the mid-myocardial cell model proposed by ten Tusscher et al. (2006) while varying the SL ( Figure S2.2). As observed by Janssen et al. (1995), the 50% relaxation level was prolonged with increased SL and force, whereas the time to reach the 90% rise level was only slightly prolonged.

S2.3 Force-Velocity Relation
To predict the dependence of the force on the shortening velocity in the trap model, the halfsarcomere excited with a constant [Ca] was quickly released from the initial SL value of 2.3µm with various constant loads ( Figure S2.3(A)). Under this isotonic condition, the shortening velocities were measured ( Figure S2.3(B)). Physiologically reasonable curves similar to the experimental results given by Piazzesi et al. (1995) were reproduced, with the maximal half-sarcomere shortening velocity close to 3 µm/s. Figure S2.4(A) is similar to the experimental results obtained by Reconditi et al. (2004). As can be seen on the left side of Figure S2.4(B), the initial plateau region 0 − 5 ms of the SL transient corresponds to the duration shifting towards the new equilibrium of the isotonic contraction after the non-equilibrium state was induced by the quick release, where a significant number of MHs in the pre-power stroke state (Pre) and the post-first power stroke state (PS 1 ) were transferred rapidly to the post-second power stroke state (PS 2 ). These rapid transitions were induced by a rapid decrease in the rod strain in the Pre and PS 1 states, as shown in Figure  S2.4(C). In the case of a smaller load (50 pN/Filament), the averaged in PS 2 decreases until the increase in the shortening velocity stops. Although the averaged in PS 2 is negative after 10 ms, the half-sarcomere produces a positive contractile force equal to 50 pN/Filament. This is due to the nonlinearity of the spring force ( Figure S1.2).

S2.4 Tension Recovery for the Step Length Change
To examine the tension recovery after quick shortening, the half-sarcomere model was quickly shortened by 3, 5, or 7 nm within 0.1 ms after the contractile force was developped under a constant [Ca] = 10 µM ( Figure S2.5(A)). The shortening effect is clearly recognizable as the sudden decrease of averaged rod strain in the pre-power stroke state (Pre) ( Figure S2.5(B)). This induced the quick recovery of tension given by the collective power stroke transition from Pre to PS 2 (Figures S2.5 C and D ).

S3. Derivation of the Active Stress Tensor and its Stiffness
The muscle power P per unit volume delivered by the contractile tension (per unit area in the reference configuration) in the direction specified by the unit vector (in the reference configuration) is given by the product of and the time derivative of the stretching = � � along : Here, the : symbol denotes the dot product of the two tensors. Thus, the active stress can be represented by the first Piola-Kirchhoff stress tensor: As a result, Equation (S3.1) can be rewritten as The application of the first Piola-Kirchhoff stress tensor act to an area element in the reference configuration, in which the normal vector points outward from the area, yields the traction force in the current configuration, as follows: Here, � = 1 is the unit vector directed along the fiber orientation direction in the current configuration, and ( • ) is proportional to the number of actin filaments that pass through the unit area perpendicular to . Thus, Equation (S3.4) indicates that the active stress tensor in Equation (S3.2) is based on the muscle power , which corresponds with the usual stress tensor defined by the traction force.
From Equation (S3.2), the active stress tensor is represented by the second Piola-Kirchhoff stress tensor: To perform an implicit time integration scheme, such as the Newmark-β method, we need to compute the derivative of the stress tensor with respect to the nodal displacement vector , as well as its time derivative ̇. In Equation (S3.5), is a function of , and is a function of and ̇ , because in Equation (38) contains ̇ (which is a function of and ̇) : Thus, the derivative of act can be computed by Now, from Equation (37), the derivative of is given by and the derivatives of and ̇ are, respectively,

S4 Passive and Viscous Parts of the Ventricle Model
The passive stress tensor is given by the deformation potential function as The potential function is determined by a macroscopic passive potential: where ̃1 is is the reduced invariant defined as with the right Cauchy-Green deformation tensor = .
is a quadratic form of the Green-Lagrange strain tensor (Usyk et al., 2000): where the components are defined based on the fiber-sheet structure of the muscle walls as Here, { , , } is the orthonormal basis of the ventricle walls that determines the fiber and laminar structures. denotes the volume ratio of the sarcomere in the ventricle muscle wall, and sar is the potential function introduced by the material properties of the sarcomere mechanics. The thick filaments are connected to the Z-line with another filamentous protein called titin. These proteins are believed to prevent overstretching of the sarcomere, and the presence of the thick filaments is assumed to prevent too much shortening. These effects are modeled by introducing the potential as a function of the stretch : Here, mf is the shortening threshold of of the thick filament, and titin is the elongation threshold of the titin.
For the viscous part, the Newtonian viscosity is given by where is the viscosity coefficient, and is the deformation velocity tensor defined as: Note that the derivatives are given with respect to the Eulerian coordinates .

S5. Circulatory System
As the atrial model, the formulations by Kaye et al. (2014) were applied. The left and right atrial pressure is related to the time-varying elasticity and the chamber volume as follows. Hereafter, the subscript "A" stands for "LA" or "RA" for the left or the right atrium, respectively.
Here, the functions , and , give, respectively, the pressure at the end of the diastolic and systolic phases for chamber volume . These functions are defined by The time-varying elastance is given as a function of the time as follows: Here, 0 is the start time of atrial excitation, max is the time to maximal chamber elastance, and is the time constant of relaxation. Together with the other part of the circuit model of Figure 5 in the main text, the overall equations for pulmonary circulation are given by Here, denotes the pulmonary venous compliance, and denotes the increase of the blood volume from zero pressure. � and � , the flow rates in the case of no rectification, are given by Similarly, the overall equations for the systemic circulation are given by with flow rates (for no rectification) The parameters adopted in our simulation are listed in Table S5.1. The parameter values were chosen to reproduce the temporal changes in ventricular blood pressure for a standard healthy heart.

S6. Newmark-beta Time Integration for the Macroscopic Nonlinear Equation
In our process, the combined system composed of the mechanical equations for the biventricular model and the blood circulation equations were simultaneously solved implicitly using the Newmarkbeta scheme. The combined system of equations is given by the following three formulas, for the biventricular FEM, pulmonary circulation, and systemic circulation models. Thus, the missing components in have no influence on the computational process, although they are also computed together with the other unknowns.
In the Newmark-beta time integration scheme, the variables ,̇, and ̈ are updated from those at time to time + ∆ so that both the interpolation relations in time (Equations (S6.8A-B)) and the equilibrium condition (Equation are fulfilled. Here, and are the interpolation weights which determine the numerical accuracy and stability of this scheme. In our case, we adopted = 0.3025 and = 0.6. To solve this system, Newton iterations are applied from the initial guess: From Equations (S6.8A-B), the updates of ̇ and at the first iteration are given from that of ̈ as follows.
The linearization of Equation (S6.8C) at the initial guess (Equation (S6.9)) with the increments in Equation (S6.10) gives the linear equation for ̈( 1) : where the matrices are given by the derivatives of , as follows ( = 0).
For the first iteration, we may solve Equation (S6.11) to obtain ∆̈( 1) , then we update the variables according to: Although the initial guess in Equation ( In our simulation, the iterations were terminated when the residual (the right-hand side in Equation (S6.15)) was sufficiently small. Figure S7.1 compares the elapsed times for simulating 1 ms at the peak of the systolic phase for various numbers of cores, ranging from 240 to 3840. Here, 32 elements were assigned to each core in the minimal configuration (#cores=240), while 2 elements were assigned to each core in the maximal configuration (#cores=3840). During the computation required for 1 ms, calculation of the macroscale portion occurred 250 times, and required about 65 s, regardless of the numbers of cores used. Because the elapsed time for the macroscale part was sufficiently small, even in the maximal configuration, a further reduction can be expected if a much larger number of cores are available. Figure S7.1(B) shows the effectiveness of using a large macroscale time step size, ∆ . If we choose a smaller time step, for example, ∆ = 312.5 ns ( = ∆ ∆ ⁄ = 1250), the macroscale computational portion occupies more than half of the time. Thus, an MTS scheme that admits a large macroscale time step size is necessary for this application.