Biophysically detailed mathematical models of multiscale cardiac active mechanics

We propose four novel mathematical models, describing the microscopic mechanisms of force generation in the cardiac muscle tissue, which are suitable for multiscale numerical simulations of cardiac electromechanics. Such models are based on a biophysically accurate representation of the regulatory and contractile proteins in the sarcomeres. Our models, unlike most of the sarcomere dynamics models that are available in the literature and that feature a comparable richness of detail, do not require the time-consuming Monte Carlo method for their numerical approximation. Conversely, the models that we propose only require the solution of a system of PDEs and/or ODEs (the most reduced of the four only involving 20 ODEs), thus entailing a significant computational efficiency. By focusing on the two models that feature the best trade-off between detail of description and identifiability of parameters, we propose a pipeline to calibrate such parameters starting from experimental measurements available in literature. Thanks to this pipeline, we calibrate these models for room-temperature rat and for body-temperature human cells. We show, by means of numerical simulations, that the proposed models correctly predict the main features of force generation, including the steady-state force-calcium and force-length relationships, the length-dependent prolongation of twitches and increase of peak force, the force-velocity relationship. Moreover, they correctly reproduce the Frank-Starling effect, when employed in multiscale 3D numerical simulation of cardiac electromechanics.


S2 Appendix. Parameters calibration
To calibrate the parameters of the proposed models, we will restrict ourselves to experimental measurements from the literature coming from intact cardiac cells since the skinning procedure alters in a significant (and only partially understood) manner the activation and force generation dynamics [1][2][3]. Moreover, thanks to the technique of flura-2 fluorescence, it is nowadays possible to measure the intracellular calcium concentration without depriving the cell of its membrane, and it is also possible to inhibit the sarcoplasmic reticulum calcium uptake by cyclopiazonic acid, so that the calcium level can be controlled without the need of skinning the cells [3,4].

Calibration of the XBs rates
In [5] we have shown that the parameters of the distribution-moments equations describing the XB dynamics can be calibrated starting from five quantities, that are described in the following. Under isometric conditions, we consider the isometric tension T iso a := a XB µ 1 and the fraction of attached XBs µ 0 iso := µ 0 , where µ 0 and µ 1 are the steady-state solution for v = 0. The force-velocity relationship is characterized by the maximum shortening velocity v max , and by v 0 , its intersection with the axis T a = 0 of the tangent of the curve in isometric conditions, defined as: where T a (v) denotes the steady-state tension for velocity v. Finally, as discussed in [5], the response to fast inputs is characterized, in the small-velocity regime, by the normalized slope of the tension-elongation curve after a fast step. Such quantity is defined as follows: by applying a step in length ∆L to an isometrically contracted muscle in a small time interval ∆t, we define by T a (∆t) the tension recorded after the step is applied, and we define the normalized stiffness as: .
As discussed in [5], when the small-velocity regimes is considered, the quantityk 2 corresponds to the slope of the T 2 -L 2 relationship.

