Multiscale Modeling of Cardiovascular Function Predicts That the End-Systolic Pressure Volume Relationship Can Be Targeted via Multiple Therapeutic Strategies

Most patients who develop heart failure are unable to elevate their cardiac output on demand due to impaired contractility and/or reduced ventricular filling. Despite decades of research, few effective therapies for heart failure have been developed. In part, this may reflect the difficulty of predicting how perturbations to molecular-level mechanisms that are induced by drugs will scale up to modulate system-level properties such as blood pressure. Computer modeling might help with this process and thereby accelerate the development of better therapies for heart failure. This manuscript presents a new multiscale model that uses a single contractile element to drive an idealized ventricle that pumps blood around a closed circulation. The contractile element was formed by linking an existing model of dynamically coupled myofilaments with a well-established model of myocyte electrophysiology. The resulting framework spans from molecular-level events (including opening of ion channels and transitions between different myosin states) to properties such as ejection fraction that can be measured in patients. Initial calculations showed that the model reproduces many aspects of normal cardiovascular physiology including, for example, pressure-volume loops. Subsequent sensitivity tests then quantified how each model parameter influenced a range of system level properties. The first key finding was that the End Systolic Pressure Volume Relationship, a classic index of cardiac contractility, was ∼50% more sensitive to parameter changes than any other system-level property. The second important result was that parameters that primarily affect ventricular filling, such as passive stiffness and Ca2+ reuptake via sarco/endoplasmic reticulum Ca2+-ATPase (SERCA), also have a major impact on systolic properties including stroke work, myosin ATPase, and maximum ventricular pressure. These results reinforce the impact of diastolic function on ventricular performance and identify the End Systolic Pressure Volume Relationship as a particularly sensitive system-level property that can be targeted using multiple therapeutic strategies.


Additional details relating to the MyoSim contraction model
The MyoSim model used in this work was described in detail by Campbell et al. (Campbell et al., 2018). Key details are reproduced in this supplement to enhance clarity.
As shown in Fig S1, binding sites on the thin filament transitioned between an inactive state termed N off (which myosin heads could not attach to) and an active state N on (which was available for myosin binding). Activated binding sites could not switch back to the inactive state if a myosin head was attached. Myosin heads transitioned between an OFF state (that could not interact with actin), an ON state (that could potentially bind to actin), and a single attached force-generating state. Sites on the thin filament switch between states that are available (N on ) and unavailable (N off ) for cross-bridges to bind to. Myosin heads transition between an OFF detached state, an ON detached state, and a single attached force-generating state. J terms indicate fluxes between different states.

Thin filament transitions
The fraction of binding sites in the active state was defined as N on. The flux for the N off to N on transition (that is, the number of sites per unit time switching from N off to N on ) was defined as where k on is a rate constant, N overlap is the fraction of binding sites that are in range of myosin heads and depends on the prevailing half-sarcomere length as described in Campbell (Campbell, 2009), and k coop is a constant that defines the strength of thin filament cooperativity.
The flux of binding sites through the off transition was defined as where k off is a rate constant and the other terms are as defined previously. The (N 0 -N bound ) term means that only unbound sites can deactivate (that is, myosin heads have to detach before a binding site can revert to the off state). The k coop ((N overlap -N on )/N overlap ) term causes binding sites in the off state to induce further deactivation.

