Various mutations in the structural proteins nebulin and titin that are present in human disease are known to affect the contractility of striated muscle. Loss of nebulin is associated with reduced actin filament length and impairment of myosin binding to actin, whereas titin is thought to regulate muscle passive elasticity and is likely involved in length-dependent activation. Here, we sought to assess the modulation of muscle function by these sarcomeric proteins by using the computational platform muscle simulation code (MUSICO) to quantitatively separate the effects of structural changes, kinetics of cross-bridge cycling, and calcium sensitivity of the thin filaments. The simulations show that variation in thin filament length cannot by itself account for experimental observations of the contractility in nebulin-deficient muscle, but instead must be accompanied by a decreased myosin binding rate. Additionally, to match the observed calcium sensitivity, the rate of TnI detachment from actin needed to be increased. Simulations for cardiac muscle provided quantitative estimates of the effects of different titin-based passive elasticities on muscle force and activation in response to changes in sarcomere length and interfilament lattice spacing. Predicted force–pCa relations showed a decrease in both active tension and sensitivity to calcium with a decrease in passive tension and sarcomere length. We conclude that this behavior is caused by partial redistribution of the muscle load between active muscle force and titin-dependent passive force, and also by redistribution of stretch along the thin filament, which together modulate the release of TnI from actin. These data help advance understanding of how nebulin and titin mutations affect muscle function.
During muscle contraction, myosin binding to actin is primarily regulated by Ca2+, via the thin filament–associated proteins tropomyosin (Tm) and troponin (Tn). There are many models of thin filament regulation, but the kinetics of interactions between the Tm–Tn complexes are, at present, best described by the three-state model of McKillop and Geeves (1993) and the long-range cooperative continuous flexible chain (CFC) model (Geeves et al., 2011; Mijailovich et al., 2012a) that includes structural constraints between Tm–Tn regulatory units (Fig. 1). The interaction of Tn with actin is regulated by Ca2+ binding to troponin C (TnC). TnC has two binding sites in skeletal muscle and one in cardiac muscle (Potter and Gergely, 1975; van Eerd and Takahashi, 1975; Leavis and Kraft, 1978; Leavis et al., 1978; Holroyde et al., 1980). Coupling these regulatory models with cross-bridge models in the explicit 3-D sarcomere lattice (Daniel et al., 1998; Chase et al., 2004; Mijailovich et al., 2016) can provide a general tool for studying the effects of mutations in regulatory and auxiliary structural proteins on muscle function.
The structural proteins nebulin and titin (Fig. 1) are essential for proper physiological function of muscle. The degree of expression and the precise roles of these proteins vary between different types of muscle. Genetic mutations in these proteins are associated with many diseases (Pelin et al., 1999; Gerull et al., 2002; Wallgren-Pettersson et al., 2004; Lawlor et al., 2011; Taylor et al., 2011; Herman et al., 2012; Ceyhan-Birsoy et al., 2013; Begay et al., 2015; Oates et al., 2018; Ware and Cook, 2018). Several hypotheses have been proposed to explain the roles of nebulin and titin in muscle function and disease. For example, nebulin isoform size, associated with different types of skeletal muscle, has been proposed to regulate the length of thin filaments (Kruger et al., 1991; Labeit et al., 1991, 2011; Witt et al., 2006). Thus, in muscle deficient in nebulin, thin filament length is reduced and active tension is lower due to a reduction in thick–thin filament overlap (Witt et al., 2006; Ottenheijm et al., 2012; Li et al., 2015). Additionally, nebulin enhances myosin binding (Bang et al., 2009; Chandra et al., 2009; Ottenheijm et al., 2010), similarly to Tm (Maytum et al., 1999), and might play a role in regulating muscle contraction. In cardiac muscle, there are only low levels of nebulin present, marginalizing its effect on contraction (Kolb et al., 2016). Similar to nebulin-deficient muscles, cardiac muscle shows variable lengths of thin filaments (Robinson and Winegrad, 1977, 1979; Burgoyne et al., 2008) and lower isometric tension at full activation. On the other hand, titin elasticity, modulated by various mutations and splice variants, can strongly affect cardiac muscle function (Granzier et al., 2009; Chung et al., 2013; Methawasin et al., 2014, 2016; Hinson et al., 2015), for example length-dependent activation of cross-bridges (Wang and Fuchs, 1994; Metzger, 1995; Cazorla et al., 2001; Fukuda et al., 2003; Fukuda and Granzier, 2005; Ait-Mou et al., 2016).
The key contractile unit of striated muscle, denoted as the sarcomere (Fig. 1), integrates all the above structural elements and, depending on the activation level, defines changes in muscle tension and length. The muscle contractile tension per unit amount of actin–myosin overlap is determined by the number of attached cross-bridges and their state in the force-generating cycle. For a well-defined actomyosin cycle, the kinetics of myosin binding and detachment determines the density of bound cross-bridges. However, the number of attached cross-bridges per unit length of thin/thick filament overlap is tightly regulated by the calcium concentration and the associated kinetics of the interactions of regulatory proteins with actin. Higher concentrations of calcium and higher degrees of filament overlap strongly correlate with contractile force.
The degree of filament overlap is determined by the sarcomere length (SL) and the lengths of myosin and actin filaments. The length of myosin filaments is well conserved at ∼1.6 µm, but the length of actin filaments can vary from muscle type to muscle type. Mutations in nebulin or lack of nebulin will also produce large variations in actin filament length. For example, Witt et al. (2006) observed a large variation of actin filament lengths in tibialis cranialis muscle, isolated from a knockout (Neb-KO) mouse model, having a range from ∼0.4 to 1.2 µm, whereas in the WT muscles, the actin filament lengths are approximately constant at ∼1.2 µm. This large variation in length in Neb-KO muscle can explain, at least in part, the reduced isometric force generation and the leftward shift in the descending limb of the force–SL relations of Neb-KO fibers (Witt et al., 2006). Similar behavior is observed in patients with nebulin-based nemaline myopathy (Neb-NM; Ottenheijm et al., 2009). Additionally, force–pCa relations showed a lower sensitivity to calcium concentration in both Neb-KO and Neb-NM fibers (Chandra et al., 2009; Ottenheijm et al., 2010). However, it is unlikely that variable and, typically, shorter lengths of thin filaments in Neb-KO muscle can quantitatively explain all observations, and it may be necessary to consider other mechanisms such as reduced myosin cycling rates and reduced thin filament activation (Ottenheijm et al., 2012).
In cardiac muscle, nebulin is found only in miniscule amounts (Kazmierski et al., 2003; Kolb et al., 2016). Similarly to the range of actin filament lengths in nebulin-deficient skeletal muscle and the different actin filament length in different skeletal muscle types (Ottenheijm et al., 2012), in cardiac muscle the observed range of thin filament lengths is in atrial muscle ∼0.6–1.1 µm (Robinson and Winegrad, 1977) and in rat papillary muscle ∼0.9–1.1 µm (Robinson and Winegrad, 1977, 1979; Burgoyne et al., 2008).
The elastic properties of the sarcomere can also be modulated by mutations in structural proteins. The elastic properties of heart muscle are important for both diastolic filling and energy storage. In cardiac muscle, titin plays important roles in the structure and elasticity of the sarcomere (Trombitás et al., 1998; Granzier et al., 2002; Granzier and Labeit, 2006; King et al., 2011), and titin may affect the Frank–Starling relation in early diastole (Helmes et al., 2003). Notably, titin mutations can compromise elasticity and modulate cardiac cell function (Anderson et al., 2013; LeWinter and Granzier, 2013, 2014; Hinson et al., 2015). The cellular basis of this behavior is strongly associated with the change in calcium sensitivity with SL, resulting in a leftward shift of the force–pCa relations with an increase in SL (Wang and Fuchs, 1994; Metzger, 1995; Cazorla et al., 2001; Fukuda and Granzier, 2005). This increase in calcium sensitivity at longer SL, known as length-dependent activation, has been primarily attributed to the effect of changes in radial and longitudinal titin-based forces (Cazorla et al., 2001; Fukuda and Granzier, 2005; Fukuda et al., 2005; Methawasin et al., 2014; Ait-Mou et al., 2016). Recent studies have provided evidence that changes in the stiffness of titin can indeed play an important role in modulating calcium regulation via altering the structure of the thick filament (Ait-Mou et al., 2016). Several theoretical studies address the role of the elasticity of sarcomeric filamentous proteins in muscle contraction: (1) in a parametric model, Campbell (2009) estimated the effect of titin’s nonlinear elasticity on the heterogeneity of SLs along muscle fibers in response to a lengthening protocol; (2) in a spatially explicit 3-D model of the sarcomere, Fenwick et al. (2017) estimated the effects of cross-bridge and thin filament compliances on isometric tension and force–velocity relations; and (3) Powers et al. (2018) determined how titin elasticity modulates force–SL relations and ATP utilization rates. However, from all these studies, it is difficult to quantitatively separate the role of each of the possible contributing mechanisms to functional behavior at the organ level, or whether there are other mechanisms that have not been evaluated. Here we focused on the work by Cazorla et al. (2001), because this study contains the most complete data required for the simulations we present here. It includes multiple passive tension levels at the same SL (achieved by varying the history of prestretch before activation), it measured interfilament spacing at different SLs, and it reports the effect of lattice compression by dextran.
To quantitatively estimate the effect of nebulin deficiency on muscle function and the effect of titin passive elasticity and change in interfilament lattice spacing on length-dependent activation, we used the computational platform muscle simulation code (MUSICO). MUSICO takes into account the cross-bridge cycling kinetics, the explicit 3-D geometry of the myofilament lattice, including variation in actin filament lengths and titin based passive elasticity, and thin filament regulation involving a CFC for the Tm–Tn system (Smith, 2001; Smith et al., 2003; Geeves et al., 2011; Mijailovich et al., 2012a). In the calculations shown below, we determined the isometric active force (per myosin filament) for the physiological range of intracellular calcium concentrations and all conditions of geometry and muscle type. Specifically, the simulations include the lattice spacings between the thick and thin filaments, constant or variable thin filament lengths, the passive stiffness of titin molecules, and the observed slack lengths in skeletal and cardiac muscles. The roles of nebulin in decreasing isometric force and affecting its sensitivity to calcium are quantitatively tested by MUSICO simulations. The predictions differentiate the effects of changes in thin filament lengths, cross-bridge kinetics, and TnC sensitivity to calcium and match observed changes in nebulin-deficient muscles relative to WT muscles. Similarly, the MUSICO simulations distinguish the effects of mutation-induced changes in titin elasticity, interfilament spacings at different SLs, and the lattice compression by dextran. It is important to note that variable thin filament lengths are observed in skeletal muscles lacking nebulin due to mutations in the NEB gene and in cardiac muscle where nebulin is normally not expressed. It is warranted, therefore, to evaluate how thin filament length variation affects the contractility of both muscle types. The passive elasticity derived from titin is more relevant in the physiological SL range of cardiac muscle than in skeletal muscle because of the aforementioned shorter molecular spring region of cardiac muscle (Powers et al., 2018); thus, for simplicity, we thoroughly analyze only the role of titin in cardiac muscle. In so doing, the MUSICO simulations provide an opportunity to quantitatively assess and separate the multiple effects associated with structural and kinetics changes caused by mutations in nebulin and titin that can be experimentally measured only in the aggregate.
Results show that only part of the tension deficit observed in nebulin-deficient muscle can be explained by the effect of shorter lengths of actin filaments, but matching of the observed tension deficit also requires that nebulin enhances the rate of myosin binding to actin and alters the calcium sensitivity at different SLs. Similarly, the decrease in active tension and calcium sensitivity seen with decreased titin-based passive tension in cardiac muscle can be attributed to a partial redistribution of the muscle load between active muscle force and titin-based passive force, and by redistribution of strain along the thin filament. The nonuniform strain per se may cause significant changes in TnI-actin bond forces, modulating the release of TnI from actin, and hence, the degree of activation of the thin filament.
Materials and methods
The MUSICO simulation platform (Mijailovich et al., 2016) has a modular organization, allowing it to incorporate a spatially explicit 3-D sarcomere structure, any one of several actomyosin cycles (Mijailovich et al., 2016; Mijailovich, S.M., et al., 2015. Biophysical Society 59th Annual Meeting. Abstr. 337a; Mijailovich, S.M., et al., 2017. Biophysical Society 61st Annual Meeting. Abstr. 115a), and various models of thin filament regulation by calcium (Mijailovich et al., 2012a,b). This makes the MUSICO simulation platform well suited for investigating the effects of changes in sarcomeric proteins, their 3-D arrangement in the sarcomere, protein–protein interactions, thin filament regulation by Ca2+, and the roles of structural proteins such as nebulin and titin (Fig. 1). The basic sarcomere geometry minimally includes the fixed lengths of myosin and actin filaments, the specified slack length of a sarcomere, and the transverse interfilament spacing (Mijailovich et al., 2016). Each actin filament includes tropomyosin–troponin (TmTn) units aligned along each of two actin monomer strands, forming two continuous regulatory chains (Mijailovich et al., 2012a). The MUSICO platform now incorporates a coupling of thin filament regulation (Geeves et al., 2011; Mijailovich et al., 2012a) to our published 3-D explicit model of cross-bridge kinetics (Mijailovich et al., 2016) and the implementation of new features that include variable lengths of thin filament and (nonlinear) elastic contributions of titin to sarcomere contraction. In this section, we describe the broader concepts and put the detailed descriptions in the appendices. Appendix A describes the geometrical aspects of the 3-D interactions between myosin heads with their binding sites on actin. In Appendix B, we describe the formulation of strain-dependent rates in the actomyosin cycle. In Appendix C, we describe the basic concepts of thin filament regulation by calcium (Mijailovich et al., 2012a,b), and in Appendix D, the coupling between strain-dependent cross-bridge cycle and thin filament regulation by calcium.
3-D myofibril lattice structure
The minimal contractile unit in vertebrate muscle can be considered as a myofibril. For a typical myofibril of 1.2 µm in diameter, the estimated the number of myosin filaments from structural data (Millman, 1998; Irving et al., 2000) is ∼500–650. A myofibril runs the entire length of the muscle fiber and consists of a large number of sarcomeres in series. In a half-sarcomere, each half myosin filament with ∼150 myosin molecules faces six actin filaments; each actin filament has from 360 to 440 myosin binding sites, depending on the actin filament length, arranged in a double-stranded helix. The maximum overlap region is defined by the number of myosin crowns where 50 crowns interact with ∼512 active actin sites. The actin filament length as well as the number of actin monomers may vary due to mutations or presence or absence of nebulin.
Actin–myosin interactions are determined by the discrete lattice structure of interdigitated actin and myosin filaments (Fig. 1). Each myosin filament is decorated with crowns of myosin dimers, spaced by ∼14.3 nm along the filament where each crown consists of three myosin dimers with transverse orientations spaced by 120°, and successive crowns along the filament are rotated by +40° looking toward the Z-line or Z-disc (Fig. 2). An important parameter is the average interfilament spacing, d1,0, which may be measured from the 1,0 equatorial reflection in small-angle x-ray diffraction patterns (Higuchi et al., 1995; Millman, 1998).
The actin monomers in a thin filament form a double-stranded helix structure associated with the regulatory proteins Tm and Tn, with binding sites separated by ∼5.47 nm on each strand, with a half-period of ∼35.55 nm in the relaxed state (Huxley et al., 1994; Wakabayashi et al., 1994; Bordas et al., 1999; Prodanovic et al., 2016). The difference in periodicities between actin binding sites (∼35.55 nm) and myosin crowns (42.9 nm) creates a range of longitudinal spacings between myosin heads and actin binding sites, and this range of spacings is strongly modulated by the extensibility of myofilaments. In the 3-D sarcomere lattice, the actin binding domains on the myosin heads and the binding sites on actin require both longitudinal position matching and angular matching in the azimuthal plane (Fig. 3). A myosin head and the closest actin site form the most probable pair of these molecules that can create a cross-bridge interconnecting actin and myosin filaments (Mijailovich et al., 2016).
Mutations in, or lack of, nebulin, as observed in nemaline myopathy, are correlated with a large variation in thin filament lengths and are also associated with a reduced myosin binding rate to actin. Large variations in actin filament lengths are also observed in cardiac muscle due to low expression levels of nebulin. The passive elasticity of muscle due to titin can be compromised by mutations strongly modulating muscle contractility and regulation. To quantitatively assess the effects of changes in structural proteins titin and nebulin on muscle contractility, we implemented in MUSICO both the variable actin filament lengths due to mutations in nebulin and titin as a nonlinear spring connecting the Z-disc with the tip of myosin filament (Fig. 1).
Cross-bridge rate kinetics
Matching myosin heads with binding sites on actin
In the 3-D sarcomere lattice, the relative distance between a myosin head and the adjacent binding site on actin is defined by four factors (Mijailovich et al., 2016): (1) the axial displacements along myosin and actin filaments; (2) the transverse distance between myosin and actin filaments,
where d1,0 is the average interfilament lattice spacing derived from the 1,0 x-ray reflection; (3) the angle α defining the relative position of myosin to actin filament; and (4) the angle β defining how much a myosin head needs to turn along the actin perimeter of actin to reach an actin monomer in the correct orientation. These spatial distances and angles are essential information for formulating the strain dependence of the state transition rates, reflecting the discrete geometric relationships between myosin heads and actin binding sites (Fig. 3). The definition of the myosin binding step in relation to the 3-D sarcomere lattice is described in Appendix A.
Strain-dependent rates in the actomyosin cycle
The modular structure of the computational platform MUSICO can incorporate many different actomyosin cycles. To quantitatively assess the effect of titin and nebulin, without loss of generality, we considered a minimal three-state model that includes a swinging lever arm step or power stroke (Huxley and Simmons, 1971). For simplicity, we use here a three-state scheme (Mijailovich et al., 2016), similar to that proposed by Duke (1999) and Daniel et al. (1998), that includes one detached myosin state, M·ADP·Pi, and two attached myosin states: weakly bound, A·M·ADP·Pi, and strongly bound to actin after Pi release, A·M·ADP (Fig. S1).
The state transition rate constants between actomyosin states in cross-bridge models are derived from the functions for the free energy of the cross-bridge states in terms of the axial strain component x. In the 3-D sarcomere lattice, however, the strain-dependent rate of myosin binding to actin is modulated, in addition to the axial strain x (Fig. 3 A), by the spatial position of the actin filaments and the azimuthal departure of a myosin-binding site on actin filament from the plane passing through actin longitudinal axes (Fig. 3 C). This effect is implemented in MUSICO via weight factors Cα and Cβ, for the azimuthal departures and, in addition, includes the normalization factor fsk, that takes into account the number of actin sites that a cross-bridge can reach and bind (for further details see Appendix A and Mijailovich et al. ).
Calcium regulation of myosin binding to actin filaments
The contractility of striated muscle is regulated by the calcium-dependent azimuthal movements of Tm–Tn complexes over the surface of the actin filament. In the steric blocking model, the affinity of myosin for regulated F-actin is controlled by Tm, where Tm molecules are aligned along each strand of the actin double-stranded helix. An assumption common in most regulatory models is that each Tm molecule moves between three discrete orientations modulating actin affinities for myosin (McKillop and Geeves, 1993; Pirani et al., 2005; Poole et al., 2006).
A closer look at the Tm units shows structural evidence that neighboring tropomyosins overlap and that one end of TnT is bound to a specific site on Tm, whereas its N terminus overlaps the adjacent Tm, forming linked Tm–Tm regions. Thus, on each strand of F-actin, the interconnected neighboring TmTn units form the appearance of a CFC (Lorenz et al., 1993; Vibert et al., 1997) rather than a set of the independent TmTn units. The implementation of a CFC in sliding filament models provides structural, rather than ad hoc, insight into the observed cooperativity effects in models of thin filament regulation (Mijailovich et al., 2012a).
For modeling thin filament regulation as a CFC, we follow the approach of Smith and colleagues (Smith, 2001; Smith and Geeves, 2003; Smith et al., 2003; Geeves et al., 2011), and we have developed a Monte Carlo algorithm to quantify spatially explicit myosin binding to regulated F-actin in solution (Geeves et al., 2011; Mijailovich et al., 2012a). The algorithm consists of three main steps: (1) calculation of the state transitions between actin–TmTn states as a function of calcium concentration, (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 calcium concentration. The CFC position and its fluctuations modulate the original state transition rates of TnI or myosin binding to actin, whereas the spatial position of pinning sites, at which the TnI and myosin are bound to actin, determines the mean CFC angular position and the thermally induced azimuthal fluctuations (Geeves et al., 2011; Mijailovich et al., 2012a). 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 C.
Monte Carlo simulations of rate-dependent stochastic processes
In the stochastic model, we employed the standard Metropolis algorithm where a kinetic transition in a time step Δt occurs when a random number in (0,1) lies in the range (0,kΔt), where k is the first-order transition rate constant. This algorithm generates a Markov process if , so that at most one transition occurs per Monte Carlo time step in a single subsystem, here considered as two CFCs on one actin filament for TnI transitions or one myosin filament and its associated actin filaments for cross-bridge kinetics. This approach prevents interference between multiple transitions within a single subsystem, and the interference between the systems is negligibly small. Thus, Δt must be much smaller than the inverse of the fastest rate constant of the system, kmax, and in practice, kmaxΔt 0.001 was required to achieve satisfactory statistics.
The coupling between thin filament regulatory processes and the actomyosin cycle requires two sets of Monte Carlo random number drawings, and the overall numerical procedure within each time step includes (1) the first set of drawings defining transitions of TnI–actin states; (2) calculation of chain configuration with updated TnI–actin states and myosin bound states from the previous time step; (3) the second set of random number drawings defining the changes on actomyosin states regulated by the CFC; and (4) calculation of the mechanical equilibrium with external forces and constraints. The last step includes an iterative procedure to account for nonlinear elasticities included in the mechanical system, for example, the nonlinear elasticity of titin.
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 Δt. For each TnI or cross-bridge, the probability, in the range from 0 to 1, is divided into probability bins, Pij, 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 in which bin the drawn random number falls, the fate of a particular TnI or cross-bridge is defined and set for the following time step.
The calculation 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 in Mijailovich et al. (2012a) and Mijailovich et al. (2016), respectively. Here, we briefly describe the modulation of TnI–actin and myosin–actin transitions by accounting for the current configuration of myosin bound states and position of CFC at a prescribed concentration of calcium.
In a half sarcomere, each myosin half-filament interacts with six actin filaments, but because of symmetry, each half of a myosin filament effectively interacts with two actin filaments. The interactions between myosin heads (cross-bridges) and myosin binding sites on actin in the 3-D sarcomere lattice are as described in Mijailovich et al. (2016). Furthermore, each actin filament can be viewed as a double-stranded helix where each strand is associated with one CFC. We consider here that one CFC subsystem consists of one CFC and one actin strand. Because each TmTn unit covers seven actin monomers, the number of TnI-actin binding sites per strand is equal to one seventh of the number of actin monomers per strand (Mijailovich et al., 2012a). The coupling between cross-bridge cycling and calcium regulation of the thin filament is described in Appendix D.
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 (Figs. 1 and 3 B). Each half myosin filament is associated with six titin molecules. 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.735 nm, and the half period of one strand is 35.55 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 spacing provide maximum overlap with thin filament of ∼0.7 µm. The actin radius is ra = 3.5 nm, and myosin radius is rm = 7.8 nm (Fig. 3; 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 d1,0 are specified for each set of simulations.
Variability of actin filament lengths in skeletal muscle
It has been hypothesized that nebulin, a filamentous protein extending along the length of the actin filament, maintains constant length of thin filaments, but in disease, e.g., human nemaline myopathy, mutations in, or lack of, nebulin can result in much shorter randomly distributed thin filament lengths. To measure thin filament lengths, Witt et al. (2006) labeled skinned muscle fibers by streptavidin-coated gold beads and showed a uniform thin filament length distribution in WT muscle, ∼1.2 µm measured from the Z-disc (Fig. 4 A), whereas in Neb-KO sarcomeres, they observed shorter filaments ranging in length from ∼0.4 to 1.2 µm (Fig. 4 B). We define a complementary cumulative probability function of filament lengths, , from the function used for fitting of the histograms in Figs. 4, A and B (red lines) and normalization frequency, a. The distribution of actin filament lengths is the derivative of the cumulative probability function, , defined as (Figs. 4, C and D). The inverse function of defines the actin filament length as a function of cumulative probability, . The variable lengths of thin filaments are assigned in our spatially explicit MUSICO sarcomere geometry by substituting sequentially drawn random numbers, Rnd# in expression , where Rnd# is in range from 0 to 1 (Fig. 4 D).
Variability of lengths of actin filaments in cardiac muscle
The length of actin filaments in rat cardiac atrial trabeculae muscle also varies in range from ∼0.6 to 1.1 µm (Robinson and Winegrad, 1977). Using the best-fit distribution of actin lengths from Fig. 5 A, we obtained the cumulative probability function where La is the thin filament length (Fig. 5 B). The inverse function of defines the spatially random distribution the filament lengths,
for the assigned values for from sequentially drawn random numbers, Rnd#.
Interfilament lattice spacing
In skeletal muscle, lattice volume is approximately constant under ordinary physiological conditions. Thus, the interfilament lattice spaces can be calculated at any length from lattice volume data (Kurg et al., 1982; Millman, 1998), by the relation
where for mouse extensor digitorum longus muscle = 4.03 × 106 nm3 and sarcomere length SL is in the range from 2 to 4 µm. The skeletal muscle force–pCa simulations are performed at SLs of 2.5 and 2.0 µm, at which interfilament lattice spacing d1,0 is 37.36 and 41.77 nm, respectively. However, lattice volumes in cardiac muscle are not constant (Yagi et al., 2004) and need to be independently measured.
The interfilament spacings of skinned cardiac trabeculae muscle were obtained from experimental observations (Irving et al., 2000; Cazorla et al., 2001) at SLs of 2.3 and 2.0 µm, having d1,0 values of ∼40.9 and ∼42.9 nm, respectively. The spacings in trypsin treated muscles were larger, having values of ∼43.7 and ∼47.3 nm at SLs 2.3 and 1.9 µm, respectively, whereas they were much smaller after compression with 2.5% dextran, having values of ∼35.0 and ∼35.6 nm at SLs 2.3 and 2.0 µm, respectively (see also Table 1).
Actin and myosin filaments are extensible with filament moduli (elastic modulus times cross section area) derived from x-ray diffraction or direct measurement: for actin, Ka = 0.65 × 105 pN; and for myosin Km = 1.32 × 105 pN (Huxley et al., 1994; Kojima et al., 1994). We used here the modulus, AE, rather than stiffness, AE/L, because the reported stiffness values depend on the filament length, L, and cross section area A is not well defined for myosin and actin filaments. The presence or lack of nebulin can strongly affect the elasticity of the thin filament (Kiss et al., 2018), but a fourfold decrease in thin filament stiffness only mildly affects isometric contractile force in the steady state, i.e., 3%; thus only a single value of Ka is used in all simulations (Fig. S2). Similarly, the change in calcium sensitivity is small (ΔpCa 0.025).
The elasticity of titin in the physiological range is much more important in cardiac muscle than in skeletal muscle because of the shorter molecular spring of cardiac titin (Powers et al., 2018). For simplicity, we assumed that all passive tension transmitted to a myosin filament is contributed by titin. After recalculation of passive force per half myosin filament to force per titin molecule at SLs of 2.3 and 2.0 µm, we found an excellent agreement with the observations of Helmes et al. (1999) (Fig. 5 C). Because all calculations are done at isometric SL, the compounded stiffness of six titin molecules per half myosin filament is set to reach the passive force at the prescribed length and linearly changes around that point due to small deformations of the extensible myofilaments. At a slack length of 1.9 µm, the passive force and titin stiffness are assumed to be ∼0 (Granzier and Irving, 1995).
Cross-bridge model parameters
For simplicity, we used a three-state cycle that includes a detached state, a weakly bound state, and a strongly bound state. Following the approach of Duke (1999), the state transition rate constants are as follows: WT equilibrium constant for binding at SL = 2.5 µm is , where for skeletal and cardiac muscle values for the rate constants at appropriate SL and passive tensions are shown in Table 1; forward rate constant at zero cross-bridge strain for skeletal muscle at SL = 2.5 µm kbind is 40 s−1 in WT and 16.8 s−1 in Neb-KO; at SL = 2.0 µm, kbind is 25 s−1 in WT and 7.5 s−1 in Neb-KO; for cardiac muscle, kbind is given in Table 1; for the power stroke, the equilibrium constant is defined by for skeletal (WT) and for cardiac, whereas the power stroke d = 10.5 nm is the same for both types of muscle; and for ADP release/detachment, 150 s−1 and the second power stroke δ = 1.0 nm are taken to be the same in all simulations. The tensions for nebulin mutants in skeletal muscle and under various conditions in cardiac muscle are achieved by adjusting kbind and ΔGbind, but keeping k–bind at a constant value of 2 s−1. These binding rates strongly depend on interfilament spacing, and their values are assigned for each simulation. Because the exponential forms in the expressions for the state transition rates can become very large and can generate numerical problems, they are capped to 1,000 s−1, 100 s−1, and 104 s−1. These values are chosen as optimal values to satisfy Monte Carlo statistics for time steps on the order of 1 µs. When the cap value is reached, the reverse rates are changed to decay exponentially to satisfy the equilibrium constant, . In all simulations, cross-bridge stiffness is taken to be κ = 1.25 pN⋅nm−1 (as used by Duke ) and the value for = 4 pN·nm at 15°C (T = 288°K).
CFC model parameters
For the CFC model, we used the same parameters reported in Mijailovich et al. (2012a): a Tm pinning angle, ; myosin imposed Tm angular displacement, ; and the Tm persistence length estimated in multiple experiments ranges between 45 and 170 nm (Phillips and Chacko, 1996; Li et al., 2010a; Sousa et al., 2010; Loong et al., 2012). If we take into account a curved shape of Tm, which allows helical wrapping around without any strain, this could result in less flexible Tm having longer persistence length up to 450 nm (Li et al., 2010b), thus we adopted the value of (Wolgemuth and Sun, 2006; Mijailovich et al., 2012a) as a reasonable compromise. TnT interaction with Tm could have increased bending stiffness of Tm–Tn chain and show longer persistence length, as for example used in Mijailovich et al. (2012a). For other parameters, the angular Tm bending stiffness per unit length ΚCFC = 2.0 × 104 pN⋅nm4; the persistence length of Tm–Tn confined chain
as in Smith et al. (2003) and Mijailovich et al. (2012a); the angular standard deviation of the free chain σo = 29.7°; and the strength of the chain confining potential α = 0.341 pN. For cardiac muscle, the interaction of the free CFC chain with actin can be weaker and the stiffness of Tm–Tn chain stronger, so for cardiac muscle we used , which corresponds to a slightly wider σo = 35.7°. If ΚCFC is kept constant, the strength of the chain confining potential is weaker, having a value for αcp of 0.205 pN.
The level of thin filament activation is related to the calcium concentration, [Ca2+], and was derived from the approach of Mijailovich et al. (2012a) where the level of activation is defined by the TnI detachment rate from actin, , whereas the TnI reattachment rate is held constant at 100 s−1. The detachment rate is in large part associated with the binding of calcium to TnC, but the strength of the Tm interaction with the actin surface, the elasticity of CFC, and the cooperative effects of interactions between myosin binding and Tm chain could all play a role. To uncouple the enhancement of the activation level by the cooperativity induced by bound myosins to actin from [Ca2+]-induced dissociation of TnI from actin, we constructed for each simulation case the relationship between the equilibrium constant and the calcium concentration (Tables S1 and S2).
In the model of muscle activation, calcium binding is coupled with a change of TnI affinity to actin, and this interaction is defined by . The full description of calcium binding kinetics to TnC and transitions within Tm components that include the interaction of TnI with actin is essential for modeling muscle responses to calcium transients (e.g., twitch contractions). For simplicity, however, we did not go into the details of these transitions and we kept only essential information necessary for steady-state conditions. Because both simulations (F–KB) and observations (T–[Ca2+]) follow a Hill curve, log KB and pCa are approximately linearly related, with the slope equal to the ratios of the Hill coefficients, . The log KB–pCa relationships show a slight shift between the cases and reflect the differences in calcium sensitivity (Tables S1 and S2) observed under different experimental conditions.
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 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 to muscle function, but requires a large number of model parameters to properly define muscle fiber structure, cross-bridge kinetics, and thin filament regulation by calcium. We contend that conventional fitting of the data is not ideal, because the published experimental data show large variations and fitting many free parameters that can result in unreasonable estimates. Thus we used established structural data and best estimates deduced from experimentally established kinetic data for as many parameters as possible. In this study, we consider only steady-state isometric tension, in which case the number of variable parameters is reduced to two to three depending on the particular sets of data that model predictions have to match. Specifically, we adjusted only three parameters: (1) kbind, because change of interfilament spacing at different SLs and possible myosin activation by stretching thick filament via titin at longer SLs likely will modulate the myosin binding rate rather than rates of other steps on the cycle; (2) that defines interaction of TnI to actin associated with Ca2+ binding on TnC or release of Ca2+ from TnC, and (3) confined persistence length, , which includes elasticity of CFC and Tm interaction to actin and defines the cooperativity factor related to the Hill coefficient in force–pCa relations.
A sensitivity analysis is typically needed when fitted parameters have a large range of possible values (Mijailovich et al., 2010, 2012a; Ujfalusi et al., 2018), but in this study, we adjusted only three parameters (kbind, and ), and they were sufficient to match all experimental data. It is important to note that we adjusted these parameters using data from different experiments: kbind to achieve isometric tension at full activation where are irrelevant, and to achieve observed calcium sensitivity (pCa50), and confined persistence length to achieve Hill coefficient nH in force–pCa relations, where kbind is set from the previous fit. Because the adjusted parameters kbind and KB are decoupled, a classic sensitivity analysis seems unnecessary. Also, KB strongly affects pCa50, and the persistence length, , is primarily associated with the Hill coefficient nH, having only weak coupling; thus each parameter can be adjusted with minimal changes in the others.
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, F, 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 that is used in the tension calculations is measured.
In fibers from mouse skinned skeletal muscle, the lattice spacing, at slack SL of 2.25 µm, is d1,0 = 39.38 (Kurg et al., 1982; Millman, 1998), and the corresponding number of thick filaments per square micrometer of myofibrillar area is ∼562. In most species and skeletal muscles, the myofibrils occupy from 80 to 85% of the fiber volume (Eisenberg, 1983). Assuming the fraction of cross section occupied by myofibrils as 0.83 (Mobley and Eisenberg, 1975), the number of thick filaments per square micrometer in cross section of skinned muscle cell is ∼466. The scaling factor provides a relationship between isometric tension, To, and average force per myosin filament, Fo, where, for example in skeletal muscle,To = 125 kPa corresponds to per myosin filament.
The lattice spacing in skinned 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 square micrometer myofibril is ∼458. Taking into account the fraction of the cross-sectional area in cardiac cell occupied by myofibrils as 68% (Cazorla et al., 2000), this reduces the number of myosin filaments over the muscle cell cross section to ∼311/μm2. In this case, the scaling factor provides that To = 40 kPa corresponds to Fo = 128.5 pN/myosin filament or 64.25 pN/actin filament.
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 up to ±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, for the number of myosin filaments used in simulations that matches the number of the 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, ∼1/3 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 having to run the simulation multiple times. In the final simulations, we increased the number of myosin filaments per half sarcomere to 500.
A myofibril shows an incomplete hexagonal lattice at the periphery. By keeping the ratio of the number of half myosin filaments to the number of thin filaments to 1:2, the error in calculations in an incomplete lattice is associated only with the effect of thick and thin filament extensibility. In steady state, the effect of filament extensibility is small and is related to matching-mismatching of cross-bridges with myosin binding sites on actin filaments, resulting in minor changes in the number of bound cross-bridges.
The interaction of myofibrils or muscle fibers with other components of muscle, such as collagen, and components of muscle cells aligned parallel to the fibers are not taken in account, because in the current simulations, we considered only active force and tension under isometric steady-state conditions.
The MUSICO software environment and simulation details
MUSICO software has been developed as a C++ object-oriented application that includes a LAPACK linear algebra package and deal.II finite element library. 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 takes ∼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 describes the three-state actomyosin cycle (Fig. S1), the formulation of strain-dependent state transition rates; relationship between thin filament regulation state transition rates and calcium sensitivity of isometric force in WT and nebulin-deficient skeletal muscles (Table S1); the effects of passive tension on calcium sensitivity in cardiac muscle (Table S2); the prospective decrease in isometric force and muscle stiffness (Fig. S2) and in calcium sensitivity (Fig. S3) due to the increase in actin filament compliance that appears in nebulin-deficient actin filaments (Kiss et al., 2018); and normalized force–pCa relation in WT versus Nebulin-deficient muscles (Fig. S4) and at high and low passive tension (HPT and LPT; Fig. S5).
The contractility of muscle modulated by nebulin
Using MUSICO, we evaluated whether the force deficit in nebulin-based nemaline myopathy can be explained by (a) shorter and variable thin filament lengths, (b) altered cross-bridge cycling kinetics, and (c) reduced myofilament calcium sensitivity. We quantitatively estimated the effect of nebulin deficiency on muscle function in mice using MUSICO simulations, where we compared the active force (per myosin filament) in normal (WT) and nebulin-deficient (Neb-KO) muscles. Predicted force–length relationships of WT muscle, with normal sarcomere geometry including uniform distribution of the length of actin filament filaments (∼1.2 µm), and for Neb-KO muscle, with shorter and variable length of actin filaments, are shown in Fig. 6 A. These predictions are compared with observations shown as green filled circles for WT and red filled triangles for Neb-KO (Ottenheijm et al., 2009). The predicted force–length relation for WT (pink line with triangle symbols) closely follows the observations (green filled circles). Ottenheijm et al. (2009) did not measure tension of either WT or Neb-KO at SLs 2.2 µm, because at these lengths, thin filaments overlap with those from the adjacent half sarcomere and the striation pattern is compromised, making SL measurement unreliable. In the MUSICO simulations, we took into account the decrease of tension at SLs where thin filaments overlap each other, partially preventing myosin binding. In WT fibers, at SLs 2.25 µm, the decrease in tension or force per myosin filament is caused by thin–thin filament overlap, where thin filament passing the M-line partially blocks the interaction of thick filaments, reducing the effective myosin binding rate. The resulting tension at a SL of 2.0 µm is ∼12.5% less than the plateau tension.
The nebulin-deficient sarcomere with variable lengths of the thin filaments showed a leftward shift of the force–length relation (gray line with square symbols), but the shift is insufficient to reach the experimentally determined values in Neb-KO fibers (red filled triangles). Thin filaments also overlap for SL 2.6 µm in Neb-KO fibers because the longest actin filaments can reach lengths up to 1.4 µm. However, because of the small number of thin filaments involved in this process, the decrease in tension at 2.0 µm is only ∼5% compared with the tension with no thin–thin filament overlap. Interestingly, the increase of the average overlap between myosin filament and shorter thin filaments increases the corresponding tension by a much larger amount than the decrease caused by thin–thin filament overlap; thus the tension continues to increase with further decreases in SL. If only the variable thin filament lengths of Neb-KO are taken into account, the predicted tension with the same actomyosin kinetics (gray line with filled square symbols) slightly exceeds WT tension at SL ∼2.0 µm. This MUSICO prediction is similar to the estimates of Chandra et al. (2009). To obtain large decreases in the magnitude of force and a much larger vertical shift of the force–length relation requires decreasing the myosin binding rate, kbind, by more than twofold (blue line with cyan filled triangles). This large decrease in kbind can be attributed, in part, to loss of enhancement of the binding rate by nebulin. Collectively, these data demonstrate that both reduced thin filament length in nebulin-deficient muscle and reduced binding rates are necessary to explain the observations in fully activated muscles (Fig. 6 A).
The effect of nebulin deficiency on calcium sensitivity is revealed by the force–pCa curves shown in Fig. 6 B. Again, the predicted force–pCa curve for WT (pink line with triangle symbols) closely follows the observations (green filled circles), whereas the Neb-KO sarcomere showed a drop of ∼25% at full activation (gray line with filled square symbols) and negligible change in sensitivity to [Ca2+] (Fig. S4). This decrease is less than half of the decrease observed in Neb-KO muscles (red triangles; Witt et al., 2006). Decreasing the myosin binding rate, kbind, by ∼2.4 times, with no change of kinetics of calcium binding to TnC and associated TnI affinity to actin, enabled matching the simulation to the observed force at full activation (blue dotted line with hollow blue triangles), but showed lower calcium sensitivity by ΔpCa 0.116 compared with the observations (Witt et al., 2006). A large part of this loss of calcium sensitivity (ΔpCa 0.090) can be attributed to the reduced cooperativity between the CFC chain and a reduced number of bound myosins (blue arrow in Figs. 6 B and S4 A). To simulate the observations (red filled triangles) it was necessary to increase the detachment rate of TnI from the thin filament, , denoted by the equilibrium constant . Because all simulations are done under isometric conditions and at steady state, the observed three- to fourfold decrease in thin filament stiffness in Neb-KO muscles (Kiss et al., 2018) results in only a small effect on the isometric tension (or force per myosin filament), i.e., a decrease of up to 3% (Fig. S2). Similarly, the net change in sensitivity for the range in thin filament compliances is also small (Fig. S3, A and B): ΔpCa50 ∼ 0.038 when thin filament stiffness changed from 65 pN⋅nm−1 reported by Kojima et al. (1994) and 34 pN⋅nm−1 estimated in WT soleus muscle thin filaments from small-angle x-ray diffraction studies (Kiss et al., 2018); the net effect of the decrease in thin filament stiffness between WT and Neb-KO deficient thin filaments (by a factor of 3) on the calcium sensitivity is even smaller, ΔpCa 0.025.
Follow-up mechanical studies performed at an SL of 2.0 µm (Chandra et al., 2009) confirmed that nebulin significantly enhances tension generation (Fig. 6 C) by a relative value similar to that reported earlier (Witt et al., 2006; Figs. 6, A and B), but the isometric tensions in both fully activated WT and Neb-KO muscles were significantly lower. This reduction in tension can be only partially explained by the tension reduction due to the thin–thin filament overlap and likely is associated with altered binding kinetics, presumably caused by the differences in experimental conditions. For example, a change of SL from 2.3 to 2.0 µm increases lattice interfilament spacing (Kurg et al., 1982; Millman, 1998), and it is likely associated with a reduced myosin binding rate or, alternatively, with other mechanisms associated with length-dependent activation. Chandra et al. (2009) showed that nebulin decreases calcium sensitivity, unlike the study by Witt et al. (2006) that found no difference. To quantitatively explore the origins of these differences, we performed MUSICO simulations at an SL of 2.0 µm and adjusted the kinetics to match observations of Chandra et al. (2009) (Fig. 6 C).
The MUSICO-predicted tensions in WT muscle fibers, taking into account only the variability of thin filament lengths at full activation, showed tension values similar to those in WT fibers (Fig. 6 C). These similar tensions parallel MUSICO predictions for the force–length relation at SL = 2.0 µm in Fig. 6 A, where a significant decrease in WT tension due to the thin–thin filament double overlap is compensated by an increase in Neb-KO tension due to the increase of overall overlap between variable-length thin filaments and thick filaments. This effect is reflected in an almost identical tension–pCa relation for WT and Neb-KO if the only differences considered are in the thin filament length distributions (Fig. 6 C, pink line with triangle symbols vs. gray line with filled squares).
Because at an SL of 2.0 µm there is almost no effect of thin filament length variability, to achieve the observed drop in tension in Neb-KO requires a decrease in myosin binding rate, kbind, by ∼3.3 times to match the observed tension that is ∼45% of the tension at full activation (Chandra et al., 2009). However, if there is no change in the kinetics of calcium binding to TnC and associated TnI affinity to actin, the observed force at full activation is matched by the above decreased kbind (blue dotted line with hollow blue triangles), but this results in higher calcium sensitivity by ΔpCa 0.075 relative to the observations (red filled triangles). The loss of calcium sensitivity by decrease of kbind is due to reduced cooperativity between CFC chain and a lower number of bound myosins of ΔpCa 0.083 but only covers about a half of the observed loss of sensitivity (see also Fig. S4 B). This is opposite to that shown in Fig. 6 B, indicating that to explain the observations, it is necessary to decrease the calcium sensitivity of Neb-KO by, for example, a decrease in TnI detachment the rate from the actin filament, . In this case, the Hill coefficient, nH, of ∼3.55 for Neb-KO is slightly lower than ∼3.9 for WT, indicating a moderate loss of cooperative interactions between actin filaments, myosin, Tm, Tn, and calcium. This behavior could be affected by the use of reconstituted fiber bundles with fast skeletal muscle recombinant Tn in Chandra et al. (2009). This major structural change could in part contribute to the decreased sensitivity of Neb-KO to calcium and loss of cooperativity.
The contractility of muscle modulated by titin
Using MUSICO, we quantitatively determined the effect of titin passive elasticity and change in interfilament lattice spacing on myofilament length-dependent activation. The geometrical model for all cardiac muscles includes the observed variation in actin filament lengths and nonlinear (passive) elasticity of titin. We compared the calcium sensitivity of active force (per myosin filament) at various levels of passive force by changing SL (2.3 vs. 2.0 µm), using a large prestretch to 2.5 µm and after holding, releasing to 2.3 or 2.0 µm, and then activating, denoted here as the mechanical protocol. Alternatively, the passive tension was reduced by trypsin treatment. We simulated the effect of lattice spacing–induced changes in the force–pCa relation. All MUSICO parameters affected by the variation in protocols are shown in Table 1.
Fig. 7 shows that reduction in titin-based passive tension, by either mechanical prestretch and release (Fig. 7 A) or lowering titin passive tension by trypsin (Fig. 7 B), reduces the magnitude of the isometric force in fully activated muscle fibers (pCa = 4.5) and also decreases sensitivity to calcium. The change in isometric force at full activation induced by reducing passive tension by decreasing SL from 2.3 to 2.0 µm is larger than the mechanical reduction of passive tension at 2.3 µm. An excellent fit can be achieved by reducing the rate of myosin binding to actin, kbind, and this change is correlated with the reduction in passive tension (Table 1). The decrease in sensitivity is caused by an increased rate of calcium dissociation from TnC and a reduced affinity of TnI to actin by passive forces, i.e., an increase of TnI detachment rate from actin, . Compression of the myofilament lattice by dextran strongly increases myosin binding rates, and the force at full activation is 10% higher (Fig. 7 C).
The input data and the predicted isometric forces in fully activated muscle are summarized in Table 1. The observed decrease in the isometric active force at a fixed SL can be achieved by decreasing the binding rate, kbind. Changing SL from 2.3 to 2.0 µm causes a decrease of isometric active force due to a decrease in passive force and an increase in interfilament lattice spacing. In this case, decreases in kbind are due to an increase in interfilament lattice spacing, i.e., radial strains. The opposite effect is shown when interfilament spacing is decreased by lattice compression by dextran. In this case, the observed increase in active tension is achieved by an increase in kbind. For simplicity, in all simulations shown in Fig. 7, we fixed k–bind to 7.41 s−1 and only changed kbind.
Effect of titin-based passive tension and interfilament spacing on calcium sensitivity
The change in sensitivity, ΔpCa50 (from the control at SL = 2.0 µm), with passive tension is correlated with the increase in active force, but the change is much smaller in relative terms, indicating that additional mechanisms contribute to the disproportionate increase in ΔpCa50 (Table S2). The LPT induced by trypsin treatment showed only minor differences compared with LPT (mechanical protocol) at SL = 2.3 µm, and the changes in ΔpCa50 and were in the same proportion, suggesting that changes in binding rates are sufficient to explain changes in sensitivities relative to the control at SL = 2.0 µm.
After compression by dextran, the change in force, ΔpCa50, and magnitude of from control at SL = 2.0 µm are well correlated with passive tension, but they have different changes in the magnitude of the intercepts at zero passive tension. Importantly, ΔpCa50 increases faster than and has a lower relative intercept at zero passive tension, indicating that several different mechanisms affect these changes, but to a much lesser degree than predicted from the mechanical protocols (Fig. S5). Comparison of the LPT (control) with LPT after compression with dextran at SL = 2.0 µm shows large shifts in ΔpCa50 of 0.258 compared with the much smaller change in , indicating that large decreases in interfilament lattice spacing by dextran produce large changes in sensitivity.
The MUSICO simulations, now including the CFC model (Geeves et al., 2011; Mijailovich et al., 2012a), predict force–pCa relationships. This makes it possible to test hypotheses related to changes in calcium sensitivity under various experimental conditions including compromised contractility due to mutations in key sarcomere proteins. The inclusion of titin and nebulin in MUSICO provides an opportunity to quantitatively separate effects that cannot be measured directly, allowing validation of hypotheses concerning the roles of these proteins and potentially explaining the functional effects of mutations in titin and nebulin. This modeling approach takes into account known geometrical factors, for example, the thin filament length variability in nebulin-deficient skeletal muscles or in cardiac muscles, changes in titin elasticity, altered interfilament lattice spacing induced by changes in SL, or lattice compression by dextran. Each MUSICO simulation includes the specific conditions associated with each experiment examined, including calcium concentration and the passive force changes that occur with altered SL.
Role of nebulin in contractility of skeletal muscle
It has been hypothesized that the principal role of nebulin is to regulate the length of thin filaments (Kruger et al., 1991; Labeit et al., 1991; Witt et al., 2006). Witt et al. (2006) observed the lengths of thin filaments in fast skeletal muscle of normal (WT) and nebulin-free mice (also denoted as Neb-KO; Witt et al., 2006). They measured the effect of this deletion on force–length relations in fully activated muscles and force–pCa relations in partially activated skinned muscles. This study shows similar behavior to that observed in the quadriceps muscles of nemaline myopathy patients (Ottenheijm et al., 2009, 2012), but from the data presented it was difficult to quantitatively separate the effects of the variability of thin filament lengths observed in Neb-KO, changes in force–pCa relations, and other functional roles of nebulin (Chandra et al., 2009; Ottenheijm et al., 2010).
Observations from multiple experimental studies in humans (Ottenheijm et al., 2008, 2012) and murine muscles (Witt et al., 2006; Chandra et al., 2009; Ottenheijm et al., 2012; Li et al., 2015) suggest three principal factors in nebulin-deficient muscles: (1) changes in thin filament lengths that strongly modulate the overlap between thin and thick filaments and therefore the generation of tension (Witt et al., 2006; Ottenheijm et al., 2012); (2) changes in cross-bridge cycling kinetics (Bang et al., 2009; Chandra et al., 2009; Ottenheijm et al., 2010, 2012) that also affect generation in muscle tension; (3) changes in muscle regulation reflected in altered sensitivity to [Ca2+] (Chandra et al., 2009; Ottenheijm et al., 2012). The murine studies in WT and nebulin KOs reveal the symptoms observed in nemaline myopathy, such as muscle weakness (Ottenheijm et al., 2012), but cannot, on their own, reveal the quantitative distinctions between the effect of sarcomere structural changes, compromised binding kinetics, and modulation of thin filament regulation by calcium. The MUSICO platform includes the changes in sarcomere geometry between WT and nebulin-deficient muscles and estimates the changes in cycling kinetics in fully activated muscle required to match the observations and quantify the changes in thin filament regulation kinetics to match observed force–pCa relations. Thus MUSICO simulations enable us to quantitatively separate the contributions of the above changes to force and calcium sensitivity that can be experimentally observed only in the aggregate.
Effect of variability of thin filament lengths and cross-bridge cycling kinetics
Nebulin-deficient muscle lacks the characteristic plateau in the force–length relationship in both murine nebulin KOs and humans with nemaline myopathies (Ottenheijm et al., 2009, 2012). The MUSICO predictions for WT muscles showed a plateau at optimal thick-thin filament overlap and a descending limb at SLs exceeding ∼2.58 µm (Fig. 6 A), reflecting the approximately constant thin filament length (Fig. 4, A and C) and the length of the bare zone in the thick filament. At shorter SLs (2.2 µm), tension is significantly reduced due to impaired myosin binding in regions where thin filaments with opposite orientations overlap. The effective myosin binding rate in this region based on simulation is reduced by ∼65%, causing a reduction in tension at SL of 2.0 µm of ∼12.5% from the tension at the plateau. In contrast, the MUSICO simulations including only observed shorter and nonuniformly distributed filament lengths in Neb-KO (as in Fig. 4, B and D) predicted a continuous increase in isometric tension with a decrease in SLs (Fig. 6 A). Moreover, because the thin–thin filament overlap interference showed only a modest decrease in force, the tension even exceeds the isometric tension at SL = 2 µm. Thus, the simulations including only the nonuniform distribution in length of thin filaments in Neb-KO were not sufficient to explain the observed reductions in the force–length relation, relative to WT muscle. These simulations showed only ∼20% drop in isometric force at full thick-thin filament overlap and a leftward shift in the descending limb of force–SL relation of about a half of that observed. Therefore, other mechanisms likely operate as well (red triangles in Fig. 6 A). For example, decreasing cross-bridge cycling rates via kbind by ∼2.4 times provides good agreement with the observation (blue line with cyan filled triangles vs. red triangles in Fig. 6 A). This decrease of kbind in Neb-KO muscles likely reflects loss of the enhancement of myosin binding to actin in presence of nebulin.
There could be other possible mechanisms, however, to explain the observations. An alternative hypothesis is that the drop in tension in nebulin-deficient muscles is caused by an increase in compliance in nebulin-deficient muscle with no change in cross-bridge cycling rates (Kawai et al., 2018). However, MUSICO simulations showed only an ∼3% drop in the isometric tension (Fig. S2) that can be attributed to the increase in thin filament compliance observed by Kiss et al. (2018). On the other hand, nebulin deficiency could also cause a reduced force per attached cross-bridge (Kawai et al., 2018), a possibility that requires future study.
Myofilament calcium sensitivity
Force–pCa relations were used to identify differences in the calcium sensitivity of force production. In most of the studies, a decrease in calcium sensitivity in nebulin-deficient muscle fibers was reported (Chandra et al., 2009; Ottenheijm et al., 2010, 2012). But other studies reported no difference in calcium sensitivity between WT and Neb-KO (Witt et al., 2006; Bang et al., 2009), thus bringing some uncertainty in the observations from the former studies. The differences in sensitivity were attributed to the differences in experimental procedures or, alternatively, to the difference in SLs used in the studies (Ottenheijm et al., 2012).
The MUSICO simulations enable a quantitative separation of the effects of variable SLs, the change in cross-bridge cycling kinetics, and thin filament regulation. The simulations reveal that although the maximum isometric forces at full activation match the predicted force for WT and, after decreasing kbind by ∼2.4 times, for Neb-KO, as shown in Fig. 6 A, the predicted force–pCa relation for Neb-KO is shifted rightward compared with Witt et al. (2006) (blue dotted line with open triangles in Figs. 6 B and S4 A). The shift, i.e., loss of calcium sensitivity, is caused by reduced cooperativity by a lower number of bound myosins interacting with CFC. The predicted lower sensitivity needs to be compensated for by an increase in activation level via an increase of to match the observations. These studies were performed at an SL of 2.5 µm, where in WT sarcomeres, thick-thin filament overlap is optimal and there is no force reduction due to thin–thin filament overlap. In studies performed at shorter SLs (Chandra et al., 2009), however, thin–thin filament interference at SL of 2.0 µm, and probably other effects, decreased the isometric force at full activation in both WT and Neb-KO. To match the tension in Neb-KO, kbind should be decreased by ∼3.3-fold. In this case, the loss of calcium sensitivity caused by reduced cooperativity has a similar magnitude (ΔpCa50 0.083 at SL = 2.0 µm vs. 0.090 at SL = 2.5 µm as shown in Fig. 6, C and B, respectively). However, the decrease in sensitivity at SL = 2.0 µm is not sufficient to match the observation, and the simulation predicted higher sensitivity (blue dotted line in Fig. 6 C) than observed by Chandra et al. (2009) (red triangles). Thus to fit the observations, in the first case, at SL = 2.5 µm, it is necessary to increase sensitivity by increasing TnI detachment from the thin filament (KB arrow to the left in Figs. 6 B and S4 A) and in the second case, at SL = 2.0 µm, the decrease in sensitivity requires a decrease in TnI detachment from the thin filament (shown as KB arrow to the right in Figs. 6 C and S4 B).
In a recent publication, Kiss et al. (2018) observed a decrease in thin filament stiffness in Neb-KO of 2.5–3 times compared with WT. Our simulations show only small changes in sensitivity with a ΔpCa of ∼0.038 for the experimentally observed changes in thin filament stiffness of threefold, assuming that the thin filament length in WT and neb-KO are uniform (Fig. S3). The lengths of thin filaments in Neb-KO however, are known to be nonuniform and shorter than in WT; thus the effect of a threefold decrease in thin filament stiffness on calcium sensitivity is even smaller, ΔpCa 0.025. Because this change of sensitivity is more than an order of magnitude smaller than observed, and it is neglected in our analysis. This change in sensitivity is much smaller than the change in sensitivity due to altered thin filament stiffness reported by Chase et al. (2004). The large change in sensitivity estimated by Chase et al. (2004) is likely a consequence of the large range in thin filament stiffness that was simulated, up to 100-fold, that far exceeds the stiffness changes observed by Kiss et al. (2018). It seems, then, that the changes in calcium sensitivity in nebulin-deficient muscles are primarily due to changes in cross-bridge cycling kinetics and changes in the detachment rates of TnI.
Role of passive tension in contractility of cardiac muscle
In cardiac muscle, the role of passive tension is important for understanding the SL dependence of calcium activation, which is the underlying basis of the Frank–Starling mechanism (Kentish et al., 1986). Titin is identified as one of the key contributors to myocyte passive force and length-dependent calcium sensitivity (Granzier and Irving, 1995; Cazorla et al., 2001). It was hypothesized that the observed increase in calcium sensitivity could be triggered by an increase in affinity of TnC for calcium at larger SLs and cooperative effects associated with increase in the number of strongly bound cross-bridges (Cazorla et al., 2001). MUSICO simulations can estimate quantitative differences in the magnitude of isometric force and calcium sensitivity at different levels of titin passive forces, changes in the interfilament lattice spacing, and modulation in cycling rates in cross-bridge cycle, allowing these factors to be evaluated separately.
Effect of cross-bridge binding rates and interfilament lattice spacing on isometric force in fully activated cardiac muscle
Cazorla et al. (2001) reported isometric tensions at full activation and force–pCa relations at HPT and LPT. This was achieved by changing SL, by a mechanical protocol; by degradation of titin by trypsin; or by changing lattice interfilament spacing with dextran. Because the changes in SL and treatment by trypsin change the lattice spacings, and the decrease of SL increases the effective thin–thick filament overlap due to variable thin filament lengths in cardiac muscle, it is difficult to extract the net effect of these multiple changes from the observations that provide aggregate data. Thus, let us first start with a comparative analysis at each SL and then examine the effect of changes in SL.
At SL = 2.3 µm, mechanical reduction of the passive tension decreases the active force by ∼8.1% from the level at HPT (Table 1). In this case, there is no change in interfilament lattice spacing, but a decrease of kbind by ∼15.2% is required to match the observations (Cazorla et al., 2001). In contrast, at SL = 2.0 µm, the tension in titin is small, so the mechanical change from HPT to LPT at this SL shows negligible change in active force and no change in either kbind or lattice spacing. This comparison suggests that the decrease of kbind at SL = 2.3 µm can be attributed solely to the change from HPT to LPT.
Reducing passive tension by trypsin at SL = 2.3 µm decreases the active force by a much smaller amount (3.5%), compared with HPT, and the decrease in kbind is also small (∼4.3%) suggesting that the same decrease in passive tension has a much greater effect on active force when achieved by the mechanical maneuver than by the trypsin protocol. Comparing the increase in active tension from mechanical (protocol) LPT to trypsin LPT, also at SL = 2.3 µm, shows an increase in active force of 5%, in kbind of 12.9%, and in lattice spacing of ∼6.8%. It seems that degradation of titin not only reduces passive tension but also enhances kbind at LPT and the magnitude of the isometric force. The enhancement of the isometric force and kbind by degradation of titin might be attributed to some other mechanism than increase in lattice spacing, because it is widely anticipated that both kbind, and the force, decrease with an increase in spacing.
Compression of the sarcomere lattice by dextran at SL = 2.3 µm showed a spacing decrease of ∼14.5%, an increase in active force by 6.8%, and an increase in kbind by ∼29.1% at HPT, with a larger increase in active force of ∼10.7% and kbind of ∼37.6% at LPT induced by the mechanical protocol (Table 1). At short SLs (SL = 2.0 µm), addition of dextran at LPT after mechanical manipulation decreases interfilament spacing by 16.9%, but the enhancement of force is much smaller, an increase of 6.1%, accompanied by relatively small increases in the binding rate of only 6.7%.
If we now return to the change in force at decreasing SL, we can anticipate that a decrease in SL increases the degree of thin–thick filament overlap, but myosin binding rates decrease due in part to an increase in lattice interfilament spacing. To dissect the contribution of the degree of overlap versus the density of bound cross-bridges to the developed tension, we quantified the magnitude of each contribution separately.
The nonuniform distribution of actin filament lengths (Fig. 5 A) results in a distribution of thick–thin filament overlap lengths and an increase in isometric active force at shorter SLs. At full activation and the same kinetic state transition rates, the force at SL = 2.3 µm is only ∼66% of the force that could be generated if the thin filament lengths were uniform and 1.1 µm long. At shorter SLs, the force will increase to ∼83% at SL = 2.0 µm and ∼87% at SL = 1.9 µm of the maximum force (Fig. 8). Note that if all thin filaments are 1.1 µm long, the active force does not change in the above SL range because all the thick and thin filaments fully overlap. On the other hand, thin–thin filament overlap appears at SLs of 2.0 and 1.9 µm, but because only a small number of thin filaments with nonuniformly distributed lengths experience this overlap, the decrease of force is small (1%) and has a negligible effect on the isometric forces.
Overall, due to the effective increase of thick-thin filament overlap with decreasing SL, the magnitude of active force is expected to increase unless it is counterbalanced by decreases in kbind to meet the observations of Cazorla et al. (2001). As we described above, the decrease in kbind can be attributed only in part to increases in interfilament spacing (Table 1), and other mechanisms likely affect kbind as well.
The decrease in tension with a decrease in SL from 2.3 to 2.0 µm, where passive tension changes from HPT to LPT, can be predicted by decreasing the binding rate kbind by 42.8% (Table 1). After the decrease in passive tension by the mechanical protocol (LPT (Mech); Table 1), the decrease in the sarcomere length decreases kbind by 32.6%. In both cases, the increase in lattice interfilament spacing is ∼4.6% and is weakly associated with large decreases in kbind. Because the effective overlap increases with a decrease in SL, and a decrease in kbind affects the active force per unit of overlap, the calculated decrease of active force is proportional to difference of the products of kbind and the degree of overlap at SLs 2.3 and 2.0 µm. Thus, to extract the net effect of the decrease of kbind on the active force, the force should be normalized to the full overlap length that can be achieved with 1.1-µm-long thin filaments at all three SLs as shown in Fig. 8. Due to the distribution of thin filament lengths, the effective thick–thin filament overlap at SLs 2.3 and 2.0 µm is 66 and 83% of the full overlap, respectively (Fig. 8). Thus, after normalization of the net decrease in active force with the decrease in sarcomere length for change from HPT to LPT is 10.5%, and the decrease in force with the same decrease in length at LPT (after mechanical protocol) is only 2.6%. These data suggest that passive tension has a much stronger effect on the change of normalized force than does the interfilament lattice spacing.
In trypsinized preparations, the passive tension is low, showing minor differences with change of length. But, in this case, lattice spacing increases by 8.4% and kbind decreases by 67.7% for a change of length from 2.3 and 1.9 µm. These much larger changes in kbind and in lattice spacing than in the absence of trypsin are associated with a much larger reduction in normalized net active force, by 66.7%. Because the change in passive tension is small, the change in normalized active force can be almost entirely attributed to the decrease in kbind that is likely strongly associated with the large increase in lattice spacing.
In contrast, in the sarcomere lattice compressed by dextran, the increase in lattice spacing with a decrease in SL from SL of 2.3 to 2.0 µm is much smaller (∼1.7%) than in the trypsin case, and it is likely only weakly associated with the large decrease in binding rate kbind by 47.7% at LPT (mechanical protocol) and by 55.0% for change from HPT to LPT by decreasing SL only (Table 1). The net active force, after normalization, decreases by 26.2% at LPT (mechanical protocol) and by 29.6% for change from HPT to LPT by decreasing SL only. This change parallels the decrease in kbind.
Effect of titin-based passive tension on myofilament calcium sensitivity
MUSICO predictions show similar values for the Hill coefficient, nH, and ΔpCa as observed (Cazorla et al., 2001) and confirm an increase of calcium sensitivity with increases in passive tension (Table S2). The same trend is also obtained from the calculated KB–pCa relations, suggesting that the kinetics of binding and detachment of TnI from actin are synchronized, in part, with calcium binding to TnC. Typically, a decrease in the rate of calcium dissociation from TnC, i.e., an increase in the equilibrium rate, KB, is associated with a leftward shift of pCa curves, i.e., toward higher sensitivity, and can be expressed in terms of ΔpCa50. However, changes of the Hill coefficient, nH, in the force–pCa relation and associated changes in ηH in the KB–pCa relations indicate that the cooperative interactions between myosin binding and Tm chain have a strong effect on calcium sensitivity. There are likely other contributions to myofilament calcium sensitivity. Ait-Mou et al. (2016) proposed a model where myofilament length-dependent activation involves, at least in part, transduction of titin-based passive tension into changes in the ordering of the myosin heads, and some communicating structure transmits strain to the thick filaments, leading to an alteration in the structure of Tn and, presumably, the degree of activation of the thick filament. In several papers from the Lombardi group, as summarized in a recent review by Irving (2017), a model is proposed whereby a thick filament strain-dependent conversion of heads from an ordered state, perhaps the so-called super-relaxed state proposed by Cooke (McNamara et al., 2015), to a disordered state is associated with going from a myosin filament “off” state to an “on” state. Reconditi et al. (2017) have also proposed a “feed-forward mechanism” based on this concept that could possibly be used to explain myofilament length-dependent activation. A recent article (Ma et al., 2018) suggests, however, that thick filament–based regulation is likely to be more complex than the previous studies would suggest. For simplicity, in this article we include an increase of effective myosin binding rate to actin at longer SLs, where this increase in rate could be attributed to decreases in interfilament spacing, or any other mechanism associated with thick filament activation.
Toward resolving the underlying mechanisms of nemaline myopathy and myofilament length-dependent activation
Various mutations in nebulin and titin present in human disease strongly affect the contractile function in skeletal and cardiac muscles. Several hypotheses have been proposed to explain these observations, but clear mechanisms of the dynamic action of these proteins in the 3-D sarcomere lattice are still lacking. As shown here, MUSICO simulations can provide valuable insight into separating the effects of nonuniform length distributions of thin filaments, modulation of binding rates by the presence of nebulin, changes in interfilament lattice spacing, and the modulation of calcium sensitivity associated with thin filament calcium regulation.
The flexibility of the MUSICO platform and simulations of plausible mechanisms have provided good fits to the data, but many questions regarding the underlying mechanisms remain. For example, what could be the possible roles of nebulin in skeletal muscle or titin-based passive tension in cardiac muscle regarding modulation of calcium sensitivity? Is it enhancing Ca2+ binding to TnC, or reducing TnI affinity to actin, or could the force that a myosin molecule generates be enhanced? What are the molecular details of the mechanisms by which passive force alters calcium sensitivity and interfilament spacings? Are there additional mechanisms for modulation of calcium binding to TnC or TnI to actin? Does titin contribute to activating thick filaments and/or does it alter TnI–actin interactions?
Resolving the mechanisms in quantitative terms for the effects of variable thin filament lengths, passive tension, titin elasticity, and interfilament lattice spacing is complex and requires a thorough analysis with consistent sets of data that go beyond maximum isometric tension and force–pCa relations. To fit multiple experiments by the same set of parameters, we primarily used relatively well-established parameters. However, a few parameters that are less well established, for example myosin and TnI binding rates to actin, were adjusted to provide agreement with the experimental data. It is important to note that the chosen sets of parameters, although providing excellent agreement between MUSICO simulations and the observations (Figs. 6 and 7), may not be unique. Nevertheless, the parameters we used in simulations are within the range of values reported in the literature, including our own observations, and use of a different set of values is not expected to alter overall outcomes from those presented here.
To establish more consistent, and constrained, parameter sets, more experiments under similar if not identical conditions are required. For example, the history of force development, isotonic shortening, and fast transients needs to be included in determining model parameters. Overall good fits are not sufficient to quantitatively explain the involved mechanisms if they are not constrained with a sufficient number of experimental observations obtained under the same conditions. For another example, the analysis of calcium sensitivity is usually based on the normalized force–pCa data, but this analysis may, in some cases, be misleading because the number of strongly bound cross-bridges modulates position and states of regulatory proteins and therefore generated force. Thus, not taking into account absolute active tension levels can skew cooperativity effects, leading to interpretations that deviate from the underlying mechanisms.
The strong coupling between modeling and experiments in the present study, however, enabled quantitative separation of the contributions from structural changes in nebulin-deficient muscles, changes in myosin binding kinetics, and modulation of thin filament calcium regulation. This takes a significant step toward understanding the underlying mechanisms leading to muscle weakness in nemaline myopathy. Similarly, we quantified the effects of titin-based passive tension, length dependence of calcium activation, and interfilament lattice spacing on cardiac muscle’s ability to generate tension and alter calcium sensitivity. These are key steps toward understanding contractile performance at larger filling volumes of the heart. These findings suggest that the Frank–Starling mechanism of the heart may be due, in part, to an effect of titin-based passive force on the length dependence of maximal active tension and calcium sensitivity.
This study supports the notion that both nebulin and titin play important roles in muscle contraction. The effect of variable thin filament lengths in nebulin-deficient muscle accounts for about half of the decrease in isometric force and the leftward shift of the descending force–length limb observed in fully activated muscle. The match is further improved by 2.4- and 3.3-fold decreases in kbind at SLs of 2.5 and 2.0 µm, respectively. Thus, the changes in binding kinetics could be attributed to an enhancement of binding rates by the presence of nebulin, but the increase in calcium sensitivity at longer SL (∼2.5 µm) and the decrease in sensitivity at shorter (∼2.0 µm) SLs could be achieved by increasing or decreasing the sensitivity of TnI detachment from actin, respectively.
In cardiac muscle, isometric active force decreases when titin-based passive tension is reduced at both long and short SLs. Increases in interfilament lattice spacing decrease myosin binding rates at shortened muscle length, and overall tension decreases despite a significant increase in thick–thin filament overlap due to the nonuniform distribution of thin filament lengths. Sensitivity to calcium is decreased by both a decrease in passive tension and an increase in interfilament lattice spacing. Similar effects are observed in lattices compressed by 2.5% dextran, although isometric force is significantly higher due to the overall reduction in interfilament lattice spacing. In these cases, the changes in binding kinetics could be attributed only in part to enhancement of binding rates by reduced interfilament spacing at longer SLs (∼2.3 µm) or by compression by dextran. The latter effect could be linked to higher calcium sensitivity at higher passive tension at longer SLs. These overall effects can be attributed to myofilament length-dependent activation (Lakatta and Jewell, 1977; Kentish et al., 1986; Cazorla et al., 2001; de Tombe et al., 2010; Ait-Mou et al., 2016), but unveiling the detailed underlying mechanisms requires further investigation.
In summary, the MUSICO simulation platform is a multiscale model of the sarcomere that is able to reveal interdependencies that cannot be resolved experimentally. To realize its full potential to unveil fundamental underlying mechanisms leading to diseases such as nemaline myopathies in skeletal muscle and the subcellular basis of the Frank–Starling mechanism, more comprehensive datasets that tightly constrain the parameters are required. In this process, MUSICO, in itself, can provide a valuable hypothesis-generating tool to provide guidance on what experimental measurements would be needed to optimally constrain fitting the model to data. In so doing, MUSICO iteratively becomes an increasingly more effective tool for revealing fundamental contractile mechanisms in health and disease.
Geometrical aspects of interactions between myosin heads with binding sites on actin
Defining the myosin binding step in the 3-D sarcomere lattice is complicated, because myosin heads can interact with only one myosin binding site on actin out of several reachable adjacent sites (Fig. 3 A). The distance between a myosin head and the binding site on actin, , that determines the net cross-bridge force is calculated from the instantaneous spatial positions of each myosin molecule. Here the subscript m denotes the myosin head and l the associated binding site on actin reachable by head m. The vector is defined by four parameters: the axial strain , the radial spacing between centers of actin and myosin filament , and the relative angles and (Fig. 3). Note that for each myosin molecule m, there would be lmax binding sites on actin in the neighborhood of m. lmax-associated sites on actin can be located on one or two actin filaments, depending on the angle , thus can include multiple s, one for each associated actin filament (for details see Mijailovich et al. ).
In almost all sliding filament models, the strain-dependent rate constants between actomyosin states, denoted as (kij), exclusively depend on the strain component x, i.e., axial change of distance between a myosin head and the binding site on actin due to bond stretch Fig. 3 A, while all other components are ignored. In the stochastic 3-D sarcomere model used here, we adopt the same approach, but for possible binding of each myosin on multiple actin sites, we calculate strain-dependent rate constants , where is a cross-bridge strain of lmax-reachable actin sites. These constants are further modulated by weight factors that take into account the lattice spacing between the filaments, (Fig. 3 B) and the azimuthal angles α and β (Fig. 3 C). The resulting lmax binding rate constants per myosin molecule are then used for construction of binding probabilities for the stochastic process as explained below. To match the overall binding flux in probabilistic sliding filament models, the binding rate distribution is scaled down by the factor fsk that takes into account the average number of sites on actin reachable by each cross-bridge for the prescribed strain-dependent binding rate function (Mijailovich et al., 2016).
Strain-dependent rates in the actomyosin cycle
The three-state actomyosin cycle is defined by six state transition rates (Fig. A1), of which three are strain dependent. The strain dependence shown below is strictly defined in terms of strain in axial direction, x, whereas the effects of azimuthal movements are incorporated into the rate constants via the azimuthal weight factors.
For myosin binding to actin, the strain-dependent rate in quadratic form is derived from a Langevin type of equation balancing thermal fluctuations of the detached myosin molecule, elastic restoring, and inertial and viscous drag forces (Kramers, 1940; Papoulis, 1991; Hunt et al., 1994; Daniel et al., 1998). At zero strain (x = 0), binding rate is defined by overall rate kbind
as a function of the elasticity of the cross-bridges, κ, and the cross-bridge (axial) displacement from its unstrained position, denoted as x. The reverse reaction occurs at a constant rate , where is reduction in free energy due to myosin binding, kB is the Boltzmann constant, and T is absolute temperature in °K.
The transitions between the two attached states A·M·ADP·Pi and A·M·ADP are rapid, and this transition includes Pi release accompanied by a large (negative) change in chemical free energy and the displacement of the lever arm carrying out the power stroke, d:
For the strain-dependent ADP release (i.e., A.M.ADP → A.M), the forward rate is defined as function of the rate of ADP release when the elastic element is relaxed, , and the displacement that the lever arm must move to open the nucleotide pocket, δ:
The full expressions for the above strain-dependent rate functions are given in the Supplemental material.
Thin filament regulation by calcium
Structurally, the Tm molecule is a coiled coil of ∼40 nm in length (Censullo and Cheung, 1994), covers seven monomers on the same strand of the actin double-stranded helix, and is associated with one Tn forming a TmTn unit (Fig. 1). The Tn complex consists of three subcomponents: TnC, TnI, and TnT (White et al., 1987). In the absence of calcium, the N-terminal region of TnC is closed and the C-terminal of adjacent TnI is bound to actin, prohibiting Tm movement, i.e., holding Tm in the ‘‘blocked’’ state. Thus, in relaxed muscle, TnI holds Tm in an azimuthal position that sterically blocks myosin-S1 binding sites on F-actin (Mijailovich et al., 2012a).
In the presence of Ca2+, binding of one or two Ca2+ ions to TnC generates a conformational change in TnI, lowering the affinity of TnI to F-actin. The release of the TnI C terminus from F-actin allows the unconstrained Tm chain to move toward the ‘‘closed’’ state, i.e., the azimuthal position , favoring myosin binding to F-actin and, therefore, muscle contraction (Vibert et al., 1997; Smith and Geeves, 2003). The closed state permits weak myosin binding as observed at low Ca2+, but the strong myosin binding requires further movement of Tm chain to azimuthal position denoted as the “open” state.
There is structural evidence that neighboring tropomyosins overlap and that one end of TnT is bound to a specific site on Tm, whereas its N terminus overlaps the adjacent Tm, forming linked Tm–Tm regions. The interconnected neighboring TmTn units, therefore, form the appearance of a CFC (Lorenz et al., 1993; Vibert et al., 1997). There are two CFCs on each actin filament following the filament double-stranded helix structure.
CFC angular position and its variance
The Tm-actin spatial interactions are modeled as the interactions of loosely confined quasicontinuous semiflexible chain (CFC) with the actin surface in the presence of Tn. Each CFC spans one strand along the whole length of actin filament, having two CFCs per F-actin. The governing equation of the CFC interacting with F-actin is defined by the expression for the energy of a distorted chain with angular displacement at position s (Smith, 2001). For simplicity, the CFC is assumed to be elastically homogeneous, with bending stiffness, κTmTn, and an angular confinement potential well, with single minimum, αcp (Smith et al., 2003).
The angular displacements of the TmTn chain are dynamic, having a range of values caused by the energy of thermal fluctuations. The functional of thermally excited chain configurations is defined via the Feynman path integral (Feynman and Hibbs, 1965; Smith, 2001) over the complete set of confined worm-like configurations displacement . For 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 standard deviation of chain angles arising from thermal excitation (for details see Smith, 2001; Smith et al., 2003; Geeves et al., 2011; Mijailovich et al., 2012a).
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 standard deviation of fluctuating angle, σφ, has a value of approximately
where , where R is the radius at which Tm 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 the conditional constraints guarantee that the blocked, open, and closed states correspond to the three orientations seen in cryo-EM images (Pirani et al., 2005; Poole et al., 2006).
The configuration of a whole CFC, constrained at the positions of bound molecules of myosin and TnI along actin filament and their respective angular displacements , , and , is obtained as a piecewise semianalytic solution constructed, along each strand of F-actin, by merging the analytic expressions for the functions and σφ of neighboring TmTn units (Smith, 2001; Mijailovich et al., 2012a). The compatibility conditions require matching of angular displacements, slope, and curvature at common pinning points between neighboring chain segments, where each segment is defined by the functions and over the length of the arc satisfying the boundary angular displacements and slopes (Smith, 2001; Mijailovich et al., 2012a). The standard deviation at a pinning site is zero, and deviation increases for myosin binding sites further away from the pinning sites, and at distances beyond the confined persistence length, , the standard deviation approaches σo.
The coupling between the calcium-regulated position of the CFC and a three-state myosin cycle involves (a) 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 (b) restriction of TnI rebinding to actin by nearby bound myosins. The biochemical Tn-A states include a bound state of TnI to actin that maintains TmTn in the position inhibiting myosin-S1 binding, and the other state where TnI is not bound to actin allowing the TmTn 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 standard deviation, , where i denotes actin site at the discrete position, si, along an actin filament strand (Mijailovich et al., 2012a). 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
TnI and actin have two states, and the transition between these states is defined by the equilibrium state transition constant where kI is binding rate of TnI to actin. In the absence of calcium, most TnIs are bound to actin, but in the presence of calcium, calcium binds TnC and, in a [Ca2+]-dependent manner, increases KB (Mijailovich et al., 2012a) or decreases KI (Geeves et al., 2011). For simplicity, we assume that k–I is only dependent on [Ca2+], whereas the binding rate of TnI to actin is a weighted function of the rate of TnI binding to actin from the unweighted closed state rate, (Mijailovich et al., 2012a). 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 CFCs, i.e., two times the number of all actin filaments. 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 three-state model are constructed so that each state can transition to two neighboring states. The transition states from the detached state M.ADP.Pi includes two attachment probabilities, P12 and P13, associated with axial strain-dependent rates, including 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 Patt 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 k12. The scaling factor, fsk, decreases the magnitude of k12 at each x by an average number of reachable binding sites to myosin heads (Mijailovich et al., 2016). The weight factors Cα and Cβ 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 β (Mijailovich et al., 2016). The weight factors are proportional to the fraction of time step, Δt, when these transitions are possible, i.e., to the weighted probabilities that the CFC is at positions , which modulates weak myosin binding, denoted as subscript (Mijailovich et al., 2012a). The binding rate k13 is very slow and we assume the probability .
The transitions from the weakly to the strongly attached states (2 and 3) are defined by the probabilities and for the prestroke state 2 (A·M·ADP·Pi), and and for the poststroke state 3 (A·M·ADP). Calculation of these probabilities is almost identical as reported in Mijailovich et al. (2016), except P23 includes the weight factor that modulates weak to strong myosin binding, denoted as subscript . The factor is proportional to the fraction of time step, Δt, when the isomerization is possible, i.e., . (Mijailovich et al., 2012a).
We thank for graceful support the Mijailovich family, especially Dragica Mijailovic, Esq. (LLM).
This project was supported by National Institutes of Health grants R01 AR048776 (S.M. Mijailovich), P41 GM103622 (T.C. Irving), R35HL144998 (H. Granzier), and R01AR053897 (H. Granzier); by the British Heart Foundation 30200 (M.A. Geeves); and partially by H2020 European Research Council project 777204 – SILICOFCM (S.M. Mijailovich), Serbian Ministry of Science grants III41007 and OI174028 (B. Stojanovic), and National Institutes of Health R01 DC 011528 (R.J. Gilbert).
The authors declare no competing financial interests.
Author contributions: S.M. Mijailovich conceptualized the study, introduced new additions into MUSICO platform including thin filament regulation, implementation of new structural elements, and changes in sarcomere geometry, designed numerical procedures, analyzed the data from simulations, wrote the first draft and completed the final draft of the manuscript, and supervised each step of the work. B. Stojanovic assisted in designing numerical procedure, monitored implementation of changes in computer code, and supervised execution of the simulations, D. Nedic implemented changes in computer code and performed most of simulations. M. Svicevic assisted in running simulations and data analysis, M.A. Geeves contributed in development of the model of cross-bridge cycle and thin filament regulation. T.C. Irving and H.L. Granzier contributed in formulating the key points in discussion and assisted in writing the final draft of the manuscript.
Richard L. Moss served as guest editor.
This work is part of a special collection on myofilament function.