1/6
Parameter Value Units Reference T iso a 120 kPa [4] Table A. List of the experimental data used for model calibration.
In conclusion, as shown in [5], by acting on the five parameters µ 0 f P , µ 1 f P , r 0 , α and a XB we can fit experimentally measured values concerning the isometric solution (T iso a and µ 0 iso ), the force-velocity relationship (v max and v 0 ) and the fast transients response (k 2 ). All the above mentioned experimental setups are such that the thin filament activation machinery can be considered in steady-state. Indeed, [Ca 2+ ] i is constant in all the cases and, concerning SL: it is also constant under isometric conditions; constant shortening experiments are typically performed in the plateau region of the force-length relationship, and thus the effect of changes in SL is irrelevant; fast transient experiments are carried out at a time-scale such that the activation variables can be considered constant, since they are characterized by a much slower dynamics [6,7]. Therefore, in these cases, the values of P i , k PN T,i and k N P T,i can be considered as fixed in Eq. (11) (and similarly in Eq. (21)).
We notice that, while for the models considered in [5] the relationship between the five parameters and the five experimentally measured values can be analytically inverted, in this case we find the values of the parameters with a numerical strategy. Specifically, to find the steady-state solution we solve Eq. (11) by setting to zero the time derivative terms; we consider the exact solution after the fast transient (the linear ODE system (11) can be solved analytically); we approximate the derivative appearing in the definition of v 0 andk 2 by finite differences. Finally, we solve, for the five parameters, the nonlinear system of equations linking the five measured values with the parameters themselves. With this aim we employ the Newton-Raphson method, by approximating the Jacobian matrix by means of finite differences.
The experimental data used to calibrate the model are reported in Table S2-1, together with a reference to the source in literature. As we mentioned before, we employ data coming from room-temperature intact cardiac rat cell, apart from µ 0 iso (acquired from skeletal frog muscle), for which we did not find measurements from cardiac muscles. However, as shown in [5], this variable only affects the value of the microscopic variables (i.e. µ p i,α ), but not that of the predicted active tension T a . We notice that the constants v max , v 0 andk 2 are normalized with respect to T iso a and are thus valid for a wide range of activation levels (see Introduction). Conversely, the value of T iso a is associated with a SL in the plateau region and to saturating calcium concentration. Therefore, when we calibrate the parameters we set [Ca 2+ ] i = 10 µM and SL = 2.2 µm.

Calibration of the RUs rates (steady-state)
The steady-state characterization of the muscle tissue activation is represented by the force-calcium and force-length relationships (see Introduction), whose main features are the behavior for SL in the plateau region (characterized by the tension for saturating calcium T iso a , the calcium sensitivity EC 50 and the cooperativity coefficient n H ) and the effect of SL (on the saturating tension T iso a and on the calcium sensitivity EC 50 ). We recall that, thanks to our strategy, the tension for saturating calcium concentrations T iso a in the plateau region of SL is automatically fitted. The effect of k d 2/6 is that of shifting the force-calcium curves with respect to the log [Ca 2+ ] i axes since it only appears in combination with [Ca 2+ ] i in the model equations. Therefore, the value of k d can be easily calibrated to match the experimental data as it only affects EC 50 . The effect of γ, on the other hand, is that of tuning the amount of cooperativity (indeed, it acts on n H ). The role of the remaining parameters (Q and µ) is more involved and cannot be easily decoupled, as they affect the cooperativity, calcium sensitivity, and the asymmetry of the force-calcium relationship below and above EC 50 [10]. Moreover, in the SE-ODE case, they also act on the SL-driven regulation on calcium sensitivity.
In the following we set µ = 10, coherently with the fact that the experimentally measured dissociation rate of Tn from calcium varies of one order of magnitude in different combinations [11]. For the SE-ODE model, we set γ, Q and k d , to fit the steady-state force-calcium measurements of [4] (referred to the two different values of SL of 1.85 and 2.15 µm) from intact rat cardiac cells at room temperature.