Thick filament transitions
The flux of myosin heads transitioning from the M OFF state to the M ON configuration was defined as where k 1 is a rate constant, k force is a parameter with units of N -1 m 2 , F total is the force in the muscle, and M OFF is the proportion of myosin heads in the OFF state.
The flux of myosin heads into the OFF state was defined as where k 2 is a rate constant and M ON is the proportion of myosin heads in the ON state.
Myosin heads bound to actin with a flux J 3 defined as where k 3 is a rate constant, k cb is the stiffness of the cross-bridge link, k B is Boltzmann's constant (1.38×10 -23 J K -1 ) and T is the temperature. Equation 5 mimics a second order reaction so that the rate at which myosin heads attach to the thin filament increases with the product of the proportion of myosin heads that are able to attach (M ON ) and the fraction of available binding sites (N on -N bound ). The Gaussian form of the myosin straindependence reflects the probability of the cross-bridge spring being extended to x by random Brownian motion when the myosin head binds to actin (Campbell and Lakie, 1998).
Similarly, the flux through the myosin detachment step was defined as where k 4,0 is a rate constant, k 4,1 is a parameter that sets the strain dependence of the cross-bridge detachment rate, x ps is the power-stroke of an attached cross-bridge, and M FG (x) is the proportion of cross-bridges attached to actin with spring-lengths between x and x+δx. (S7) Cross-bridge populations were evaluated with 0.5 nm resolution over the range -10 nm ≤ x ≤ 10 nm. n was thus equal to 41 giving a complete set of 45 equations. These were integrated numerically with adaptive step-size control. Calculations were initiated with all binding sites in the N off configuration and all myosin heads in the M OFF state.
Dynamic interfilamentary movement was incorporated into simulations when required by using polynomial interpolation to displace the distribution that described the number of heads bound with each spring length (Campbell, 2014). Filament compliance effects (Huxley et al., 1994;Wakabayashi et al., 1994) were mimicked by assuming that a half-sarcomere length change of Δx displaced each actin-myosin link by ½.Δx (Getz et al., 1998;Campbell, 2009 and N 0 is the number of myosin heads in a hypothetical cardiac half-sarcomere with a cross-sectional area of 1 m 2 . N 0 was set to 6.9×10 16 m -2 throughout this work based on the assumptions that (a) myofibrils occupy ~60% of the cross-sectional area of myocardium, (b) half-thick filaments contain 283 myosin heads, and (c) half-thick filaments have a spatial density of 4.07×10 14 m -2 within myofibrils (Linari et al., 2007;Campbell, 2014).
Passive force was calculated as where x hs is the length of the half-sarcomere, and σ, L slack ,, and L are model parameters. In this equation, σ is a scaling parameter, L slack defines the length of the half-sarcomere at which passive force is zero, and L sets the curvature of the passive-force / length relationship. These properties reflect collagen content and the isoform and posttranslational status of titin which change with disease in humans (Hidalgo and Granzier, 2013;LeWinter and Granzier, 2014).
Myosin ATPase was calculated as where ΔG` is the free energy produced by ATP hydrolysis (70 kJ mol -1 ), L 0 is the reference length of halfsarcomere (1.1 µm), and N A is Avogadro's number (6.02 x 10 23 mol -1 ). Ventricular efficiency was calculated as stroke work divided by myosin ATPase and omits the energy required to maintain the ionic gradients.

Circulation
Equation S11 defines the rate of change of the volume of each compartment in the circulatory model. Each term is simply the difference between the blood flows into and out of the compartment.
Equation S12 defines the inter-compartmental flows. The aortic and mitral valves were simulated using simple condition statements based on the pressure gradient. Other flows were defined by Ohm's law.
References for supplementary information Campbell, K.S. (2009). Interactions between connected half-sarcomeres produce emergent behavior in a mathematical model of muscle.      Table S1.

PLOS
The red lines shows the best-fit of a 5 th order polynomial to the simulated data.
Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of k 3 (equation S5) ranging from 0.1 to 10 times the value shown in Table S1.
The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of k cb (equation S5) ranging from 0.1 to 10 times the value shown in Table S1.
The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of k force (equation S3) ranging from 0.1 to 10 times the value shown in Table S1.
The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of k on (equation S1) ranging from 0.1 to 10 times the value shown in Table S1.
The red lines shows the best-fit of a 5 th order polynomial to the simulated data.
Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of k 4,0 (equation S6) ranging from 0.1 to 10 times the value shown in Table S1.
The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of σ (equation S9) ranging from 0.1 to 10 times the value shown in Table S1.
The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of k coop (equations S1 and S2) ranging from 0.1 to 10 times the value shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data.
Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of k 4,1 (equation S6) ranging from 0.1 to 10 times the value shown in Table S1.
The red lines shows the best-fit of a 5 th order polynomial to the simulated data.      Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of g to ranging from 0.1 to 10 times the base value in ten Tusscher et al's electrophysiological model (ten Tusscher et al., 2004). The red lines shows the best-fit of a 5 th order polynomial to the simulated data.  Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data.  Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data.  Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data.  Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of C aorta (see Methods in main text) ranging from 0.1 to 10 times the value shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of R capillaries (see Methods in main text) ranging from 0.1 to 10 times the value shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of R aorta (see Methods in main text) ranging from 0.1 to 10 times the value shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of R ventricle (see Methods in main text) ranging from 0.1 to 10 times the value shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of R arteries (see Methods in main text) ranging from 0.1 to 10 times the value shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of R veins (see Methods in main text) ranging from 0.1 to 10 times the value shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of C arteries (see Methods in main text) ranging from 0.1 to 10 times the value shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of C capillaries (see Methods in main text) ranging from 0.1 to 10 times the value shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for values of C arterioles (see Methods in main text) ranging from 0.1 to 10 times the value shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data. Panels A to L show values (blue circles) for 12 system-level properties (for example, maximum ventricular pressure) predicted for different combinations of σ and L (see equation S9) so that passive stress at a halfsarcomere length of 1069 nm was held constant. Other values as shown in Table S1. The red lines shows the best-fit of a 5 th order polynomial to the simulated data.