Understanding the dynamics of a cardiac muscle twitch contraction is complex because it requires a detailed understanding of the kinetic processes of the Ca2+ transient, thin-filament activation, and the myosin–actin cross-bridge chemomechanical cycle. Each of these steps has been well defined individually, but understanding how all three of the processes operate in combination is a far more complex problem. Computational modeling has the potential to provide detailed insight into each of these processes, how the dynamics of each process affect the complexity of contractile behavior, and how perturbations such as mutations in sarcomere proteins affect the complex interactions of all of these processes. The mechanisms involved in relaxation of tension during a cardiac twitch have been particularly difficult to discern due to nonhomogeneous sarcomere lengthening during relaxation. Here we use the multiscale MUSICO platform to model rat trabecular twitches. Validation of computational models is dependent on being able to simulate different experimental datasets, but there has been a paucity of data that can provide all of the required parameters in a single experiment, such as simultaneous measurements of force, intracellular Ca2+ transients, and sarcomere length dynamics. In this study, we used data from different studies collected under similar experimental conditions to provide information for all the required parameters. Our simulations established that twitches either in an isometric sarcomere or in fixed-length, multiple-sarcomere trabeculae replicate the experimental observations if models incorporate a length–tension relationship for the nonlinear series elasticity of muscle preparations and a scheme for thick-filament regulation. The thick-filament regulation assumes an off state in which myosin heads are parked onto the thick-filament backbone and are unable to interact with actin, a state analogous to the super-relaxed state. Including these two mechanisms provided simulations that accurately predict twitch contractions over a range of different conditions.
In cardiac muscle, transient increases in the cytosolic calcium concentration trigger twitchlike contractions during the systolic phase of heart function (Gordon et al., 2000). Ca2+ binding to cardiac troponin (cTn) results in movement of tropomyosin (Tpm) over the F-actin surface of the thin filament. This exposes myosin-binding sites on actin that allow for myosin–actin cross-bridge formation, cyclic interactions of myosin heads with actin, force generation, and contraction. After reaching a peak, intracellular Ca2+ concentration ([Ca2+]) is rapidly removed from the cytosol, lowering [Ca2+], enhancing additional dissociation of Ca2+ from cTn, and allowing Tpm to return to its blocking position. The increasing number of Tpms in the blocked position reduces overall binding flux, reduces the number of bound cross-bridges, and results in muscle relaxation. Understanding the dynamics of this process in a single twitch contraction is complicated because it requires a detailed understanding of the kinetics of each step of the process: the Ca2+ transient, the dynamics of the troponin (Tn)–Tpm switch, and the myosin cross-bridge chemomechanical cycle. Each of these steps is well defined individually, but how all three processes work in combination is not understood. Here we build on our previous work on modeling the different elements of the contraction cycle to define a novel multiscale model of a cardiac muscle twitch contraction.
Computational models have the potential to interconnect multiple types of experiments and protocols and can be a useful tool for generating and testing hypotheses regarding the mechanisms of contractile function in healthy hearts and contractile dysfunction in diseased hearts. Early models defined structural cooperative effects in sarcomeres, such as end-to-end interactions between adjacent Tpm–Tn regulatory units, between cross-bridges bound to actin and the thin-filament regulatory units, and between neighboring bound cross-bridges (Campbell et al., 2001; Rice et al., 2008; Tanner et al., 2012). Several of these models are derived from the original formulation by Razumova et al. (1999, 2000), where the muscle tension is proportional to the product of fractional occupancy of the myosin bound states and the mean cross-bridge distortion of each state population. The approximation of strain-dependent distributions by mean distortions is mathematically effective because the model can be defined by a system of ordinary differential equations. The newer generations of spatially explicit models, coupled with a cross-bridge cycle and thin-filament regulation by Ca2+, show significantly improved relationships between molecular interactions and muscle fiber behavior (Chan et al., 2000; Tanner et al., 2012; Mijailovich et al., 2016; Mijailovich et al., 2019). The key feature of these models is the ability to predict mechanical response to Ca2+ transients reflecting cardiac muscle behavior in a living heart. The dynamic relationships in a contraction are defined by the coupling of thin-filament regulation by Ca2+ and the strain-dependent kinetics of myosin binding to actin. However, just as with experiments studying force generation and cross-bridge cycling, this type of model has been limited primarily to simulating force generation under steady-state Ca2+ conditions (Chase et al., 2004; Kataoka et al., 2007; Tanner et al., 2007; Tanner et al., 2008; Tanner et al., 2012). Modeling twitch contractions has been approached with different levels of complexity, ranging from models based on solution of ordinary differential equations (Campbell et al., 2001; Niederer et al., 2006; Rice et al., 2008; Land et al., 2017), to a spatially detailed cardiac sarcomere model (Washio et al., 2012), to a continuous time Markov (regulatory) chain coupled with a generalized Huxley 57 model (Regazzoni et al., 2020), to more detailed models based on solution of partial differential equations (Chung et al., 2017; Campbell et al., 2018).
Moreover, the latest models include mechanisms that could alter cardiac muscle contractility by structural changes in thick filaments. Possible mechanisms could include “mechanosensing” in thick filaments (Linari et al., 2015; Irving, 2017) and interfilament mechanical signaling by myosin-binding protein C (MyBP-C; Kampourakis et al., 2016; Irving, 2017; Brunello et al., 2020). Following these ideas, Campbell et al. (2018) hypothesized that the transition rate from the parked state (PS) to the myosin on state, which permits binding to actin, is dependent on force. They successfully fitted the length-dependent activation data, where the twitch contractions only include fitting of the observed tension traces, but they did not take into account the observed significant shortening and lengthening at the sarcomere level known to occur in trabeculae held at fixed length (Janssen and de Tombe, 1997; Ferrantini et al., 2014; Caremani et al., 2016). In addition, these fits were achieved in the model using simplified cross-bridge cycle kinetics and sarcomere geometry. The simulations of twitch contractions, in which transition rate from PS state to M.D.Pi (myosin head with nucleotides ADP and phosphate, Pi) are force dependent, using a more detailed model (muscle simulation code, MUSICO), however, showed that the relaxation of tension is slow at higher levels of tension even when the [Ca2+] falls to low values. Our preliminary results suggest that the transition rate from PS state should be defined by the relationship of some process other than a linear dependence on tension because the PS should act as a sink for myosin heads during relaxation and a source for additional myosin head recruitment during contractions (Prodanovic et al., 2020).
Taken together, all of these models can more or less successfully replicate observed tension traces in response to Ca2+ transients. However, they all include a large number of simplifications, so that changes in molecular interactions within sarcomeres, such as by mutations of contractile and regulatory proteins, may not be transparent across the multiple length scales. Thus, these models may not be able to translate the changes at the molecular level to the observed changes in cardiac muscle behavior in associated diseases.
Furthermore, a major challenge in validating computational models by simulating experimental data has been the variability in conditions used by investigators. For example, some experimenters have recorded tension and sarcomere length (SL) transients but not the Ca2+ transient (Caremani et al., 2016), or they have recorded Ca2+ and tension transients at different temperatures without monitoring changes in the SL during the twitch in fixed-length trabeculae (Janssen et al., 2002). The only studies that we are aware of that reported [Ca2+], tension, and SL simultaneously were reported by Janssen and de Tombe (1997) and Ferrantini et al. (2014). Studies that do not include SL control do not account for the series elastic (SE) compliance, i.e., for significant changes in SLs during contractions of intact muscle preparations held at a fixed length. Attempts to take this SE compliance into account include the study by Kataoka et al. (2007), where the series compliance of trabeculae was exclusively attributed to the compliance of thin and thick filaments but did not account for SL changes caused by SE compliance of trabeculae. This series compliance was considered by Rice et al. (2008), who demonstrated a profound effect on the tension response during twitch contractions. However, Rice at al. used a model that showed only the sensitivity of the responses to a large range of simulated compliances, with little reference to experimental data. To more precisely simulate the steady-state and dynamic properties of striated muscle contraction and relaxation, all compliances and SL changes must be considered because cross-bridge kinetics is sensitive to the changing load as the sarcomeres shorten or lengthen.
In this study, we report a model that more precisely defines the multiple sarcomere parameters that vary during contraction and relaxation, using the methodology implemented in the MUSICO computational platform (Mijailovich et al., 2016; Mijailovich et al., 2019). This model allows simulation of the tension developed during a twitch contraction, as well as the kinetics of contraction and relaxation. Previous models have successfully simulated the kinetics but have accomplished this using a normalized force value. We initially simulated the twitch data of Janssen and de Tombe, 1997; Caremani et al., 2016 without normalizing tension using a five-state cross-bridge model. This resulted in an elevation of resting force and prolonged relaxation that were addressed by addition of two features: (1) SE compliance to account for shortening and lengthening of sarcomeres during a twitch in a trabecula kept at fixed length and (2) expansion of the cross-bridge cycle to include a myosin “off state” that acts as a source for myosin recruitment during activation of contraction and as a sink for taking myosin out of the cycling pool during relaxation (analogous to the super-relaxed [SRX] state of the thick filament; Linari et al., 2015; McNamara et al., 2015; Irving, 2017). Addition of this state to the cross-bridge model accelerated relaxation kinetics and reduced resting tension, producing simulated twitches that matched the experimental twitch data of Janssen and de Tombe, 1997; Caremani et al., 2016. Additional validation of the model was performed by simulations of a twitch for multiple datasets, including varying temperature, and the force–Ca2+ relationship. Thus, the model can predict both twitch transients and the Ca2+ sensitivity of force. Preliminary reports of this work have been published previously in the forms of abstracts and proceedings (Prodanovic et al., 2020; Mijailovich et al. 2018. Biophys. J. 114:Abstract 500A. https://doi.org/10.1016/j.bpj.2017.11.2737; Mijailovic et al. 2019. Biophys. J. 116:Abstract 116A. https://doi.org/10.1016/j.bpj.2018.11.654).
Materials and methods
Here we outline the essential definitions and relationships developed in previous studies, including (1) the multiscale structural organization of intact muscle (Mijailovich et al., 2016) and (2) Ca2+-dependent kinetics that regulate Tn–Tpm interaction with F-actin in cardiac muscle (Geeves et al., 2011; Mijailovich et al., 2012; Mijailovich et al., 2019; Prodanovic et al., 2020). We also describe new features implemented in MUSICO to model transient contractions in cardiac trabeculae, namely (3) an expanded cross-bridge cycle that contains five essential states and further expansion to six states, (4) intrinsic series compliance that permits SL changes during twitch in fixed-length trabeculae, and (5) and thick-filament activation.
Multiscale model of trabeculae
Intact ventricular trabeculae from the rat heart are composed of parallel cardiac muscle cells, or myocytes, embedded in an extracellular matrix that includes the collagen structure laterally connecting the muscle fibers. Myocytes occupy ∼70% of a trabecular volume, and the cardiac cells have a diameter ranging between 10 and 20 µm and a length between 80 and 100 µm. The myocyte contractile substructure consists of myofibrils aligned in parallel, where each myofibril is composed of 40–50 sarcomeres in a series. The sarcomeres are composed of interdigitated thick and thin filaments, where each half thick filament consists of ∼150 myosin molecules projecting on each side of the M-line and each thin filament contains from 360 to 400 actin monomers emerging from Z-discs on opposite sides of a sarcomere. The thick filaments are associated with the auxiliary protein titin, connecting it via the thick-filament backbone with the Z-disc and the M-line, and MyBP-C connecting thick and thin filaments. The thin filaments contain the regulatory proteins Tpm and Tn essential for Ca2+ regulation of contraction and other auxiliary proteins regulating thin-filament length and other functions.
Trabecular geometry and elasticity
The contractile characteristics of cardiac muscles are typically derived from isometric and isotonic twitch contractions in ex vivo preparations. However, twitches of trabeculae that are isometric at the muscle level are not isometric at the level of the sarcomere (Janssen and de Tombe, 1997; Caremani et al., 2016). During twitches, when force rises, the sarcomeres shorten; during relaxation, they lengthen. The origin of the elasticity in trabeculae is complex and may include muscle cell alignment and connective tissue serial elasticity. The simplest method is to include a trabecula SE element in series with the contractile element (CE) consisting of multisarcomere muscle fibers (Fig. 1).
Spatially explicit 3-D model of a sarcomere
In the 3-D sarcomere, an overlapping hexagonal lattice of thick and thin filaments has a well-defined interfilament spacing where the mean spacing is measured by x-ray diffraction (Higuchi et al., 1995; Millman, 1998). Explicit formulation of the sarcomere lattice is defined by the discrete position of two-headed myosin molecules and their binding sites on actin filaments. For simplicity, we consider that only one myosin head per myosin molecule is active at any one time and is denoted as a cross-bridge. Overall, the 3-D sarcomere structure is viewed as an array of thin and thick filaments connected by cross-bridges and other elastic elements (e.g., titin and MyBP-C) in a lattice network, and all of these elements can be represented as linear springs and elastic beams (Daniel et al., 1998; Mijailovich et al., 2016). The instantaneous equilibrium between actin and myosin filaments interconnected by cross-bridges is defined by a stiffness matrix, which includes the elasticity of thick and thin filaments, attached cross-bridges and titin, and load containing a vector of all external forces and internal forces generated by the action of cross-bridges. The loading conditions are defined at the boundaries of the muscle system, and they follow the experimental protocol of the relevant experiment. Typically, these include changes of length or imposed tension at the muscle end(s). The stiffness matrix is constantly changing as actomyosin bonds are created or disrupted, and active forces are also changing during conformational changes in attached myosins.
In this study, we used experimental conditions such as isometric (half) sarcomeres and fixed-length trabeculae. The latter case includes variable loading (tension) on multisarcomere (muscle fiber) structure due to serial elasticity in fixed-length trabeculae. Regarding the MUSICO simulations of trabeculae, it is possible to prescribe at the end of trabeculae arbitrary changes of length, such as in previous studies (Janssen and de Tombe, 1997; Caremani et al., 2016) where trabecular length was changed in order to keep SL constant or approximately constant or tension was applied on the trabecular ends. This kind of loading is shown in Fig. S3. The solutions of the equations defining the mechanical system of sarcomeres in series are obtained by standard finite element procedures for nonlinear systems using an iterative procedure (Bathe and Mijailovich, 1988; Mijailovich et al., 1993; Bathe, 1996).
Strain-dependent kinetics of cross-bridge cycle in the 3-D sarcomere lattice
The discrete nature of the myosin-binding sites in two helically arranged strands in each actin filament and the regular arrangement of myosin crowns along the myosin filaments require that myosin heads in a 3-D sarcomere lattice move not only axially but also radially and azimuthally to reach binding sites on actin monomers with the correct orientation (Squire, 1992). The relative spatial positions between myosin heads and binding sites on the thin filament are obtained from the current spatial positions of the sites on deformed thin and thick filaments (Mijailovich et al., 2016; Mijailovich et al., 2019), as shown in Fig. S1. The cause of deformation in the elastic filaments originates from local forces imposed by bound cross-bridges and other elastic elements (e.g., titin), but, in turn, the forces in cross-bridges depend on state and the relative positions between actin and myosin sites. The strain-dependent rate of myosin binding to actin in the 3-D sarcomere lattice is modulated by the axial strain component, as used in all cross-bridge kinetic models, and by radial and azimuthal spatial arrangement relative to binding sites on the actin filaments (Fig. S1), displayed by interfilament spacings, and azimuthal angles α and β (Mijailovich et al., 2016, 2019).
Minimal cross-bridge cycle including five essential states
Simple kinetic models include only three states (e.g., Daniel et al., 1998; Duke, 1999; Mijailovich et al., 2016; Mijailovich et al., 2019), where multiple states are lumped into a single state. A cross-bridge model with additional states is desirable to relate the underlying well-defined molecular kinetics to the simulations of, for example, twitches in cardiac trabecula and allow simulation of tension levels in units of the observations. The additional states provide closer proximity to measurements of state transition rates from assays in solution and motility assays and thus provide experimentally measured model parameters at the molecular scale rather than parameters estimated from fits of data from higher length scale experiments. Hence, we expanded the above models into a five-state cross-bridge model that introduces a separate strain-dependent ADP release to distinguish it from ATP binding associated with cross-bridge detachment and the hydrolysis step that separates M.T (myosin head with ATP) and M.D.Pi states (Fig. 2 A). The former defines the maximum cycling rate of the cross-bridge cycle as a function of load. The hydrolysis step defines the minimum lifetime of the detached cross-bridge. Separation of M.T and M.D.Pi is convenient for introduction of the concepts associated with thick-filament regulation defined as interaction of the detached M.D.Pi state and the PS (Fig. 2 B).
Strain-dependent, thermally activated chemical reaction rates are obtained from transition state theory (Eyring, 1935; Glasstone et al., 1941) as updated by Kramers (1940), assuming that the rate of a reaction is limited by the energy barrier presented by the point of highest Gibbs energy on the reaction path. The general form of the state transition rates in models with a power stroke is formulated by Hill (1974) as the ratio of forward to backward rates that must satisfy Gibbs’s thermodynamic identity. The rates are defined in terms of the Gibbs energies of the initial and final states, including the elastic strain energy derived from the cross-bridge tension (Hill, 1974; Wood and Mann, 1981).
The state transition rates between states in the cross-bridge cycle are derived from the free energy of the cross-bridge states in terms of the axial strain component The complete set of the state transition rate constants is described in Appendix A. Moreover, in the 3-D sarcomere lattice, the strain-dependent rate of myosin binding to actin is modulated, in addition to the axial strain component by the spatial position of the actin filaments and the azimuthal departure of a myosin-binding site on an actin filament from the plane passing through actin longitudinal axes (Fig. S1). This effect is implemented in MUSICO via weight factors and for the azimuthal departures and, in addition, includes the normalization factor, that takes into account the number of actin sites that a cross-bridge can reach and bind (for further details, see Mijailovich et al., 2016; Mijailovich et al., 2019).
Thin-filament regulation by calcium
Muscle contraction and relaxation are regulated by the Ca2+-dependent azimuthal movements of Tpm–Tn complexes over the surface of the actin filament. Structurally, Tpm covers seven monomers on one strand of the actin double helix and is associated with one Tn molecular complex. The Tn complex consists of troponin T, troponin C, and troponin I, denoted as TnT, TnC, and TnI, respectively. The affinity of myosin for the actin filament is controlled by the azimuthal position of Tpm, where Tpm molecules are aligned along each strand of the actin double-stranded helix. Most regulatory models assume that each Tpm molecule moves between three discrete orientations: blocked, closed, and myosin induced, originally called the “open state” (McKillop and Geeves, 1993; Lehman et al., 2000; Pirani et al., 2005; Poole et al., 2006). These Tpm positions are defined by TnI bound to actin, a free chain interacting with actin, and a chain displaced by strongly bound myosin heads, respectively. Calcium regulation of these processes is defined by calcium binding or dissociation from cardiac TnC and change in the affinity of TnI to actin surface induced by a conformational change in the Tn molecule.
Calcium-dependent kinetics of Tpm–Tn complex with actin in cardiac muscle
Ca2+ regulation of the interaction of the Tpm–Tn complex with actin is defined by open and closed allosteric states of TnC with one Ca2+ ion bound in cardiac muscle or with two Ca2+ ions in skeletal muscle (McKay et al., 2000). Ca2+ binding to TnC exposes a hydrophobic region in the N-domain of TnC, providing a binding site for the segment of TnI, strengthening the interaction of TnC with TnI, and weakening the affinity of TnI to actin. A minimal description of the allosteric mechanism in cardiac muscle includes four states (Fig. 2 C), comprising two TnC closed states, where TnC has no or one bound Ca2+ to N-terminal site II, denoted as TnC and CaTnC, and similarly two TnC open states. In the presence of Ca2+, the open state is favored over the closed state, and TnI binds preferentially to the open TnC, forming CaTnC.TnI state. In the absence of Ca2+, the closed form of TnC is favored, and TnI binds preferentially to actin, forming A.TnI state. The equilibrium rate constant of Ca2+ binding to TnC closed states, is effectively defined via the forward constant, where linearly depends on and the Ca2+-independent dissociation constant, The dissociation constant of TnI from closed A.TnI to CaTnC.TnI is defined by a first-order affinity λ when one Ca2+ is bound. With no bound Ca2+, the transition from A.TnI to TnC.TnI shows slow dissociation of TnI from actin with an attenuated rate, where allosteric fraction is On the other hand, Ca2+ binding to TnC in a TnC.TnI complex is accelerated to keeping the population of the A+TnC.TnI state low.
Thin-filament regulation of myosin interactions with actin
The kinetics of interactions between the Tpm–Tn complexes with an actin filament are, at present, best described by the long-range cooperative continuous flexible chain (CFC) model (Geeves et al., 2011; Mijailovich et al., 2012) that includes structural constraints between Tpm–Tn regulatory units. The model is built on structural evidence that neighboring Tpms overlap and that one end of TnT is bound to a specific site on Tpm, whereas its N terminus overlaps with the adjacent Tpm, forming linked Tpm–Tpm regions. The interconnected neighboring Tpm–Tn units form the appearance of a CFC on each strand of an actin (Lorenz et al., 1993; Vibert et al., 1997) rather than a set of the independent Tpm–Tn units.
For modeling thin-filament regulation as a CFC, we follow a Monte Carlo approach with spatially explicit myosin binding to regulated actin filaments in solution (Geeves et al., 2011; Mijailovich et al., 2012) that we later expanded to the model of thin filaments in a 3-D sarcomere lattice (Mijailovich et al., 2019). The coupled algorithm for thin-filament regulation and the actomyosin cycle consists of three main steps: (1) calculation of the state transitions between actin–Tpm–Tn states as a function of [Ca2+], (2) calculation of mean CFC angular positions and their azimuthal angular fluctuations along the actin filament, and (3) calculation of state transitions between actin–myosin states. These processes are interrelated and regulated by the [Ca2+] where CFC angular position is modulated by either TnI bound to actin or dissociated from actin after calcium binding to TnC. The release of TnI from the actin surface allows Tpm angular movement, exposing myosin-binding sites on actin and thus permitting myosin binding and force generation by cross-bridges. The CFC moves azimuthally at a high frequency, and the mean CFC angular position and the thermally induced azimuthal fluctuations are defined by the instantaneous spatial positions of the pinning sites, at which the TnI and myosin are bound to actin (Geeves et al., 2011; Mijailovich et al., 2012; Mijailovich et al., 2019). General definitions, calculations of CFC angular position and its variance, and the coupling between the calcium-regulated position of the CFC and myosin cycle states are described in Appendix B.
Thick-filament regulation by calcium
There is now substantial evidence for a myosin off state in the relaxed structure of the thick filament in which the two myosin heads self-associate and pack down onto the thick filament to form what is referred to here as the “PS.” This state is denoted as a disordered relaxed state (DRX) associated with the thick-filament backbone (Anderson et al., 2018; Brunello et al., 2020) or an ordered SRX state (McNamara et al., 2016). The kinetics of the equilibrium between the PS and the on state of myosin on the thick filament is necessary to define how the filament is activated and relaxed and how the transition between these thick filament states is modified by mechanosensing in thick filaments (Linari et al., 2015; Irving, 2017), interfilament mechanical signaling by cardiac MyBP-C (McNamara et al., 2016; Brunello et al., 2020), and phosphorylation of the regulatory light chain (RLC) and cardiac MyBP-C (Kampourakis et al., 2016) and, perhaps, titin.
In the absence of precise information concerning the molecular details of each of these prospective mechanisms, we began by implementing a minimal version of the PS. We treated each myosin head as an independent force generator that can isomerize from the M.D.Pi state to the PS. For reducing high levels of tension at low [Ca2+] (Prodanovic et al., 2020), we modeled the PS to be activated by Ca2+ association with a component of the thick filament so that the activation follows the Ca2+ transient. This idea is supported by the observation of Huxley et al. (1994) that, during fast Ca2+ activation, myosin heads move radially away from the thick-filament backbone much faster than they attach to the thin filament. Thus, this process could be attributed to a [Ca2+]-dependent transition from the PS to M.D.Pi state that can reduce resting tension at low levels of Ca2+ and accelerate tension relaxation when [Ca2+] falls. Alternative hypotheses are considered in the Discussion section.
Monte Carlo simulations of rate-dependent stochastic processes
In Monte Carlo simulations, we used the standard Metropolis algorithm, where a kinetic transition in a time step occurs when a random number in (0, 1) lies in the range (0, where k is state transition rate constant. must be much smaller than the inverse of the fastest rate constant of the system, and, in practice, was sufficient to achieve satisfactory statistics and avoid interference between multiple transitions within a single subsystem and the negligibly small interference between the subsystems.
The coupling between thin-filament regulatory processes and the actomyosin cycle requires two sets of Monte Carlo random number drawings within each time step. The first set of random number drawings defines transitions of TnI–actin states from the calcium kinetics scheme in Fig. 2 C, depending on [Ca2+]. The position of bound TnI to actin after updating and bound myosins from the previous time step set the chain configuration for obtaining the mean CFC angular position along the actin filament strand and its azimuthal variation. The second set of random number drawings defines the changes in actomyosin states regulated by the CFC, setting the configuration for the calculation of the mechanical equilibrium with external forces and constraints.
For each TnI or cross-bridge, we use one Monte Carlo random number drawing to define whether the TnI or cross-bridge remains in its current state or will change its state into one of the other possible states within the current time step The change of state is defined by the probability, in the range from 0 to 1, that is divided into probability bins, in a specified order, including the set of probability bins associated with a TnI or cross-bridge changing state and a bin associated with the probability of remaining in the current state. Depending which bin the drawn random number falls into, the fate of a particular TnI or cross-bridge is defined and set for the following time step. The calculations of the transition probabilities from the transition rates between TnI–actin states and between actomyosin states in the 3-D sarcomere lattice are as described previously (Mijailovich et al., 2012; Mijailovich et al., 2016; Mijailovich et al., 2019). The coupling between cross-bridge cycling and calcium regulation of the thin filament is described in Appendix C.
Stochastic process, model size, and myofibril edge effects
Due to stochastic transitions in the cross-bridge cycle, the force per myosin filament can vary between the filaments by ±15% at any instant of time and also by about the same amount in the same filament over time. In contrast, the variation in average force per myosin filament is small; for example, if the number of myosin filaments matches the number of filaments in an average myofibril, the fluctuations in isometric force are minimal (<0.1%) and reflect observed tension variations.
In the parameter exploration phase, we limited the stochastic simulations to a half-sarcomere with 200 myosin and 400 actin filaments that is approximately one-third of the number of filaments in a cross-section of a typical myofibril. This number of filaments is sufficient for stable simulations and statistical averaging without requiring running the simulation multiple times. In the final simulations, we increased the number of myosin filaments per half-sarcomere to 500 in order to get realistically smooth tension responses.
Procedure for estimation of model parameters
The MUSICO platform is formulated on three principal pillars: the 3-D explicit (multi)sarcomere geometry, a strain-dependent cross-bridge cycle, and calcium-dependent thin- and thick-filament regulation. This model provides the precision necessary to quantify the effects of kinetics and the structural changes, associated with mutations in contractile and regulatory proteins, on muscle function, but it requires a large number of model parameters to properly define muscle fiber structure, cross-bridge kinetics, and thin-filament regulation by Ca2+. We contend that conventional fitting of the data is not ideal, because the published experimental data show that large variations and fitting many free parameters can result in unreasonable estimates. Thus, we used established structural data and best estimates deduced from experimentally established kinetic data to define as many parameters as possible (Table 1). In this study, we then varied only a few parameters depending on the particular sets of data that model prediction had to match. In most of the cases, we adjusted (1) the myosin–actin binding rate, to account for changes of interfilament spacing, high speed of sarcomere shortening, and possible myosin activation by stretching thick filaments via titin at longer SLs; (2) the power stroke Gibbs free energy change, and the ADP release rate, because these parameters are load dependent and the degree of shortening and or sarcomere disorder is difficult to define, particularly during relaxation; (3) the kinetics of Ca2+ binding, and dissociation, from TnC that reciprocally defines interaction of TnI with actin; and (4) the [Ca2+] for the half-maximal rate from PS, [Ca2+]50, a key parameter defining the state transition rate from PS, due to large changes in the cytosol from intact cells to the bathing solution in demembranated muscles. The parameters used in simulations are described in detail in Appendix D, and the complete set of parameters used in simulations of trabecular twitch responses to calcium transients is shown in Table 1. In a few simulations, different parameters were used to accommodate large differences in experimental conditions, such as a large change in interfilament spacing, due to demembranation of trabeculae, or to achieve good fits of the relaxation phase of the twitch contractions that must account for inhomogeneous lengthening between sarcomeres.
A sensitivity analysis is typically needed when fitted parameters have a large range of possible values (Mijailovich et al., 2010; Mijailovich et al., 2012; Ujfalusi et al., 2018). In this study, we also performed a sensitivity analysis of the parameters and and b in which were necessary to be estimated (usually no more than three, depending on specific case) to match all experimental data. The sensitivity analysis of these and other important parameters is shown in the supplemental text (see bottom of the PDF and Figs. S8–S13 and S16–S23). We used the sensitivity analysis of a set of free parameters to estimate parameters for the next simulation by methodology described by Mijailovich et al. (2010). From those estimates, the changes in parameters were filtered, using our experience, in order to avoid imposing changes in parameters that marginally affect the predicted response and therefore move the next model prediction more quickly toward the observation. Using this approach, we achieved good fits with a minimal number of long simulation runs.
The adjusted parameters that differ from data in Table 1 are shown in comparative tables and in the figure legends associated with the specific simulations where the adjustments were necessary.
The MUSICO software environment and simulation details
The MUSICO software was developed as a C++ object-oriented application that includes the LAPACK linear algebra package and the deal.II finite element library. In parts of MUSICO, where parallelization provides a significant reduction in computational time, OpenML is implemented. Typical run times for these simulations depend on the number of actin and myosin filaments. The simulation of 200 myosin filaments over 1 s with a time step of 1 µs requires ∼10 h on the AEG IS04-KG grid site, consisting of six nodes, each equipped with two AMD Opteron 6276 16-core processors and 96 GB RAM, totaling 192 processors.
Online supplemental material
The supplemental text includes (1) the spatial positions between myosin heads and binding sites on the thin filament in 3-D sarcomere lattice; (2) the correction of Caremani et al. (2016) observed tension under isometric half-sarcomere length due to the loss of trabecular length control during the sarcomere-lengthening phase of a twitch; (3) how the accelerated relaxation can be achieved by variation of parameters in the five-state model; (4) measurements of passive SE compliance during twitch in isometric trabeculae; (5) ATPase consumption during twitch contractions in fixed-length trabeculae; (6) sensitivity analysis of twitch tension to model parameters in fixed-length trabeculae; (7) estimation of the active tension increase due to decrease of titin passive tension during twitch contractions in fixed-length trabeculae where large changes in SL are observed (Caremani et al., 2016); and (8) sensitivity analysis of twitch tension transients to model parameters in the absence of a PS in isometric half-sarcomere. The following supplemental figures are available online in a ZIP package. Fig. S1 shows the interaction between myosin heads and actin filaments in three dimensions defined by the triple-helical arrangement of myosin molecules along the myosin filament and the double-helical arrangement of monomers (myosin-binding sites) along the actin filaments. Fig. S2 shows the effect of the five-state cross-bridge model rate on relaxation of twitch tension in isometric half-sarcomeres (IoSarc.). Fig. S3 shows the SL controlled isometric tension of a half-sarcomere. Fig. S4 shows the significance of the PS in twitch relaxation phase. Fig. S5 shows the comparative tension–displacement loops obtained from Janssen and de Tombe, 1997; Caremani et al., 2016. Fig. S6 shows a comparison of ATP consumption rate and ATPase during twitch contractions of isometric half-sarcomere and fixed-length trabeculae. Fig. S7 shows that ATPase is strongly correlated to sarcomere-shortening velocities in fixed-length trabeculae during force development (in SL per second). Fig. S8 shows the changes in twitch tension from Fig. 5 B by ±20% variation in binding rate. Fig. S9 shows the changes in twitch tension from Fig. 5 B by ±20% variation in detachment rate. Fig. S10 shows the changes in twitch tension from Fig. 5 B by ±20% variation in ADP release rate. Fig. S11 shows the changes in twitch tension from Fig. 5 B by ±40% variation in Hill coefficient of the rate of sigmoidal rise. Fig. S12 shows the changes in twitch tension from Fig. 5 B by ±20% variation in [Ca2+] at [Ca2+]50. Fig. S13 shows the changes in twitch tension from Fig. 5 B by ±40% variation in baseline rate from the PS. Fig. S14 shows active tension corrected for the decrease in titin force due to shortening of a half-sarcomere. Fig. S15 shows the effect of allosteric TnC–TnI–actin–Ca2+ interaction parameter, on resting tension at observed baseline calcium level and when [Ca2+] was reduced to 0.05 µM. Fig. S16 shows changes in twitch tension predicted by the five-state model for change in for up to ±50% from those used for Fig. 3, denoted as original Fig. S17 shows changes in twitch tension predicted by the five-state model for change in or up to ±50% from those used in Fig. 3, denoted as original Fig. S18 shows changes in twitch tension predicted by the five-state model for change in for up to ±50% from those used in Fig. 3. Fig. S19 shows changes in twitch tension predicted by the five-state model for change in for up to 24-fold from those used in Fig. 3. Fig. S20 shows changes in twitch tension at isometric half-sarcomere obtained from simulations by the five-state model. Fig. S21 shows changes in twitch tension relaxation rates with increase of predicted by the five-state model. Fig. S22 shows the best fit to the observations with the five-state model requiring the change in power stroke rate Fig. S23 shows the effect of changes in baseline calcium level on resting and twitch tension. Fig. S24 shows a comparison of tension–pCa relationships between the cross-bridge cycle models with or without PS (i.e., the six-state and five-state models).
Isometric fixed end twitch contractions in cardiac muscle result from Ca2+ transients that lead to transient tension responses and SL changes. In intact trabeculae, this involves the combined effects of multiple underlying mechanisms acting together. Thus, to test the effects of two novel features of our new model, namely thick filament activation during contraction and the series elasticity observed in cardiac muscle, we first address them separately. To establish the kinetic parameters of the Ca2+ activation scheme and the cross-bridge cycle, we first simulated twitches in a single isometric half-sarcomere. Then, using these parameters, we simulated twitches of a fixed-length trabecula using a multisarcomere system in series that contains an elastic component to predict changes in the SLs.
Why a PS is necessary for twitch contractions of an isometric sarcomere
We tested the ability of the five-state cross-bridge model (Fig. 2 A), coupled to a Ca2+ regulatory scheme (Fig. 2 C), to simulate twitch contractions in an isometric half-sarcomere using MUSICO. There are only a few published trabecular twitch experiments where SL was kept approximately constant (Janssen and de Tombe, 1997; Caremani et al., 2016). The experiments of Janssen and de Tombe recorded Ca2+ and tension transients at 22°C, but the data have an unusual flat top to the tension transient, which is not typical of the transients recorded in the literature and is hard to describe in simulations (see Fig. S4; and Niederer et al., 2006). Caremani et al. recorded tension but not Ca2+ at 27.2°C. To provide the Ca2+ transient for the Caremani et al. experiment, we adapted the calcium transient of Janssen et al. (2002) collected at similar temperature (i.e., 27.5°C) during twitch contractions in fixed-length trabecula experiments. The transients typically achieve a higher peak [Ca2+] in fixed-length trabeculae, where internal shortening can occur, than in isometric sarcomere twitches where feedback length is controlled (Janssen and de Tombe, 1997). We rescaled the adapted Ca2+ transient for the Caremani et al. trabecula to the isometric half-sarcomere using the ratio of peak Ca2+ recorded in both protocols by Janssen and de Tombe, 1997 at the lower temperature (for details, see Appendix D).
In the isometric half-sarcomere twitches at 27.2°C observed by Caremani et al. (2016), the length of half-sarcomeres was well controlled during the rising phase of tension (Fig. S3). However, it was far more difficult to control tension during relaxation, where inhomogeneity in the sarcomeres can occur, leading to differential shortening and lengthening of sarcomeres. We considered whether a loss of SL control could induce a discrepancy in the recorded tension. Specifically, the Caremani et al. tension transients showed some lengthening at the sarcomere level during the relaxation phase (Fig. S3). This lengthening increased tension in a passive parallel elastic (PE) component that is not present if the half-sarcomere is kept isometric. Consequently, the adjusted tension for the effect of PE (dark gray line in Fig. 3) had faster relaxation and lower resting tension than the observed tension in the Caremani et al. experiment.
The obtained prediction of the twitches using our five-state model shows a good match to the peak tension, but tension rise was slightly slower, and there was a poor overall fit to the relaxation phase (dark green line in Fig. 3), with the fall in tension being too slow, compared with the observed tension transient (dark gray line in Fig. 3). The parameters used in the simulations of the isometric half-sarcomere are listed in Table 1, except parameters and which were adjusted for the isometric conditions having values 50 s−1, 21 s−1, and 150 s−1, respectively. In addition, the five-state model also predicted a high level of resting tension, of the order of ∼40 kPa, which was corrected by lowering baseline [Ca2+] (Fig. 3, pink dashed line) to the observed tension values of <10 kPa. Both the high resting tension and the slow relaxation were the result of a high population of cross-bridges in the M.D.Pi state that can bind to actin to generate tension even at low [Ca2+]. Thus, there are too many bound, tension-generating cross-bridges at baseline [Ca2+] levels and during relaxation after [Ca2+] becomes low.
To resolve this problem, we explored the possibility that changes in a few key rate constants of the five-state model could increase relaxation rate and closely replicate the observed twitch. A reduction in the number of available cross-bridges could potentially be achieved by weakening the association equilibrium constant for cross-bridge binding by increasing the value of An increase of to 350 s−1 brought the rate of relaxation closer to that observed but also began to slow the rate of tension rise (see Fig. S2, green line versus dark gray line). A 350 s−1 decreased from 2.38 to 0.14, which reduced the probability of cross-bridge binding to an available actin site from 70% to 12%. Such a low probability would be incompatible with the tension levels observed at full calcium activation. Similarly, assuming an effective actin concentration of 1 mM in the sarcomere, the predicted affinity of the cross-bridge for actin would be 420 µM for KA of 2.38 and ∼7 mM for a KA of 0.14. The value of [A]/KA = 7 mM is far too high, resulting in low estimates of actin affinity comparing to observed values in ATPase assays.
We further explored the possibility of various combinations of parameters (see Figs. S16, S17, S18, S19, S20, S21, S22, and S23) so that model predictions with the five-state model match the observed tension responses. The best match of MUSICO five-state model predictions with experimental observations was achieved by varying values of and shifting forward power stroke rate for nm (Fig. S22). The shift (see inset in Fig. S22) effectively increased the number of cross-bridges that can complete the power stroke and go through a full ATPase cycle. The latest step increased the cross-bridge detachment and rate of tension relaxation, enabling the model predictions to match the observations (Fig. S22). However, this was only achieved by a reduced level of baseline [Ca2+], as used in Fig. 3 for simulations with the five-state model. Moreover, the simulations with the observed Ca2+ transients showed the high resting tension and slightly higher tension peak, and they slowed the rate of tension relaxation (Fig. S23). We further explored other values of parameters, but we were not able to reduce the resting tension. The reason is shown in Fig. S24, where the tension–pCa relationship, obtained using the parameters that achieved very good fit to the observation (displayed in Fig. S22), shows increased Ca2+ sensitivity at low [Ca2+]. Importantly the tension at baseline [Ca2+] shown in Fig. S24 is approximately the same as the resting tension in Fig. S23 at the same [Ca2+].
An alternative way to reduce the number of available M.D.Pi cross-bridges, accelerate relaxation, and lower resting tension is to incorporate a regulated thick-filament relaxed state (see Fig. 2 B). Such a relaxed state of the thick filament, called the “PS,” has been widely discussed in the literature in recent years, and the background to this is introduced in the Materials and methods section (Linari et al., 2015; McNamara et al., 2015; Irving, 2017).
There are several ways in which a relaxed thick filament could be regulated and implemented in our model, as outlined in the Materials and methods section. The PS could be regulated via strain in the thick filament, via phosphorylation of the RLC and MyBP-C, or Ca2+-dependent association/dissociation from the thick filament backbone. To implement any of these in the spatially explicit 3-D model, we would need to know details of myosin head packing onto the thick filament, the degree of cooperativity in the transition between the parked and active states, and any necessary differences between, for example, the C and D zones of the thick filament, as suggested by Brunello et al. (2020). In the absence of direct experimental evidence for a detailed thick-filament activation mechanism, and for simplicity, we made the PS Ca2+ sensitive with a Ca2+ affinity of 1 µM. This allowed thick-filament regulation to mirror that of the thin filament (see Eq. 1 in the Materials and methods section for a definition of Ca2+-dependent state transition rates). The simulations using the six-state model (red line in Fig. 3) substantially improved matching of the observed experimental resting tension and tension transients, and, importantly, the relaxation time course was close to observations. The parameters used in our simulations of isometric sarcomere experiments that are different from those listed in Table 1 are shown in Table 2.
We also did attempt to fit the transient of Janssen and de Tombe (see Fig. S4). We could generate a reasonable description of tension rise and peak tension; however, the relaxation phase was inconsistent with the observed tension trace by Janssen and de Tombe. This could be consistent with a problem controlling the SL during relaxation similar to that of Caremani et al. The parameters used in these simulations were adjusted for the 5.2°C-lower temperature than that used by Caremani et al., and they are listed in Table 2.
Importantly, simulations with the cross-bridge model with the PS showed significantly reduced resting tension and had little effect on the rate of tension rise or peak tension, but they accelerated the relaxation phase and resulted in a near-perfect fit to the adjusted tension transient (Fig. 3, red line). Thus, this six-state model was able to effectively simulate relaxation, which has been a major challenge in modeling data where tension is not normalized, because detaching myosin heads entering the PS was not available to reattach and contribute to the maintenance of tension.
Tension and SL responses to Ca2+-evoked twitch transients in fixed-end trabeculae
As noted above, sarcomeres are not isometric during an observed “fixed-end” trabecular twitch, as shown by Janssen and de Tombe, 1997; Caremani et al., 2016, who simultaneously measured the tension and half-sarcomere displacement during the twitch. The sarcomere-level shortening and lengthening due to trabecular compliance affects the time course and magnitude of the observed tension (Caremani et al., 2016). The observed tension–displacement loop shows higher tensions during sarcomere shortening and lower tensions during sarcomere lengthening (Fig. 4 A, orange circles). Accordingly, sarcomeres shortened during tension rise and lengthened during relaxation. The next step in development of the model was to incorporate a compliance originating from the series elasticity of cardiac muscle structural components (Fig. 1).
For simplicity, we derived a nonlinear serial spring-like elasticity (SE) from the average loop tension–displacement relationship (gray line in Fig. 4 A). The difference of tension (orange circles) to the mean (gray line) in Fig. 4 A showed an approximately symmetric shape along the path of the loop (Fig. 4 B), reflecting the viscoelastic effects due to shortening and lengthening velocities. A similar relationship was derived from the Janssen and de Tombe experiments (see Fig. S5), where the magnitudes of the peak tensions and displacements show about the same values, but viscoelastic effects were a little different, more so during relaxation.
Caremani et al. did not observed [Ca2+] transients for the fixed-length trabecula at 27.2°C. Therefore, as for the isometric sarcomere twitches above, we used the time course of the [Ca2+] transient from Janssen et al. (2002) at 27.5°C, but with a reduced peak from 1.51 to 1.20 µM to match the lower tension at the lower temperature. The corrected [Ca2+] trace, proportionally scaled down from Janssen et al., is shown as a pink line in Fig. 5 B.
The tension and [Ca2+] transients shown on Fig. 5, A and B, reflect the change in temperature between the two conditions. In all simulations, we used the six-state model with the same parameters as for the isometric sarcomere data in Fig. 3, and the nonlinear trabecular SE stiffness (Fig. 4 A and Fig. S5). The MUSICO simulations did not produce a good match to the observed transients (not shown). To produce a match to the tension transient, three parameters needed to be adjusted: the forward and reverse rates for myosin attachment to actin and the power stroke free energy change, Their values used are listed in Table 2. The simulations, shown in Fig. 5, A and B (green line), describe the tension transients (red line) quite well in both the rising and relaxing tension phases. However, during tension relaxation in both experiments, the predicted half-sarcomere lengthening was slower than experimentally observed.
A principal cause of the inability to fit both tension and displacement could be the oversimplified nonlinear stiffness used for the SE component, which neglects the hysteresis in the observed tension–displacement data caused by the viscoelastic effect shown in Fig. 4 A and Fig. S5. However, the model predictions were able to describe tension and the size of the SL changes in fixed-length trabeculae (Fig. 5). Because measurements of tension are more precise and more reliable than measurements of displacements, in the following analyses and comparisons, we focus on the tension transients.
The adjustments required to simulate the fixed-length trabecula versus isometric sarcomere tension transients were similar at both temperatures used (Table 2). The values of and increased approximately twofold in each case (range, 1.8–2.2-fold), whereas remained unaltered and increased by ∼15%. This adjustment can be accounted for to a significant extent by the shortening/lengthening of sarcomeres that occurs during the twitch. This change in structure affects the number of accessible actin-binding sites for myosin through changing the alignment of the thick and thin filaments and alterations to the strain-dependent rate constants (k+A, k−A). Note that because both and change by the same amount, there is no change in the equilibrium constant of M.D.Pi binding to actin. For a sensitivity analysis of the adjusted parameters, see Figs. S8, S9, S10, S11, S12, and S13.
It is interesting to note that lengthening and shortening also result in significant changes in the cross-bridge cycling rates (ATPase rates; Figs. S6 and S7), as expected from the Fenn effect (Fenn, 1924). Length effects on Tn Ca2+ release/uptake rates have also been reported (Kataoka et al., 2007). Moreover, it is worth noting that correcting for the 5.2°C temperature change between the two datasets was identical for the isometric sarcomere and the fixed-end twitches. This is consistent with the temperature effect being due to the temperature sensitivity of the underlying molecular events (actin binding and ADP release).
To test our model further, we took a third independent fixed-length rat trabecular twitch dataset from Janssen et al. (2002) collected at three different temperatures (see Fig. 6). These data were collected under conditions similar to the earlier measurements and indicate a slight increase in peak tension, as well as marked acceleration of both tension rise and relaxation with a temperature increase from 22.5°C to 30°C (Fig. 6, red lines). Note that these transients are not identical to those of Fig. 5 because the [Ca2+] transient has a slightly higher peak and hence also a higher peak tension. This provides a reasonable test of the precision of our model. For simulations of tension transients at the three temperatures (Fig. 6, green lines), only three parameters needed to be adjusted (Table 2). Our simulations used the identical set of parameters, including the serial elastic element parameters, used for the 22°C (Janssen and de Tombe, 1997) and 27.2°C (Caremani et al., 2016) datasets (Table 1 and Table 2). The only exception is a slightly increased rate of from 140 to 149 s−1 for an increase in temperature of 0.5°C (from 22 to 22.5°C). The simulations at temperatures paralleling the previous two experiments (Fig. 6, green lines) show an excellent match to tension traces for the Janssen et al. (2002) experiments (Fig. 6 A, red lines). Also, the predicted tension–[Ca2+] loops (solid lines) matched the observations (circle symbols) at all three temperatures (Fig. 6 B). Notably, further increasing temperature to 30°C also produced a modest increase in peak tension and in the rates of tension rise and relaxation. The simulations required a further increase of parameters and by factors of 1.13, 1.4, and 1.5, respectively, for the 2.5°C temperature increase (from 27.5 to 30°C). This is in line with those required for the ∼5°C temperature change from 22 to 27.2°C. Thus, this independent dataset with higher Ca2+ transients and higher peak tensions required no special adjustments from experiments used in Figs. 3 and 5.
Predictions of tension–pCa relationship in intact muscles with parameters from twitch contraction simulations
A major challenge in simulating cardiac muscle contraction is to predict the [Ca2+]50 and Hill coefficient in the force–pCa relationship of intact muscle. We took the parameters from Table 1, and the cross-bridge kinetics at 22.5°C listed in Table 2, used to generate fits of the twitch data of Janssen et al. (2002) at 22.5°C (Fig. 6 A), to simulate values of the force–pCa relationship (Fig. 7) for fixed-length trabeculae. The simulations with the six-state model (Fig. 7, red line) matched very well the force–pCa data of Janssen et al. (2002) (Fig. 7, black filled squares), including the pCa50 and Hill coefficient (nH) values. The simulations also fit the intact muscle tension–pCa data from Gao et al. (1994) (Fig. 7, gray triangles), which were almost identical to those of Janssen et al. (2002). It is notable that the simulations, with no change in values from those used in Figs. 5 and 6, generated not only the correct pCa50 but also the observed high value of the Hill coefficient. Models that rely solely on Ca2+ regulation of the thin filament via CFC (i.e., intrinsic cooperativity) do not typically produce values of the Hill coefficient above 3 in cardiac muscle (Mijailovich et al., 2019). The higher Hill coefficient, shown in Fig. 7, is a consequence of the cooperativity assumed in the regulation of the thick filament, in addition to thin-filament mechanisms. It is important to note that other models based on cooperativity factors can achieve higher cooperativity factors (Campbell et al., 2001; Rice et al., 2008), but these are parametric models, and these kind of models cannot be compared with structurally based models such as CFC (Geeves et al., 2011; Mijailovich et al., 2012; Mijailovich et al., 2019).
The MUSICO platform, which includes explicit 3-D sarcomere geometry, several multistate cross-bridge cycles, and Ca2+ regulation of thin-filament activation by a CFC (Mijailovich et al., 2016; Mijailovich et al., 2019), has now been further developed to allow simulations of twitch contraction in cardiac trabeculae. The key new implemented features reported here include a cross-bridge cycle with a PS that enables regulation of thick-filament activation and an SE component observed in twitch experiments in fixed-end trabeculae. In this study, we showed that both features are necessary for our model to simulate twitches and predict the tension–pCa relationship from twitches.
Our simulations show that we can describe three independent datasets for fixed-end trabecula twitches (Janssen and de Tombe, 1997; Janssen et al., 2002; Caremani et al., 2016) using the same parameters after corrections for the different temperatures. The corrections required to the three parameters and listed in Table 2 were modest and in line with expectations of the temperature dependence of protein reactions. Similarly, we have simulated the Caremani et al. (2016) data in which SL was held constant during the rising phase of tension at 27.2°C and the corrected tension during twitch relaxation phase for the distortions of the PE elements caused by imperfect sarcomere control during relaxation. The corrections between the isometric sarcomere and the fixed-end twitches were twofold decreases in the rates of attachment and detachment of myosin to actin, which may reflect changes in accessibility to actin sites during lengthening/shortening.
In addition to the sensitivity analysis in the supplemental material, two other features came out of the model that were not part of the initial objective. The same parameters used to describe the twitches in Figs. 5 and 6 were able to generate a force–pCa relationship observed in two previous studies (Gao et al., 1994; Janssen et al., 2002). Similarly, the same sets of parameters were able to predict the change in the ATPase/cross-bridge cycling rate during the rising and falling phases of the tension transient. The cycling rates, as expected, change due to the shortening and lengthening of the sarcomeres (Fenn, 1924).
Why a PS may be necessary for cardiac muscle and normal ventricular function
Simulations using a minimal five-state cross-bridge model, coupled with thin-filament activation by Ca2+ in intact muscle, presented three apparently unrelated problems: slow tension relaxation during the twitch (Fig. 3), high levels of resting tension, and unrealistic values for the myosin (cross-bridge) detachment rate, to bring the tension relaxation rate toward the observations (Fig. S2). High resting tension was also predicted in simulations of demembranated cardiac muscle fibers and at low Ca2+ (Prodanovic et al., 2020). These observations suggest a missing component in the underlying mechanisms of activation and relaxation. This conclusion is supported by recent studies indicating that, in addition to thin-filament activation, thick-filament activation may be an important process (Linari et al., 2015; McNamara et al., 2015; Irving, 2017; Reconditi et al., 2017). This is especially true when considering length-dependent activation (Ait-Mou et al., 2016; Caremani et al., 2019) and the contractile abnormalities that occur in cardiac myopathies associated with mutations in thick-filament proteins (Hooijman et al., 2011; Spudich, 2015; Nag et al., 2017; Anderson et al., 2018). The thick filament may contain one or more “off” states where myosins are not available to bind to the thin filament. For example, both structural (Zoghbi et al., 2008) and biochemical energetic studies suggest an ordered state, the so-called SRX state proposed by Cooke (McNamara et al., 2015), and a disordered state still associated with the thick-filament backbone (Anderson et al., 2018; Brunello et al., 2020) that has extremely low probability for myosin heads to bind actin. One possibility for the transition from a myosin PS to an “on” (i.e., M.D.Pi) state could be triggered by thick-filament strain (Linari et al., 2015; Irving, 2017). This kind of mechanism has been implemented in a model of cardiac contraction where the rate of transition from the myosin PS state to M.D.Pi state is proportional to force (Campbell et al., 2018). This model provides reasonable predictions of length-dependent activation and twitch contractions in cardiac muscle, but the model used simplified cross-bridge cycle kinetics and sarcomere geometry and did not take into account the observed significant SL changes (Janssen and de Tombe, 1997; Ferrantini et al., 2014; Caremani et al., 2016). In addition, these good fits are achieved by an empirically defined flux of thin-filament regulatory units to an on state tightly coupled to the Ca2+ transient and inclusion of a cooperativity factor that amplifies the number of on states during calcium rise and attenuating during the fall of [Ca2+] levels. The cooperativity factor was necessary to close units during tension relaxation because the thick filament is strongly activated at high tension, and relaxation was prolonged without this factor. The concept of thick-filament activation rates coupled to force should be further investigated to provide possible mechanisms. Moreover, a recent paper (Ma et al., 2018a) suggested that thick filament–based regulation is likely to be more complex than the previous studies implied and that other mechanisms of thick-filament activation could be operating during twitch contractions, as suggested, for example, in our preliminary studies (Prodanovic et al., 2020) and justified here.
In our model, the PS acts as a Ca2+-dependent buffer for myosin cross-bridges, increasing the available myosin heads on Ca2+ activation and reducing the population of M.D.Pi state during relaxation. This was achieved in the model by the transition rate out of the PS, having a slow rate at low [Ca2+] and a fast rate at high [Ca2+] concentrations (Eq. 1). Such a model may be supported by x-ray diffraction observations that myosin heads rapidly move radially away from the thick-filament backbone during fast Ca2+ activation of frog skeletal muscle (Huxley et al., 1994), much faster than the increase in the number of attached cross-bridges and subsequent development of tension. These unbound myosin heads in the proximity of actin filaments increase the pool of cross-bridges available to bind to the thin filament. A Ca2+-dependent transition rate from PS is one of multiple alternative possibilities for activation of the thick filament that needs to be further investigated experimentally before it is fully confirmed.
Our implementation of the PS illustrates how the presence of this PS can slow the rate of tension rise, reduce resting tension, and accelerate relaxation by limiting the number of myosin heads available to interact with actin. It follows, therefore, that if the occupancy of the PS can be manipulated by phosphorylation events (e.g., RLC, MyBP-C), small molecules (e.g., mavacamten, deoxyadenosine ADP), or mutations in myosin or MyBP-C, then alterations in all three twitch parameters can be expected. Indeed, several recent studies have suggested that many hypertrophic cardiomyopathy–linked mutations in myosin or MyBP-C cause hypercontraction by reducing the occupancy of the SRX state, allowing more myosin heads to participate in a contraction, and that this can be reversed by molecules such as mavacamten. The SRX and DRX states observed in experiments is equivalent to our PS in the model presented here. Our model may therefore be able to predict the results of mutations and small molecules on the cardiac myocyte twitch. Effects on twitch relaxation times similar to those anticipated in the absence of PS here were reported by Toepfer et al. (2020) for cardiac induced pluripotent stem cells expressing hypertrophic cardiomyopathy–linked mutations in myosin that were reversed by mavacamten.
Series elasticity of trabeculae
Fixed-length contractions of cardiac muscle can result in relatively large changes in SLs during twitch contractions (∼7%; Janssen and de Tombe, 1997; Caremani et al., 2016), indicating significant compliance in intact muscle preparations. Current simulations of the cardiac twitch contraction typically use models of the half-sarcomere (though, on rare occasions, multiple sarcomeres in series) where the thick- and thin-filament compliances account for only a small fraction of the observed compliance in the cardiomyocytes or tissue. The first attempt to model the SE component was reported by Parmley and Sonnenblick (1967), who proposed schemes of CE in series and in parallel with the elastic elements (SE and PE). This kind of model has been used in simulations of twitch contractions in trabeculae; however, prediction of a highly simplified model only showed the effect of the variation of SE and PE elements but not in milieu or experimentally observed data (Rice et al., 2008). In a more detailed model (Campbell et al., 2018), SE is apparently incorporated, but the changes in SLs reported are small and largely underestimate the observed values. It is likely that they considered only compliance of thick and thin filaments and not compliance of trabeculae. Our results presented here and elsewhere (Prodanovic et al., 2020) suggest that analysis of subtle differences in cardiac muscle dynamic contractions associated with, for example, mutations in sarcomere proteins requires a thorough description of the series elasticity observed in cardiac muscle in order to be most accurate and informative.
A simplified definition of the series elasticity, derived from the observations of Caremani et al. (2016), was used in our simulations and provided reasonable but not exact fits of tension transients and sarcomere displacements during twitch contractions of fixed-end trabeculae (Fig. 5). Specifically, it was difficult to achieve good fits for both tension and displacements, particularly during relaxation. The reasons for this may come from both experimental and model limitations. On the experimental side, measurements of SLs are obtained by different techniques, often with different outcomes. For example, Caremani et al. (2016) used a striation follower to monitor SL in the middle portion of the trabeculae, whereas Ferrantini et al. (2014) used laser diffraction. In the first case, the change in a half-sarcomere length, denoted as displacement, represents an average of half-sarcomere length changes in the middle portion of trabeculae, where the SLs are more uniform. In the laser diffraction case, the half-sarcomere length is determined as an average of all sarcomeres reflected from the laser beam. The use of different measurement techniques resulted in large variations in recorded displacement traces, such as between observations of Caremani et al. (2016) and Ferrantini et al. (2014). In addition, either of these techniques cannot fully control the inhomogeneity within these SLs during tension relaxation, particularly in the chaotic phase of relaxation (Brunello et al., 2009).
The use of simplified elastic characteristics in the model instead of the observed hysteretic loop (Fig. 3) is in part responsible for imperfect matches of displacements in simultaneous traces compared with the observations of sarcomere lengthening during tension relaxation (Fig. 5, A and B). The exact tension–length relationship (orange circles in Fig. 4 A) can be implemented in MUSICO, but we prefer the simplified model (thick gray line in Fig. 4 A). This is because the observed force loop depends on the magnitude and rate of developed tension and needs to be recorded separately for each experiment, information that is usually not available. Thus, we believe that fitting the tension traces is a more general and more robust approach than fitting displacements.
The effect of changes in calcium sensitivity on twitch contractions
The majority of studies on cardiac muscle contraction have focused on force generation and the Ca2+ sensitivity of steady-state force that can be relatively easily assessed as a function of [Ca2+] in the activation solutions. Accordingly, an abundance of experimental data in the form of force–pCa relationships provides evidence for functional changes in cardiac muscle associated with cardiomyopathies and their responses to treatment by prospective drugs. However, these data cannot be directly translated to physiological transient contractions in the living heart. Our simulations predicted the tension–pCa data of Janssen et al. (2002) and Gao et al. (1994) (Fig. 7) with the same model parameters as for the twitches at the same temperature. This is a promising step toward the possibility of interrelating the Ca2+ sensitivity from steady-state measurements to tension responses to Ca2+ transients. We believe that the reverse process is possible: to predict twitches with parameters estimated from tension–pCa relationships of cardiac muscles under similar experimental conditions. However, the problem is that most force–pCa data originate from the measurements in demembranated cardiac muscle.
The experimentally observed difference in the Ca2+ sensitivity of tension between intact and demembranated cardiac muscle is well recognized, but there are only a few studies that have quantitatively assessed the magnitude of these differences (Gao et al., 1994). The [Ca2+] used to determine the force–pCa relationship observed in demembranated muscle is well controlled, but interfilament spacing and experimental conditions such as SL, pH, phosphate, and ionic strength can be significantly different from those in intact muscle cells. To the best of our knowledge, only Gao et al. (1994) have measured the force–pCa relationship for both conditions in the same muscle (Fig. 7). The sensitivity in demembranated muscles was significantly lower than in intact muscles, and it is strongly dependent on experimental conditions (Fig. 7). The differences in experimental conditions between these preparations affecting the Ca2+ sensitivity, including pCa50 and cooperativity, nH, are shown in Table 3. Bringing the experimental conditions in demembranated muscles close to those observed in intact muscle reduces the differences in both pCa50 and nH between the preparations. The strongest effect is reduction of free Mg2+ to the physiological value of 0.5 mM that reduces the difference between intact and demembranated muscle to only ∼0.2 pCa units. Interestingly, this difference in sensitivity can be attributed to the effect of the large increase in interfilament spacing, from d10 = 33.83 nm in intact muscle to 41.66 nm in demembranated muscle (both at SL = 2.2 µm). This increase in sensitivity is similar to that observed experimentally when the lattice of demembranated muscles is compressed back to intact muscle values with dextran (Cazorla et al., 2001; Martyn et al., 2004). Our simulations of these experiments predicted a change of 0.15–0.25 pCa units (Mijailovich et al., 2019) that can bring the sensitivity values of demembranated muscle at 0.5 mM free Mg2+ to the values observed in intact muscle. This is a promising result, but it is based on a single experiment, and to fully develop methodology to translate the force–pCa data from demembranated muscle to twitches, it is necessary to confirm it by stronger experimental evidence. This is additional evidence of how strongly coupled experimental data and computational simulations could translate valuable data to physiologically relevant cases.
Active tension and role of titin in length-dependent calcium sensitivity
Tension in cardiac muscle is measured as the aggregate of active and passive tension. The active tension reflects forces generated by bound cross-bridges and the passive tension as the collective contribution of titin and collagen (Cazorla et al., 2001). The passive tension is ∼0 at slack length, but at longer SLs the tension progressively increases nonlinearly (Mijailovich et al., 2019). During twitch contractions of the isometric trabeculae, a significant shortening of sarcomeres, up to 80 nm per half-sarcomere, can be observed and is due to series elasticity (Fig. 6). The shortening during a twitch starting from SL = 2.2 µm decreases force in titin from ∼16.5 pN per myosin filament to ∼0 at SL = 1.9 µm; thus, the active tension has slightly higher peak magnitude than reported in the literature, where initial (passive) tension is simply subtracted from total tension measured (Fig. S14).
Potential model limitations: variation of the model parameters, uncertainties of experimental evidence, and parameter estimation methodology
Fitting multiple sets of data presents numerous challenges. The model is imperfect (as are all models), and multiple sets of parameters can fit the desired set of data. On the other hand, there is a need to reconcile the large variability in the data obtained in the same kind of experiments from different laboratories with slight or no obvious changes in experimental conditions. There are also differences in preparations (e.g., intact versus demembranated trabeculae) and in modes of contraction, such as transient twitch responses and steady-state contractions usually displayed as force–pCa relationships. Here, we minimized the number of fitting parameters by using experimentally determined values whenever possible. Examples include SL and lattice spacing and other well-established experiments such as ATPase measurements and motility assays. The remaining parameters are estimated as best guesses based on experience and experimentation. However, some parameters can inherit hidden deviations due to specific conditions of the experiments that could be quite different in intact or demembranated muscles. For example, demembranated muscles show a large increase in sarcomere lattice spacing that is associated with decreases in myosin-binding rate and Ca2+ sensitivity compared with intact muscle. Experiments often have their own large variation, with differences in magnitude and shape of tension responses during twitches that may reflect seemingly minor differences in experimental conditions. A second large hurdle is the incompleteness of most datasets, where simultaneously measurements of tension and changes in SL during twitch contractions of trabeculae with fixed lengths and Ca2+ transients are not available. Thus, with these multiple uncertainties on the experimental side, it is exceedingly difficult to fit all experiments with the same set of input parameters. The strength of the detailed modeling approach presented here is the ability to distinguish the contributions of each effect individually to the overall recorded responses where multiple effects operate in aggregate. For example, the contribution of changes in interfilament spacing between demembranated and intact trabeculae to the Ca2+ sensitivity of tension development can be estimated from MUSICO simulations (Mijailovich et al., 2019). The estimate of such an increase in sensitivity of demembranated muscle due to an increase in lattice spacing is only a small fraction of the total increase in Ca2+ sensitivity. The remaining differences need to be attributed to other mechanisms when they are known.
Developing an accurate model of a twitch contraction is a serious challenge that requires multiple theoretical concepts to be integrated into simulations to account for diverse experiments. Our current model includes the Ca2+ transient, mechanisms of cooperativity between thick- and thin-filament activation processes, series elasticity of cardiac tissue, and the force-generating cross-bridge cycle. This was sufficient to provide an accurate description of a twitch contraction of rat cardiac muscle, enabling simulations of tension and SL changes in fixed-length trabeculae as a response to intracellular Ca2+ transients. To accomplish this, two new features needed to be introduced: an SE element and a process for thick-filament myosin activation and deactivation. An elastic component that is adjustable during tension rise and fall is essential to account for the internal sarcomere shortening that occurs in fixed-end “isometric” cardiac muscle contractions. Using thin-filament regulation alone was not sufficient to match the observed relaxation patterns. The main reason for simulations having slow relaxation was the large population of myosin heads in the M.D.Pi state. These heads can rebind to actin filament and therefore impede the relaxation process. Inclusion of a PS (either DRX or SRX state), where myosin heads interact with the thick-filament backbone, provided a mechanism to remove heads from the cycling cross-bridge pool. Our simulations also account for the effect of change in sarcomere interfilament spacing and differences in local experimental conditions (e.g., ionic strength, SL, and temperature). Thus, our simulations demonstrate the need for spatially explicit models to understand the multiple interacting processes that occur within cardiac muscle during twitch contractions. Furthermore, in future studies, the MUSICO platform developed here may provide a powerful tool to explore the effects of point mutations in sarcomeric proteins on twitch transient responses, the effect of the degree of penetrance of the mutated protein(s), and the potential influence of small molecules on rebalancing altered twitch kinetics.
Strain-dependent rates in the actomyosin cycle
For ATP binding and cross-bridge detachment rates (i.e., transition between A.M. and M.T states), we assume for now that constant, although it depends on ATP concentration, and due to a large increase in
The state transition rates between detached states M.T and M.D.Pi, and are independent of strain. The forward rate, involves a conformational change in the head following the opening of the 50-kD cleft leading to the recovery stroke and hydrolysis of ATP transitioning into the M.D.Pi state.
The transition of the M.D.Pi state into the PS is defined by the rate constant that could potentially depend on [Ca2+], myosin light chain, and MyBP-C phosphorylation. For simplicity, we first consider only the effect of [Ca2+] by setting, for example, constant and in the form of a Hill curve described by Eq. 1.
Thin-filament regulation by calcium
CFC angular position and its variance
The angular displacements of the Tpm–Tn chain are defined by the function of thermally excited chain configurations defined by the Feynman path integral (Feynman and Hibbs, 1965; Smith, 2001), where the current configuration of bound TnIs and myosins at defined positions, along F-actin, the path of minimum energy represents the mean chain angle and the SD of chain angles arising from thermal excitation (Smith, 2001; Smith et al., 2003; Geeves et al., 2011; Mijailovich et al., 2012).
If there are no constraints, the mean configuration of the CFC is in the closed state sitting on the bottom of the confined potential at the angle and the SD of fluctuating angle, has a value of about where R is the radius at which Tpm sits on the actin filament, and is the confined persistence length of the CFC (Smith et al., 2003). When the CFC is constrained, the chain is pinned at the positions where TnIs are bound to actin at an angle forming a local blocked state, whereas the chain is assumed to be conditionally constrained at positions where myosins are strongly bound to actin, preventing the CFC from reaching angles forming a local open state. At positions where myosins are weakly bound, the chain is conditionally constrained, preventing the CFC from reaching angles The pinning and conditional constraints guarantee that the blocked, open, and closed states correspond to the three orientations seen in cryo-EM images defined by angular displacements and (Pirani et al., 2005; Poole et al., 2006).
Chain-regulated cross-bridge kinetics
The coupling between the Ca2+-regulated position of the CFC and the state myosin cycle involves (1) modulation of myosin binding and the transition between weakly and strongly bound myosin states by the current position of the CFC to the actin site and (2) restriction of TnI rebinding to actin by nearby bound myosins. The biochemical TnI-A states include a bound state of TnI to actin that maintains Tpm–Tn in the position inhibiting myosin binding and in the other state, where TnI is not bound to actin, allowing the Tpm–Tn chain to move azimuthally along the actin surface. The CFC chain on the actin surface is not static but dynamically moves azimuthally, except at locations where TnI is bound to actin, and permits myosin binding for a fraction of time when a binding site on actin is available. Because these azimuthal fluctuations are much faster than myosin binding, the fraction of time that an actin site is available for weak myosin binding is proportional to the probability that the local position of CFC is at and the transition from weak to strong binding, when These probabilities are calculated from the local mean angular position of the chain, and its SD, where i denotes the actin site at the discrete position, along an actin filament strand (Mijailovich et al., 2012; Mijailovich et al., 2019). Conversely, bound myosins reduce the mobility of the CFC and modulate TnI rebinding to actin. In this case, the fraction of time that a TnI can reach its actin site is proportional to the probability that the local position of CFC is at
The coupling between cross-bridge cycling and thin-filament regulation by calcium
Modulation of TnI–actin transitions by bound myosins
In the absence of Ca2+, most TnIs are bound to actin, but in the presence of Ca2+, calcium binds TnC in a concentration-dependent manner. For simplicity, we assume that has a fixed value, whereas the binding rate of TnI to actin is a weighted function of the rate of TnI binding to actin from the unweighted closed (or open) state rate, (Mijailovich et al., 2012). The transition probability of attachment of TnI to actin is where is a weight factor and the transition probability for the detachment of TnI from actin is simply The first Monte Carlo drawing is performed over all TnCs along all CFCs, where each actin filament has two CFC strands. If any change of TnI state is drawn, the TnI state is updated for calculation of the CFC configuration in the current time step.
Modulation of myosin–actin transitions by the CFC
The probabilities of changing state in the five-state (and six-state) cross-bridge model are constructed so that each state can transition to two or three neighboring states. The transition from the detached state M.D.Pi, associated with axial strain-dependent rates, includes attachment probability, adjusted for binding to multiple binding sites on actin. For the cross-bridges in the detached state, the attachment probability is shared between all reachable actin states defined as where and Because is the sum of the probabilities of attaching myosin heads to each of reachable sites to the equivalent axial strain–dependent binding rate is set to provide the same flux as the probabilistic binding rate The scaling factor, decreases the magnitude of at each by an average number of reachable binding sites to myosin heads (Mijailovich et al., 2016). The weight factors and are associated with the azimuthal position of actin filaments in the sarcomere lattice relative to the myosin head, angle and azimuthal angle of actin site (see Fig. S1 and Mijailovich et al., 2016). The weight factors are proportional to the fraction of time step when these transitions are possible (i.e., to the weighted probabilities that the CFC is at positions that modulates weak myosin binding, denoted as subscript (Mijailovich et al., 2012).
The transitions from the weakly to the strongly attached states (2 and 3) are defined by the probabilities for the prestroke state 2 (A·M·D·Pi) and for the poststroke state 3 (A·M·D). Calculation of these probabilities is almost identical to that reported by Mijailovich et al. (2016), except that includes the weight factor that modulates weak to strong myosin binding, denoted as subscript The factor is proportional to the fraction of time step when the isomerization is possible (i.e., Mijailovich et al., 2012).
The 3-D sarcomere lattice in vertebrate striated muscle is composed of myosin filaments, each surrounded with six actin filaments and each actin filament with three myosin filaments. Each half-length myosin filament is associated with six titin molecules (Mijailovich et al., 2019). Actin filaments can vary in length, depending on species, muscle type, and presence or absence of nebulin in skeletal muscle. The actin filaments have a monomer spacing of 2.7359 nm, and the half period of one strand is 35.56 nm under relaxed conditions (Huxley et al., 1994; Wakabayashi et al., 1994; Prodanovic et al., 2016). The length of a myosin filament is ∼1.58 µm, having 50 crowns (i.e., 150 myosin molecules per half-thick filament), with a crown spacing of 14.3 nm (Luther et al., 2008). This number of crowns and amount of spacing provides maximum overlap with a thin filament of ∼0.7 µm. The actin radius is and the myosin radius is (Mijailovich et al., 2016). The lattice interfilament spacing depends on SL, the specificities of the muscle type, and experimental conditions. The values of the spacing d10 are specified for each set of simulations.
Interfilament lattice spacing
The interfilament spacings in intact trabeculae are much smaller than in demembranated trabeculae having values of ∼33.83 nm at SL 2.2 µm. All intact trabecula fixed-length simulations used this value. In contrast, the interfilament spacings of demembranated cardiac trabecula muscle were obtained at SL 2.25 µm, having values of ∼41.36 nm (Irving et al., 2000; Cazorla et al., 2001) and used in simulations of force–pCa relationships observed by Kreutziger et al. (2011). The interfilament spacing is an important factor in modulation of myosin binding to actin and possibly other rates in the cross-bridge cycle (Williams et al., 2013; Mijailovich et al., 2019).
Actin and myosin filaments are extensible with filament moduli (elastic modulus times cross-sectional area) derived from x-ray diffraction or direct measurement: for actin, for myosin, (Huxley et al., 1994; Kojima et al., 1994). We used the modulus, rather than stiffness, because the reported stiffness values depend on the filament length, and cross-sectional area is not well defined for myosin and actin filaments. The only nonlinear element included in calculations is titin, whose nonlinear elastic characteristics were described previously (Mijailovich et al., 2019). In filaments, the deformation in contracting muscle is small: <1% (Huxley et al., 1994; Wakabayashi et al., 1994; Ma et al., 2018b). The elasticity of cross-bridges originates from multiple sources and can be nonlinear (Kaya and Higuchi, 2010) but, for simplicity, here we used a linear elastic model, as is usually done in sliding filament models.
Cross-bridge model parameters
To achieve the observed rate of twitch relaxation, we used a six-state cycle that includes a PS (Fig. 2, A and B). Following the approach of Duke (Duke, 1999; Mijailovich et al., 2016; Mijailovich et al., 2019), the state transition rate constants are defined by specifically designed experiments, estimated from indirect observations, or deduced from multiple sources. The simulations used a large number of the same parameters, except a few, that can strongly vary with the experimental conditions. For clarity, we show the complete set of cross-bridge model parameters, where variable ones are taken from simulations of trabecula twitch experiments at 27.5°C (see Table 1). When different parameter values are used for simulations of the experiments performed under different conditions, such as at what temperature, SL, and interfilament lattice spacings and whether the muscle is demembranated or intact, they are shown in the main text or in figure legends.
The cross-bridge cycle parameters for myosin binding to and detachment from actin (Eq. A2) are and for the power stroke, the equilibrium constant (Eq. A3) is defined by the power stroke for ADP release, and the second power stroke Because of the exponential forms in the expressions for the state transition rates, Eq. A3 and Eq. A5 can become very large and can generate numerical problems; they are capped to and When the cap value is reached, the reverse rates are changed to decay exponentially to satisfy the equilibrium constant, The remaining rates of the six-state cycle (Fig. 2) are defined as follows: ATP-binding and myosin detach rate, hydrolysis forward rate and backward rate the transition rate to PS and the calcium-dependent transition rate from PS (Eq. 1), where baseline rate is the rate at high [Ca2+] Hill coefficient of the rate sigmoidal [Ca2+] rise and calcium concentration, when is equal to is In all simulations, cross-bridge stiffness is taken to be and the value for at 27.5°C or calculated at the temperature of the experiment.
CFC model parameters
For calcium-binding kinetics to TnC and interaction of TnI with actin (Fig. 2 C), the equilibrium rate constant of calcium binding to TnC along TnC closed states, is effectively defined via forward constant where and calcium-independent dissociation constant TnI–actin interaction equilibrium rate constant and TnI attachment rate and cooperativity coefficient
For the CFC model, we used the same parameters reported previously (Mijailovich et al., 2012): a Tpm pinning angle, myosin-imposed Tpm angular displacement, the persistence length of Tpm–Tn confined chain and angular SD (Pirani et al., 2005). The strength of the chain confining potential is weaker than used for skeletal muscle (Mijailovich et al., 2012), having a value for of 0.149 pN.
Calcium transients: measurements, conversions, and derivations when they are not available
Used in twitch simulations are [Ca2+] transients reported by Janssen et al. (2002) (shown in Fig. 6) or converted from 340-nm/380-nm fluorescence ratio data reported by Janssen and de Tombe, 1997 using the calibration equation (Backx and Ter Keurs, 1993; Gao et al., 1994; Backx et al., 1995), where and (Fig. 5 A and Fig. S4), or deduced from these measurements when calcium transients are not reported (Figs. 3 and 5 B). In the latter case, we adapted the Janssen et al. (2002) Ca2+ transient recorded at 27.5°C in the fixed-length trabecula to the Caremani et al. (2016) experiment at similar temperature (27.2°C) by reducing the peak concentration from 1.51 to 1.20 µM to match lower observed maximum twitch tension (lower by 27.3%) in the Caremani et al. (2016) experiment.
Because the calcium transients typically achieve a higher peak [Ca2+] in fixed-length muscles than in isometric sarcomere twitches (Janssen and de Tombe, 1997), we rescaled the adapted calcium transient for the Caremani et al. fixed-length trabecula experiment (from the peak of 1.2 µM to the isometric half-sarcomere using the ratio of the peaks of [Ca2+] recorded by Janssen and de Tombe at 22°C).
Conversion between observed tension and isometric force per myosin filament
Due to the stochastic process of myosin interacting with actin, the forces in the myofilaments fluctuate in time and each filament experiences somewhat different force. For comparison with observed isometric tensions, we include in all plots both the tension in kilopascals and the average force per myosin filament, in piconewtons. The scales are related by a factor that takes into account how many myosin filaments there are per unit of the fiber cross-sectional area (Linari et al., 1998; Mijailovich et al., 2016). Because the total number of thick filaments in a myofibril or muscle fiber does not change with experimental conditions, for the estimation of the number of thick filaments per unit of cross-sectional area, we used the lattice spacing, d1,0, at slack length (i.e., at the same length where the muscle cross-sectional area used in the tension calculations is measured).
The lattice spacing in demembranated rat cardiac trabeculae cells at slack SL (1.9 µm) is d1.0 = 43.5 nm (Irving et al., 2000) and the number of thick filaments per μm2 myofibril is ∼458. Taking into account the fraction of the cross-sectional area in cardiac cells occupied by myofibrils as 68% (Cazorla et al., 2000), this reduces the number of myosin filaments over the muscle cell cross-section to ∼311 per μm2. In this case, the scaling factor provides that corresponds to per myosin filament or 64.25 pN per actin filament.
In intact trabeculae, the lattice spacing is smaller and at 2.2 µm, where the cross-sectional area was measured, is d1,0 = 33.83 nm (Irving et al., 2000), and the number of thick filaments per μm2 myofibril is ∼757. Taking into account the fraction of the cross-sectional area in rat trabeculae occupied by myofibrils to be ∼61%, as observed in previous studies (Schaper et al., 1985; Barth et al., 1992), reduces the number of myosin filaments over the trabecula cross-section to ∼462 per μm2. In this case, the scaling factor provides that kPa corresponds to pN per myosin filament or 43.35 pN per actin filament.
Henk L. Granzier served as editor.
We gratefully acknowledge the help of Prof. Thomas C. Irving with editing of the final version of the manuscript.
This project is supported by National Institutes of Health grants R01 HL128368 and RM1 GM131981 (to M. Regnier) and by the European Union’s Horizon 2020 research and innovation program under grant agreement 777204 (to S.M. Mijailovich, M. Regnier, M.A. Geeves, and C. Poggesi). This article reflects only the authors’ view. The European Commission is not responsible for any use that may be made of the information the article contains.
The authors declare no competing financial interests.
Author contributions: S.M. Mijailovich conceptualized the study, introduced new additions into MUSICO platform, including formulation of thin filaments calcium regulation, implementation of updated cross-bridge cycle, series elastic element in the model of trabecula contraction, new structural elements and changes in sarcomere geometry, designed numerical procedures, analyzed the data from simulations, wrote the first draft of the manuscript, and supervised each step of the work. M. Prodanovic assisted in designing numerical procedures and monitored implementation of changes computer code, performed all simulations and data analyses, contributed in graphical design, and assisted in the writing of the manuscript. C. Poggesi provided interpretation of experimental data and contributed in the comparative analysis of simulations and experimental data. M.A. Geeves contributed to development of the model of cross-bridge cycle and thin filament regulation, to formulating the key points in the Discussion, and in writing the final draft of the manuscript. M. Regnier contributed by providing experimental data, formulating the elucidating important points in the Discussion, and assisted in writing the final draft of the manuscript.
This work is part of a special collection on myofilament function and disease.