Calibration of the RUs rates (kinetics)
We consider the isometric twitches of intact rat cardiac muscle fibers reported in [12] for different values of SL, and with [Ca 2+ ] o = 1.0 mM. Since the corresponding calcium transients are not reported, we consider the calcium transient reported in [13] for the same muscle preparation. As the reported trace is much affected by noise, we fit it with the following idealized transient [14]: with the constants values c max = 1.35 µM, t c 0 = 0.05 s, τ c 1 = 0.02 s, τ c 2 = 0.19 s. Then, we sample the candidate parameters space (k basic , k off ) ∈ [0, 80 s −1 ] × [0, 300 s −1 ] by a MC strategy, for each parameters value we compute the tension transients corresponding to SL = 1.90, 2.05 and 2.20 µm and we compare them with the experimental measurements from [12]. We consider the following discrepancy metrics, where T i,mod a (t) denotes the tension predicted by the model for the i-th value of SL and T i,exp a (t) denotes the experimentally measured one: The first metric coincides with the L 2 error over time, while the second one evaluates the error of the predicted force peak. We notice indeed that the parameters k basic and k off also determine the force peak attained during a sarcomere twitch: the most rapid the tissue activation is, the more the force-calcium curve stays close to the steady-state curve and thus it reaches higher force values before the relaxation begins. As criterion to select the best parameters values, we consider as overall metric a weighted combination between the two, given by E tot = (E 2 L 2 + w 2 peak E 2 peak ) 1/2 , where we set w peak = 5. The obtained values of the discrepancy metric E tot in the parameters space for both the SE-ODE and the MF-ODE models are reported in the below picture. We notice that the level curves do not clearly identify an optimal pair (k basic , k off ), rather these exhibit a wide region in the parameters space producing very similar results. Given the uncertainty in the measurements of both force and, mostly, calcium, it makes no sense to select the best parameters by blindly selecting the pair that realizes the smaller discrepancy from experimental results. Therefore, we supplement the results of Fig. S2-1 with direct measurements of calcium binding rates to Tn, showing that k on = k off /k d lies between 50 and 200 µM −1 s −1 [11]. On this basis, we restrict the region of candidate values and we select the parameters reported in Tab. 3.

MF-ODE SE-ODE
The predicted isometric twitches obtained with the selected values of the parameters are compared with the experimental data in Fig. 15. We notice here that the MF-ODE model accomplishes a better fit of experimental data than the SE-ODE model. This is an effect of the phenomenological tuning of k d of Eq. (29), that allows for a significant increase of calcium sensitivity and, consequently, of twitch duration, for higher values of SL. Nonetheless, also the SE-ODE model predicts, even if to a lower extent, the experimentally observed prolongation of twitches at higher SL, without any phenomenological tuning of the calcium sensitivity.
We notice that there is room for improvement in the calibration of the kinetic parameters k basic and k off , which could be better constrained in presence of more abundant and more reliable experimental data and when a deeper understanding on the determinants of the kinetics of activation and relaxation will be available. Nevertheless, the calibration of k basic and k off for the rat model does not affect the quality of the human model, since those two parameters are the only ones to be completely re-calibrated for the human model.

Human model at body temperature
In order to adapt the parameters calibrated from rat data at room temperature to a body-temperature human model, we first focus on the steady state, to reflect a higher calcium sensitivity. For this purpose, we employ the data reported in [15], which, however, refer to skinned cells. In order to estimate the effect of skinning on k d , we compare the calcium sensitivity measured for room-temperature rat cardiac cells in skinned [16] and intact preparations [4] at SL = 2.15 µm and we assume that the same relationship holds for skinned versus intact, body-temperature human cells. Finally, we rescale the values of k d to reflect the estimated change in calcium sensitivity between intact, body-temperature human cells and intact, room-temperature rat cells, obtaining Since the RUs kinetics may depend on both the species and the temperature, we re-calibrate the parameters k off and k basic on the base of the kinetic metrics reported in [17] (the data are referred to body-temperature human cells). These metrics include the peak tension T peak a , the time-to-peak T T P (defined as the time separating the beginning of force raise and the tension peak) and the relaxation times RT 50 and RT 95 (defined as the time between the tension peak and 50% and 95% of relaxation, respectively). Since, as to the best of our knowledge, detailed calcium transients measurements for intact human cells at body temperature are not currently available, we employ the synthetic calcium transient predicted by the ToR-ORd model [18]. The metrics reported in [17] are referred to different values of SL, expressed as fraction of the optimal sarcomere length (i.e. the length for which an increase of developed force is compensated by an increase in resting tension), corresponding, according to the authors, to nearly SL = 2.2 µm. For the calibration, we employ the value associated with 95% of the optimal length (i.e. 2.09 µm).
Finally, for the calibration of the parameters ruling the XBs cycling, we use the same values of Tab. S2-1. Therefore, since the calibration depends on the parameters previously set for the RUs activation, the resulting values of the parameters are slightly different. We provide in Tab. 3 the full list of parameters for both species (room-temperature rat and body-temperature human) and for both models (SE-ODE and MF-ODE).