A spatially detailed myofilament model as a basis for large-scale biological simulations

Computational methods and models enable the modern biologist to investigate complex molecular mechanisms that are not easily accessible for analysis with experimental tools. The increased availability of low-cost computing power has enabled researchers to computationally model complex, emergent behaviors in biological systems. Such detailed models help researchers to bridge the gap between experimentally derived data and the underlying complex phenomena, which are not yet fully quantifiable. For example, the gross physiological response of striated muscle depends on the interactions of actin and myosin, which have yet to be fully characterized at the molecular level. Models continue to serve as quantitative tools to bridge the observable physiology with the underlying biophysical basis of complex muscle responses.

In this paper, we discuss our effort to model the actin–myosin interaction in cardiac muscle. In this system, complex behaviors, such as high sensitivity to activator Ca, are thought to arise from several cooperative interactions at the molecular level. (Ca stands for both calcium and the Ca²⁺ ion, a common convention in the field.) Here, cooperative interactions refer to coupling between proteins or protein complexes that modify behaviors that increase sensitivity to activator Ca. For computational efficiency, most current models are simplified and attempt to use lumped-parameter or mean-field approaches, two methods that are well known to researchers involved with mathematical modeling. For example, a single number that indicates concentration is used to represent large populations of proteins, and the model does not track individual proteins. Hence, these models lack spatially explicit representations of the cooperative interactions of component proteins. We choose the term “spatially explicit” to suggest that the spatial interactions between proteins are represented in the model; however, the full set of physical interactions between proteins and atoms are not computed, as is done in classical molecular dynamics simulations. In our modeling approach, which is classified as mesoscale or multiscale, higher levels of abstractions are used for molecular structures in order to speed calculations. For example, the model presented here can be executed with a time step of the order of microseconds instead of femtoseconds, which is a typical timescale for molecular dynamics simulations.

While simplified models are computationally efficient, such representations fail to account for basic responses such as the steep dependence of force on activator [Ca] and the proper temporal responses in muscle twitches [1]. (Hereafter, we use the symbol [Ca] to refer to the concentration of calcium.) In an effort to address these limitations of simplified models, we have developed a computational model that includes spatially explicit representations of cooperative interaction at the molecular level. Because biological modeling is currently an “art” as much as a routine procedure, we have attempted to choose a level of abstraction that balances biophysical detail with parsimonious representations that are computationally tractable. We demonstrate how simple spatially explicit models can produce complex Ca-activated responses that are similar to those measured in real muscle. Preliminary results of this modeling work have been discussed in [2].

Figure 1(a) is a diagram of a sarcomere, the basic repeating molecular-level structure of striated muscle (specifically skeletal and cardiac muscle). In fact, the striation in muscle arises from the high regularity of the sarcomere arrangement along the z-line direction. The end-to-end arrangement of sarcomeres produces regular and repeating striations along the length of the muscle. While the sarcomere is generally assumed to be the basic, repeating unit of striated muscle, we chose to model at the level of a half sarcomere. While our construction appears to exploit the geometric left–right symmetry, recent experimental characterizations have actually suggested some left–right asymmetries that can occur during contraction [3]. Hence, the half sarcomere may be a more logical functional unit than the full sarcomere. In the initial effort presented here, we model only a single pair of interacting thick and thin filaments. In reality, the sarcomere has a 3D structure [not visualized in the two-dimensional slice shown in Figure 1(a)], so that each thick filament is surrounded by six thin filaments that form a hexagonal lattice. Extension to larger, three-dimensional models will be a natural outgrowth of this modeling effort. Interested readers are referred to [4] for further details on how to incorporate multiple filaments to better capture the 3D properties of the half sarcomere.

Figure 1

In Figure 1(a), activated muscle generates force in a direction so as to pull the thin filaments to the center of the sarcomere. The magnified inset shows half of a thick filament (orange), composed of myosin that interacts with a single thin filament (green), composed primary of actin and regulatory units. Regulatory units (dark blue) can bind Ca, which in turn biases them to take a permissive conformation (shown schematically as raised locations) to allow myosin to bind. Myosin heads can bind to specific sites on actin (see gray head) and then rotate to stretch the extensible neck region (shown as a stretched string, or zigzag line, attached to the solid head in the diagram). Rotated heads generate a net force between thick and thin filaments.

Figure 1(b) shows details of the components of the thin and thick filaments. The thin filament is composed of a double helix of actin monomers. The regulatory proteins, troponin and tropomyosin, reside in the groove of the helix and serve to allosterically block interactions between actin and myosin. (The term allosteric refers to a change in shape and activity of a protein that can result from molecular binding with a regulatory substance.) The important feature of troponin is that it can bind one Ca ion at a single binding site that controls the allosteric switching behavior. As shown, troponin is actually composed of three subunits: TnC, which binds Ca; TnT, which serves as a linker to tropomyosin; and TnI, which has a modulatory role in Ca-based regulation. Tropomyosin spans seven actin monomers and blocks the myosin binding sites on actin. Through the interactions with troponin, mainly via TnT, the tropomyosin also makes an allosteric shift in response to changes in the Ca level. As shown, tropomyosin molecules overlap in an end-to-end fashion and are thought to communicate with their neighbors via this physical communication (for a review, see [5]).

Myosin has three major structural subunits—the head, neck, and tail. The head attaches to a binding site and rotates using the energy from the conversion of adenosine triphosphate (ATP) to adenosine diphosphate (ADP), a common energy-liberating reaction in cells. The bound linkages between actin and myosin are commonly termed crossbridges to reflect bridging linkages between thick and thin filaments as first revealed in pioneering electron microscopy studies [6]. Rotating the head can stretch the extensible neck region as shown schematically in Figure 1(a). The tail regions of myosin assemble together to form the thick filament as shown in Figure 1(b).

Thomas Daniel and colleagues [7] have developed a model of two filaments in a half sarcomere that includes the spatial interactions of myosin and actin binding sites. The main goal of their work was to study the effects of compliance (i.e., extension or displacement of a loaded structure per unit load) on crossbridge formation and force development. We have used this approach to form the spatial layout of myosin and actin binding sites with the appropriate compliances between the elements. We assume that the thick filament has myosin heads with appropriate orientations at an intrinsic spacing of 43 nm, as shown in Figure 2(a). On real myosin, the heads extend in a helical fashion such that only a subset will appropriately align to interact with the single thin filament assumed in this model. Note, however, that the 43-nm intrinsic spacing of myosin is slightly larger than the 37-nm spacing of appropriate binding sites on the thin filament, as given in [7]. Similarly to the case with myosin, the helical nature of the actin will restrict binding to a subset of actin monomers that appropriately face the thick filament in the two-filament model presented here.

Figure 2

In our model, the thick filament, thin filament, and crossbridges are treated as a system of springs that operate linearly around the rest lengths of the springs. The spring constant for the thick filament (k_m) is 2,020 pN/nm (piconewtons per nanometer), that of the thin filament (k_a) is 1,743 pN/nm, and that of the crossbridges (k_xb) is 1 pN/nm [7]. (Hereafter, “xb” is used as an abbreviation for crossbridge.) The rest lengths are set by the intrinsic spacing of myosin sites (m_s = 43 nm) and actin binding sites (a_s = 37 nm). The forces generated by the crossbridges are not computed using rest length, despite the fact that the rotational state of the head can modify the effective length by 7 nm. (This effect is further explained shortly.) We also assume that no external forces act on this system, and that other effects, such as viscosity, are negligible. Figure 2(b) shows a segment of this system of linear springs.

The new positions of the thick filament, thin filament, and crossbridges can be determined using a pair of force-balance equations for each crossbridge. Equation (1) is a coupled set of linear equations for the interaction of myosin head i and binding site j:

The set of equations (1) for all n myosin heads and p binding sites can be transformed into a position vector X = [x₁, x₂, …, x_n, y₁, y₂, …, y_p], a (p + n) × (p + n) matrix K of spring constants, and a vector A, described shortly. For the model implemented here, we assume that the thin filament is 1.1 μm long and half of the thick filament is 0.825 μm long. These values are similar to experimental characterizations, although some controversy exists as to exact lengths in cardiac muscle [1]. The lengths given correspond to 32 actin binding sites (p = 32) and a total of 19 myosin heads (n = 19).

The matrices and vectors X, K, and A are related as shown in Equation (2):

The products of the spring constants and the rest lengths of the crossbridges define the A vector. Changing the value of a component of the A vector is equivalent to an external force acting on the corresponding myosin head, actin binding site, and possibly a crossbridge, depending on the configuration of the K matrix. In our realization of the model, for each attached crossbridge in the PreF¹ state, a force of k_xb × 7 nm = 7 pN is added to or subtracted from the component of the respective myosin or actin binding site in the A vector. This change of forces represents the force correction for an extended crossbridge head before rotation to the F state, which reduces apparent strain by 7 nm. However, since this correction force is not added to the components when the crossbridge is in the F state, the transition of the crossbridge to the post-rotation state affects the developed force of the system.

Equation (3) is used to determine the new positions of the crossbridges, myosin heads, and binding sites during each step as the K matrix and rest lengths vary on the basis of the state and interactions of the myosin heads and binding sites:

Changes in the interactions among crossbridges, myosin heads, and binding sites result in changes to the appropriate elements of the K matrix. Specifically, an attached crossbridge is incorporated as extra spring-constant (k_xb) terms in the K matrix. If the myosin head i is bound with binding site j, −k_xb is added at K[i][j] and +k_xb at K[j][i]. If the crossbridge detaches, these terms are removed.

For all simulations presented here, we assume a fixed sarcomere length. This assumption is implemented by fixing the positions of the first myosin (x₁) and last actin (y_p) sites, and one can envision the thick and thin filaments as pulling between two immovable points. The developed tension of the ensemble can be most easily calculated from the developed tension at either of the non-overlapping ends of the filaments (the first myosin site or the last actin site). In our implementation of the model, we determined the value of the strain experienced by the segment between the locations of the first and second myosin heads. Hence, force for the whole ensemble is calculated as F_total = k_m(x₂ − x₁ − m_s).

The state transition rates for crossbridges are derived as a function of the strain, , which is a measure of the distance between the actin binding site and a reference point on the thick filament for a given myosin. As in the classic modeling efforts (such as those in [8]), the reference point is also known as the equilibrium point of an attached crossbridge, so that k_xbΔx = 0 at this point. Thus, an attached crossbridge can generate positive or negative force depending on its direction from the equilibrium position. In particular, positive strains produce positive forces in the direction of muscle shortening, whereas negative strains produce negative forces that inhibit muscle shortening.

Our approach to develop strain-dependent crossbridge cycling rates is inspired by the work of Pate and Cooke [9], who define the energy profiles for each crossbridge state as a function of its relative position. Pate and Cooke assumed a five-state crossbridge model that included the role of metabolite concentrations (ATP, ADP, and inorganic phosphate). We consider a reduced three-state model without metabolite dependencies to be a reasonable compromise between parsimony and biophysical detail. Daniel et al. propose a similar three-state crossbridge scheme [7], but differences exist in the state definition and transition rates. More detailed crossbridge cycling schemes with more states can be incorporated by using analogous methodology in future versions of the model.

The model defines the energy profile of each crossbridge state as a function of its relative position, and equations are defined to govern the transitions between the crossbridge states. Equation (4) provides the ratio of reversible transition between two states, and the Gibbs free energy for the states:

(4)

where R_{i_j}() is the transition rate from state i to state j at strain , G_i() is the Gibbs free energy of state i at strain , R (with no subscript) is the gas constant, and T is the temperature.

Figure 3 shows the energy profiles and transition rates as functions of distortion. Note that Gibbs free-energy profiles define only the ratio of rates, not the absolute values. Therefore, the absolute values of rates are adjusted so that transitions between states will produce realistic net ATP consumption rates. Given that the three-state scheme is simplistic, the strain-dependent transition rates are reasonable but approximate guesses. Full validation of rates under all conditions would entail a research paper in itself, and hence is beyond the scope of this preliminary report.

Figure 3

Note that upper and lower limits on transition rates are included to prevent extremely large variation in magnitudes. Such a large variation in the rates produces problems in Monte Carlo (MC) simulations because required time increments become extremely small (a situation analogous to the use of stiff equations in other numerical methods). Also, large variations in rates require state transitions in MC computations that are based on evaluation of low-order bits from pseudorandom number generators. In typical pseudorandom number generators, such low-order bits are less random than high-order bits, and hence the simulation results become unreliable. The complete mathematical formulation of free-energy levels and the state transition equations are given in the Appendix.

A critical feature of Ca-based activation is a steep force–Ca (F–Ca) relationship. As shown in Figure 4(a), real muscle shows high sensitivity to Ca, as quantified by a Hill coefficient in the range of 7–10. Figure 4(b) illustrates the physiological effects of the steep Ca sensitivity. During each heartbeat, the intracellular [Ca] changes by about an order of magnitude. In contrast, the resulting force can change by three or more orders of magnitude. The presumed rationale for such a biological design is that large changes in developed force are required for effective pumping of blood. Incomplete relaxation leaves a residual force that inhibits proper filling during diastole, which in turn leads to the ejection of only a small volume of blood during systole. In contrast, cardiac cells generate relatively small changes in [Ca] because substantial energy in the form of ATP is required to actively raise and lower [Ca] on each heartbeat. As a result, the steep nonlinearity of the F–Ca relationship bridges the disparity between the small change in [Ca] and the large changes in force.

Figure 4

To understand the model of Ca-based activation, we consider a functional unit of seven actin monomers spanned by one troponin/tropomyosin (T/T) unit and with one associated myosin that can form a crossbridge. As seen in Figure 1(b), the actual arrangement of linked tropomyosin proteins follows the two opposite grooves of the actin double helix. Hence, researchers should model two T/T units per seven actin monomers with two separate one-dimensional (1D) lattices of regulatory units per thin filament. We have not included this level of spatial detail, but instead we assume a single 1D lattice of regulatory units that controls all actin binding sites. A more precise spatial arrangement of regulatory units and actin binding site is planned for future versions of the model with multiple thick and thin filaments.

Figure 5 shows a diagram of the states in Ca-based activation. With respect to Ca activation only, a functional unit can be in one of four states (0N, 1N, 0P, 1P). As noted in the figure legend, the states with a 0 prefix indicate the state in which a Ca ion is not attached to the binding site, and the states with a 1 prefix indicate the attachment of a Ca ion to the binding site. The N states indicate that the binding site is in a nonpermissive conformation, so that no crossbridges can attach to that site. The P states indicate that the binding site is in a permissive state and myosin heads can attach. For the moment, the explicit crossbridge states are not considered; instead, only Ca activation steps are included. The forward Ca-binding rate is assumed to be diffusion-limited and independent of the permissive state of the T/T unit. Hence, the forward rate constant for a nonpermissive T/T unit (k_on) is equal to the forward rate constant for a permissive T/T unit (k’_on). The reverse rate is assumed to depend on the permissive state of the T/T unit. The model assumes k_off > k’_off which is in agreement with experimental evidence that the affinity of troponin for Ca increases when activated in the presence of cycling crossbridges [10]. The relationship between k_off and k’_off is set by a parameter μ > 1 such that Equations (5) and (6) hold:

(5)

(6)

(The parameter μ is related to binding affinities for T/T units.) For all runs in this paper, we set μ = 15, which is a plausible value based on several experimental estimates [11].

Figure 5

An important cooperative mechanism thought to exist in cardiac muscle involves nearest-neighbor cooperativity between adjacent T/T units. Several lines of experimental evidence suggest that cooperative behavior results from the overlap of several amino acids allowing communication from tropomyosin to tropomyosin [5, 12]. The model implements such a cooperative mechanism by making the rates between nonpermissive and permissive states depend on the states of the two nearest-neighbor T/T units. The nearest-neighbor interactions are set by a parameter  ≥ 1, whose physical interpretation is an energetic penalty for neighbors being in different permissive conformations. The neighbor dependencies appear as exponents on and are computed according to the number of neighboring units in the permissive conformation, as shown in Figure 5(b). The exponent n can take on the values of 0 for no permissive neighbors, 1 for a single permissive neighbor, and 2 for both neighbors permissive. The net effect of the ⁿ terms is to increase the nonpermissive-to-permissive transition rates when the neighbors are also permissive. Similarly, the ⁻ⁿ terms decrease the reverse rates from permissive to nonpermissive states. Therefore, an individual T/T unit is more likely to make the transition to permissive when its neighbors are permissive. Similarly, an individual T/T unit is more likely to make the transition to nonpermissive when its neighbors are nonpermissive. Hence, the ⁿ and ⁻ⁿ terms promote uniformity along the thin filament, so that T/T units tend to assume the permissive states of their neighbors. Interested readers are referred to [1] for a detailed discussion.

As shown in Figure 1, the half-sarcomere contains regions in which thick filaments overlap thin filaments; also, there are regions in which thin filaments do not overlap the thick filaments. While one can envision a case in which thick filaments do not overlap the thin filaments, such a case does not typically occur in cardiac muscle because the extension of the sarcomere is physically limited. To represent the overlap of thick and thin filaments, the model represents two cases:

• Case 1: If the thin filament does not overlap the thick filament, no crossbridges can form. In this case, each binding site on the thin filament is governed by the four-state model shown in Figure 5. Specifically, while T/T units can bind Ca and switch between nonpermissive and permissive states, no myosin can bind.

• Case 2: If the thin filament does overlap the thick filament, crossbridges can form according to the crossbridge kinetic scheme in Figure 3(a). However, a slight complication exists because myosin cannot bind until the T/T unit is in a permissive state. Once it is in a permissive state, a nearby myosin can transition from the P state to the PreF state (see Figure 6). The P state corresponds to a detached or weakly bound crossbridge and a permissive T/T unit as intermediate states. Once in states PreF or F, the myosin is strongly bound and prevents the T/T unit from returning to a nonpermissive conformation. This construct derives from the allosteric interactions thought to occur in such a way that an attached myosin head physically prevents T/T from shifting back into the groove on actin. Moreover, experimental evidence shows that the presence of one or two crossbridges can prevent the thin filament from switching back from a permissive conformation to a nonpermissive conformation [13, 14].

Figure 6

To account for the interactions between T/T units on binding sites and crossbridge formation, a larger coupled kinetic system must be employed, as shown in Figure 6. In this model, a binding site can be in any of the eight states, the associated crossbridge can be in any of the four states (0Pref, 1Pref, 0F, 1F), and the myosin head can be in any of the six states (0P, 1P, 0Pref, 1Pref, 0F, 1F). The states with a 0 prefix indicate the state in which a Ca ion is not attached to the binding site, and the states with a 1 prefix indicate the attachment of a Ca ion to the binding site. Once in a permissive state, a nearby myosin can attach. Initially, the myosin is assumed to be in a detached P state, which is mapped to 0P or 1P in the combined eight-state model. Note that Ca is assumed to be able to bind or detach independently of the permissive conformation or myosin attachment events, so these states have mirror Ca-unbound (0x) and Ca-bound (1x) versions. For example, the 0PreF state corresponds to a strongly bound myosin with a pre-rotation head region with Ca unbound, and the 1F state is a strongly bound myosin with a post-rotation head region with Ca bound.

We used a typical Monte Carlo (MC) method to simulate the state transitions of binding sites, myosin heads, and crossbridges. In such an MC method, time is divided into time steps of Δt_i, where i denotes the ith time step (the length of time steps may vary, as described below). The simulation starts with all crossbridges in the detached state and all binding sites in the 0N state. The state transition simulation is carried out as two steps. During the first step, the state transitions of the T/T units of binding sites are determined as follows. For each binding-site state, there are two possible transitions to new states, denoted B1 and B2, with respective rates Br₁ and Br₂. For example, in Figure 5(a), if the initial state is 1P, the possible transitions are to 0P and 1N with respective rates k’_off and ⁻ⁿk_{pn_1}. For small values of Δt_i, the probability of transition to states B1 and B2 is approximated by Δt_iBr₁ and Δt_iBr₂, respectively. A pseudorandom number (x_random1) in the [0, 1] range is generated and is compared with the transition probabilities Δt_iBr₁ and Δt_iBr₂ for that state. If x_random1 is less than Δt_iBr₁, the binding site transits to state B1. If the first case is not satisfied, and if x_random1 is less than (Δt_iBr₁ + Δt_iBr₂), the binding site transits to state B2. If x_random1 is greater than (Δt_iBr₁ + Δt_iBr₂), the binding site continues in the current state for the time step.

Once the state transitions for the binding sites are completed, the crossbridge states are similarly updated. A crossbridge may be in one of the three states—detached, pre-rotation, and post-rotation [see Figure 3(a)]. For each crossbridge state, there are two possible transitions to new states, denoted C1 and C2, with respective rates Cr₁ and Cr₂. A second pseudorandom number (x_random2) in the [0, 1] range is generated and is compared with the transition probabilities Δt_iCr₁ and Δt_iCr₂ for that state. If x_random2 is less than Δt_iCr₁, the crossbridge transits to state C1. If x_random2 is less than (Δt_iCr₁ + Δt_iCr₂), the crossbridge transits to state C2. If the first case is not satisfied, and if x_random2 is greater than (Δt_iCr₁ + Δt_iCr₂), the crossbridge continues in the current state for the time step.

Note that the coupled system in Figure 6 could be solved using a single call to a random number generator. We instead chose to update the actin binding-site states and the crossbridge states independently in order to simplify implementation. However, systems are coupled by appropriately changing transition rates to reflect the interaction of T/T units and crossbridges in Case 2 discussed above. Specifically, as we mentioned, myosin cannot bind until the T/T unit is in a permissive state (see Figure 6). Hence, for a crossbridge near a nonpermissive T/T unit, the rate R_P→PreF() equals 0 [see Figure 3(a)]. Once a crossbridge is in state PreF or F, the myosin is strongly bound and prevents the T/T unit from returning to a nonpermissive conformation. This feature is incorporated by effectively setting the transition rates from permissive to nonpermissive states to zero [k_{pn_0} = k_{pn_1} = 0, Figure 5(a)] for a T/T unit with a strongly bound crossbridge. As one final note, the longest possible time step is used to decrease the computation time. The largest appropriate time step can be computed as Δt_i = (3 × r_max)⁻¹, where r_max is the global maximum rate found after computing the transition rates for each T/T unit and crossbridge in the model at time step i.

In our model, the high sensitivity of cardiac muscle to activator Ca can be attributed to the properties of the regulatory proteins. The T/T units transition from nonpermissive to permissive states in response to Ca. After becoming permissive, myosin can attach to form crossbridges and produce force; thus, in a simple scheme, the regulatory proteins act as switches to modulate the amount of force. While this is not quite so simple in the complete system, we can investigate the Ca sensitivity by first considering only the switching behavior of the regulatory proteins. To demonstrate Ca sensitivity, the number of permissive T/T units is plotted as a function of the activator [Ca] in Figure 7. Here, crossbridges are prevented from attaching [R_P→PreF() = 0] in Figure 3(a) to isolate the Ca-activation effects of the regulatory units. Figure 7 shows a parameter variation with , the energetic penalty assigned to neighboring T/T units with different permissive conformations. Hence, greater values correspond to greater strength of the nearest-neighbor cooperativity. Using minimization of squared error, true Hill functions (not shown in the figure) are fit to the simulation data. The best fits are found with the following values:  = 1, Ca₅₀ = 2.1 μM, N_H = 1.1;  = 5, Ca₅₀ = 2.5 μM, N_H = 2.5;  = 20, Ca₅₀ = 4.4 μM, N_H = 6.0; and  = 45, Ca₅₀ = 4.7 μM, N_H = 8.1.

Figure 7

As shown in Figure 7, increased nearest-neighbor cooperativity leads to steeper transitions, beginning with most T/T unit transitions being nonpermissive and transitioning to most units being permissive as [Ca] increases. Specifically, when Hill functions (defined in the legend for Figure 4) are fit to the number of permissive T/T units, the steepness, as quantified by Hill coefficients, increases as a function of the value. For example, the  = 1 case corresponds to no nearest-neighbor cooperativity. This case produces a Hill coefficient very close to 1, as could be predicted by the single regulatory Ca binding site on each troponin unit. However, as increases to 5, 20, and 45, the apparent cooperativity increases, as respectively quantified by the Hill coefficients of 2.5, 6.0, and 8.1.

Similar results were shown in a previous study [11] using an Ising approach that produces an analytic result for equilibrium conditions such as those shown in Figure 7. Note that analytical results show minor but systematic deviations from true Hill functions and that the deviations are similar to those found in real muscle (see [11] for details). Such deviations are much harder to discern in MC methods, in which stochastic fluctuations can easily obscure subtle features of muscle responses. This becomes more apparent in the next section, in which cycling crossbridges are added. While the MC approach in the current study is considerably more noisy and computationally costly than the Ising approach, the subsequent data show nonequilibrium conditions that cannot be calculated using the Ising approach.

After becoming T/T-permissive, myosin can attach to form crossbridges and generate force. Hence, to a large extent, the amount of developed force should reflect the upstream switching events of the regulatory proteins; however, important differences can also exist. Unlike the data in Figure 7, the data in Figure 8 are generated with a normal attachment rate for myosin, so that crossbridges can form and generate force. Developed force can be computed from the attachment of crossbridges and the corresponding strains and can be compared to experimental data such as that shown in Figure 4(a); however, the addition of attaching crossbridges produces more fluctuation than the corresponding plot in Figure 7 with  = 45. The simulation in Figure 8 is run with  = 45 and μ = 15. The Hill fit shown has Ca₅₀ = 4.2 μM, N_H = 10.2.

Figure 8

Figure 8(a) shows the number of binding sites in permissive states as a function of activator [Ca], similar to the plots in Figure 7. Note, however, that the addition of attaching crossbridges produces more fluctuation than is seen in the corresponding trace in Figure 7 with  = 45. The Hill fit shown corresponds to Ca₅₀ = 4.2 μM, N_H = 10.2. The Hill equation is often used by biochemists to describe the fraction of the enzyme saturated by a ligand as a function of the ligand concentration. Here, the Hill equation has two parameters: N_H, which is known as the Hill coefficient, and Ca₅₀, which refers to the [Ca] that produces 50% of maximum activation (see caption for Figure 4). A larger value of N_H indicates a larger degree of cooperativity and produces a steeper sigmoidal curve. A smaller value of Ca₅₀ indicates that a smaller [Ca] produces half activation and therefore indicates a greater overall sensitivity to activator Ca.

As shown in Figure 8(b), developed force exhibits much more variation than the corresponding trace for the number of permissive T/T units. The increased variation is an obvious byproduct of the smaller number of attached crossbridges compared with a total of 32 T/T units on the thin filament. Moreover, the attached crossbridges generate force primarily in the F state, a feature that further reduces the effective number of force-generating events and increases the level of fluctuation. A Hill function fit has [Ca₅₀] = 4.0 μM and N_H = 11. Interestingly, the plots in Figure 8 show somewhat higher apparent cooperativity than is seen in the corresponding trace in Figure 7 with  = 45. This larger N_H value is an artifact of the noisy data, the coarse increment in [Ca], and the fitting procedure. Rerunning with a finer increment in [Ca] over a concentration range around [Ca₅₀] (data not shown) produces estimates of [Ca₅₀] = 3.7 μM and N_H = 8.7, which correspond more closely to the results in Figure 7.

The F–Ca relationship is convenient for demonstrating Ca sensitivity. However, fixing [Ca] at a constant level is an experimental manipulation that does not correspond to physiological conditions in the actual heart, in which [Ca] is constantly changing throughout each heartbeat. To better match these dynamic conditions, we simulated a twitch response for a transient increase in [Ca] similar to that which occurs during a typical heartbeat. Figure 9(a) shows a simulated twitch activated by simulated Ca transients that are generated to resemble those from real muscle. The use of simulated Ca transients is typical in modeling papers (e.g., [15]) because experimentally measured transients contain artifacts from a nonlinear response of the Ca-sensing dyes. Also, developed force should not depend too closely on the exact Ca transient shape. In Figure 9(b), similar data are shown from the study of Backx et al. [16]. Note that in Figure 9, the [Ca] has peaked and is falling when the peak of the developed force occurs. The simulated twitches in part (a) have activator [Ca] that is higher by about a factor of 10 than in the corresponding experimental data in part (b). This discrepancy occurs because the model results shown in this figure correspond to skinned fiber data in Figure 4(a) that produce a lower Ca sensitivity than observed in intact cell data. (The term skinning refers to removing the cell membrane; see [17] for more details on this effect.)

Figure 9

The twitch responses illustrate some important features of the system. Specifically, the Ca transient is relatively fast compared to the force transient, which shows slower response times. One common interpretation is that relatively slow crossbridge dynamics predominate to produce the slow force response (e.g., [18, 19]). These preliminary results also show some other complex behaviors. For example, the final relaxation of force is slower for larger peak force than for smaller peak force, as has been seen experimentally [18, 20]. Such an effect is greatly diminished for lower levels of nearest-neighbor cooperativity (e.g., near 1, data not shown). However, the comparison is not completely straightforward. Specifically, decreased strongly affects Ca sensitivity, so that developed twitch force is not as strongly modulated by the different levels of activator [Ca] shown in Figure 9. While one could change the levels of activator [Ca] to more closely match the developed force shown in Figure 9, the results can still be difficult to interpret, because the relaxation phase of twitches involves a complex interplay of Ca-activation and crossbridge kinetics (e.g., see [17]).

One interpretation of the prolonged twitches is that higher levels of force entail the binding of a greater number of crossbridges, which in turn can hold the thin filament in permissive conformation so that neighboring crossbridges are less likely to detach. While the force is generated by a relatively few crossbridges that may be separated by long distances, the nearest-neighbor interactions along the thin filament will extend the influence. Recall that in the model construction, a strongly bound crossbridge can hold the adjacent T/T unit in a permissive conformation even after Ca has dissociated. From previous modeling work using similar nearest-neighbor cooperativity ( = 45, μ = 15), correlations between T/T units are seen to spread up to 13 units along the thin filament [11]. Eventually, stochastic variation causes enough crossbridges to detach so that the thin filament can become nonpermissive, and complete relaxation ensues. In other words, generated force is somewhat self-sustaining, so that larger force transients are prolonged and relax more slowly compared with lower force transients that relax quickly.

The slow response of force is demonstrated in phase loops in which the relationship of force and activator [Ca] is plotted during the twitch. As shown in Figure 10, the resulting plots are loops that are traversed in a counterclockwise direction during the time course of the twitch. Also shown are the steady-state F–Ca relations that can be plotted on the same axes. The dashed trace is the Hill function fit to the steady-state F–Ca data with [Ca₅₀] = 3.7 μM and N_H = 8.7. Early in the twitch, [Ca] has peaked while force lags, placing the phase loop below the steady-state F–Ca relation. In contrast, later in the twitch, force is still large while [Ca] is decreasing, placing the phase loop above the steady-state F–Ca relation. These results demonstrate that the complete model has dynamic properties that cannot be predicted from the steady-state behavior. Hence, the F–Ca relation is a reasonable characterization of Ca sensitivity but is an incomplete description for dynamic responses in which slower crossbridge kinetics also play a role in shaping temporal responses.

Figure 10

Typical runs to compute steady-state force require 160 s of simulation time (approximately 1.6 × 10⁹ time steps) per [Ca] level in Figure 8. The programs are written using C++ and are executed on an IBM POWER5* platform. Each 160 seconds of simulation time requires about 2.4 × 10⁴ seconds in real time to complete. As shown in Figure 11(a), the standard deviation is generally larger as the mean force increases. Such a result is perhaps not too surprising, given the many sources of stochastic variation in crossbridge binding and myosin head rotation from the PreF to the F state. Hence, under the default-model parameter considered in this paper, generated force results from a relatively small number of stochastically controlled crossbridges, and hence variability is quite large. A small number of force events result partly from the different spacing of actin and myosin sites that produce only a small number of crossbridge bindings. This effect is an artifact that is explored in the next section.

Figure 11

As shown in Figure 2, the intrinsic spacing of actin and myosin positions is different. The difference in spacing results in certain pairs of myosin and actin binding sites having higher probabilities of binding than other pairs. This phenomenon is explored in Figure 11(b), in which the probability of an attached crossbridge is shown for every actin binding site on the thin filament (numbered 0 to 31 on the abscissa). The probabilities for having a strongly bound crossbridge are computed over time for each actin site [strongly bound crossbridges correspond to states PreF and F in Figure 3(a)]. Note that some binding sites have large probabilities for strongly bound crossbridges, while other binding sites have near-zero probabilities. For example, at the default length of 2.2 μm, peaks exist at actin sites 5, 11, and 18–19.

Figure 11(b) shows a total of five runs in which the sarcomere length is increased by 8 nm for each subsequent run. Note that the peaks in probability change for small variations in length. This effect suggests that reporting results at a given length may produce artifacts in that the system has preferred myosin and binding-site pairs. These pairs reflect the near-perfect spacing and alignment in the simulation, which is not likely to occur in a real muscle. Specifically, real muscle is more likely to show some variation in binding-site spacing and sarcomere lengths, so that peaks in probability are likely to be smeared along the abscissa. Moreover, in real muscle, each thin filament is surrounded by three thick filaments, each of which has a different phase in the actin and myosin site spacing. Therefore, one can reasonably predict that the precise peaks in the probability are artifacts of the artificial simulation geometry and do not reflect the situation observed in a real system. This limitation is addressed further in the discussion.

The model described here represents an attempt to formulate a mathematical representation of a thick and thin filament in the sarcomere. The model builds on a number of previous modeling efforts, but seeks a level of spatially explicit detail that was generally considered intractable for most earlier studies. We propose that two classes of advancements have made such detailed modeling efforts feasible. The most obvious advancement is the decreased cost of computation, even on the terascale level, in which teraflops of computation power can be obtained in a single rack of hardware. Previously, terascale computation was obtainable only at the largest supercomputer facilities. The second advancement involves the theoretical underpinning of the models. Myofilament modeling continues to be an active area of research, with new or refined models still being proposed. While the models exist at very different levels of abstraction and spatial detail, taken as a whole, the level of understanding of the biophysical mechanisms has increased. We suggest that the construction and validation of spatially explicit and biophysically detailed myofilament models, when coupled with newly available experimental characterizations, has become feasible.

A primary goal of the modeling described above is to serve as a basis for higher-level models of full sarcomeres or groups of sarcomeres. The preliminary results presented here are promising. The model produces F–Ca relations and twitches that resemble those from experimental characterizations. Moreover, many previous models fail to recapitulate even these basic responses, and the failure has been attributed to a lack of explicit consideration of cooperative mechanisms [11].

Despite some success, we must also consider important limitations; we can broadly define several such limitations. First, much work remains to be done in refining and validating the current formulation. For example, the three-state crossbridge kinetic scheme is simplistic, with rough guesses for the strain-dependent transition rates. Further studies will refine the rates and could possibly incorporate additional states as well as metabolite concentrations that also strongly affect rates [9]. Also, the data shown in the current work have all been collected under isometric conditions, with a single fixed sarcomere length (that is, all data shown were collected at sarcomere lengths equal or very close to 2.2 μm). Further studies must address variations to different fixed sarcomere lengths, which are known to have a strong effect on Ca sensitivity and maximal force in cardiac muscle [21]. Moreover, allowing active shortening to produce a net velocity between the thick and thin filaments will involve another level of complexity. For example, increased shortening velocity will produce constantly changing strains that will affect myosin transition rates as well as the force generated by attached crossbridges. We do not currently know of any models that adequately recapitulate the range of physiological responses to the varying Ca activation, velocity, and load conditions seen in muscle. In general, the existing models have a smaller scope and can be applied to limited sets of conditions. Given that the model in this paper is built by combining several models of more limited scope, one should realistically expect some difficulties when attempting to model very diverse phenomena.

The second broad set of limitations concerns spatial issues. The model presented here makes use of a single pair of thick and thin filaments. We have not attempted to modify the filament compliances (Figure 2) in our model presented here. However, previous results with a similar two-filament model by another group have suggested that changing filament compliance can strongly affect force [7]. In this previous study, attached crossbridges produce compliant realignment of myosin and actin binding sites, which can cause additional attachment events and force. In general, a more compliant thin filament can increase such realignment and developed force. However, a compliant thin filament will also transmit less force from the myosin head rotation events.

The use of a single pair of thick and thin filaments has additional limitations. In mammalian muscle, the cross section reveals a hexagonal lattice of thin filaments surrounding each thick filament. In particular, each thick filament is surrounded by six thin filaments, while each thin filament is surrounded by three thick filaments. This spatial arrangement typically increases the number of interactions between all pairs of adjacent filaments. In addition, the effects of preferred pairs of binding sites and crossbridges, as seen in Figure 11(b), should be decreased in models that go beyond the use of single pairs of thick and thin filaments. The helical nature of thick and thin filaments ensures that each thick filament will have slightly different alignment and register with each of its six adjacent thin filaments.

The larger number of interacting filaments in more complex models also increases the spatial averaging of force. As shown in Figure 11(a), the single pair of filaments generates a noisy response, even if time-averaged over 120 s. In contrast, a steady-state response in real muscle is obtained very quickly, because a high level of spatial averaging exists for crossbridge interactions spread over the whole lattice of thick and thin filaments. We can reasonably hope that subsequent models with a multitude of thick and thin filaments in the appropriate lattice arrangement will show both decreased effects of preferred binding sites [Figure 11(b)] and lower levels of noise [Figure 11(a)]. This leads us to the next section, which discusses how such models might be constructed to run on terascale supercomputer platforms.

A primary goal of the modeling we describe is to serve as a basis for higher-level models of full sarcomeres or multiple sarcomere structures. While the availability of relatively inexpensive computational power would seem to make such models possible, the construction of full working models is more difficult than simply replicating a single pair of thick and thin filaments. Specifically, while the current model could be replicated and simplistically distributed to a large cluster of independent workstations, the net result would only be to speed the collection of sufficient data to do averaging of the stochastic MC results. A more valuable simulation requires us to build communications between adjacent thick and thin filaments to generate a full sarcomere model. Moreover, if multiple sarcomeres can be assembled, one can build a myofibril, which is a small unit of muscle that is currently used for fundamental experimental characterization [21].

A myofibril can be constructed from the half-sarcomere model as follows. A 3D lattice of sarcomeres is created so that the lattice models the real muscle geometry. An activated Ca signal is provided at each lattice space at each time step. The activating signal may range from a simple, homogeneous waveform to a more complex signal such as a propagating 3D Ca wave that could more closely resemble a physiological condition in the cell. Each lattice point is then individually simulated for the Ca level at its position, and the generated force is determined. The shortening velocity for each sarcomere is next determined on the basis of the force generated in its neighborhood. This velocity is communicated to individual sarcomere units, which reconfigure accordingly and proceed to the next time step. Such a system may permit the investigation of inhomogeneous relaxation responses, as has been recently characterized at the level of myofibrils in real muscle [3, 22].

As an example of how such a system could be implemented on a terascale-level computer, Figure 12 shows a possible mapping of a myofibril model of 32 sarcomeres onto one rack of a Blue Gene*/L (BG/L) supercomputer. (The Blue Gene/L supercomputer was developed through a partnership between IBM and Lawrence Livermore National Laboratory.) The simulation of two thick and eight thin filaments is executed on a dual-core processor. The ratio of thick and thin filaments models the ratio in a half sarcomere. Each thick filament is surrounded by six thin filaments, while each thin filament is surrounded by three thick filaments. This gives a filament ratio of two to one, which must be increased to a ratio of four thin filaments to one thick filament in order to account for the left and right sides of a sarcomere. Note that the thick filaments are double-ended and require roughly twice the computation of the half-thick filament in the preliminary model. We anticipate a 64-thick and 256-thin filament to represent a full sarcomere at the level of a computer node card. A myofibril model can then comprise 32 full sarcomeres modeled at the level of a full rack of Blue Gene/L. The mapping is approximate because the final implementation may require some redistribution of the model among computation units at the level of processors or node cards.

Figure 12

This fuller simulation has an additional computational load due to the communication load that arises from the interaction among sarcomeres. Such models are characterized by a large number of local communications and scale well on toroidal interconnects. (In computing systems, a toroidal network is one in which nodes are connected circularly in more than one dimension. The resulting network topology is a torus, and the network is called toroidal.) On the Blue Gene/L supercomputer, a three-dimensional interconnect torus connects each node to its six nearest neighbors with 1) a link bandwidth of 175 MB/s (bidirectional, 2 bits per cycle); 2) a physically separate global combining/broadcast tree with a bandwidth of 350 MB/s (4 bits per cycle) and a 1.5-ms one-way latency on a 64K-node partition; and 3) a global barrier/interrupt network to allow hardware-based synchronization of large numbers of parallel processors. Such a communication backbone reduces the idle time of the compute nodes due to communication latencies which would be characteristic of such models, and it should drastically reduce the throughput time of the simulations.

The model described here represents an attempt to formulate a mathematical representation of a pair of thick and thin filaments in the sarcomere. The model combines a number of previous modeling efforts to seek a more complete representation than any of the previous models. While the preliminary results compare reasonably well with experimental characterizations of real muscle, much work remains to be done in refining parameter values and verifying the current formulation. Moreover, the current results are generally noisy, with substantial stochastic variation, at least for one pair of filaments. Future extensions of the model to multiple sets of filaments should alleviate much of the stochastic variability because numerous fibers will be working in parallel. While extension of the model to multiple sets of filaments in one sarcomere, and to a multiple-sarcomere structure, will be very computationally demanding, such large-scale calculations have been made possible on recently developed massively parallel platforms such as the Blue Gene supercomputer.

The work of P. P. de Tombe is supported by NIH grants HL62426 and HL75494. The following artwork is used with permission. Lei Jin contributed artwork for the sarcomere proteins in Figures 1 and 2. Dr. Roger C. Wagner contributed the sarcomere artwork in Figure 12. The plots in Figures 9(b) and 10(b) have been reproduced from The Journal of General Physiology, 1995, vol. 105, pp. 1–19 by copyright permission of The Rockefeller University Press.

Table 1 presents equations for free-energy functions and transition rates as functions of the distortion (the distance from the equilibrium position of myosin and binding site on actin). Since the P state is assumed to be the reference, G_P() = 0.0. Hence, free energy is independent of distortion as the myosin is detached. The free energy of P after hydrolysis of ATP is computed assuming ΔG_ATP = −23.0RT. Thus, in Figure 3(b), G_PΔATP is a constant 23RT lower than state P.

Table 1 Equations for free-energy functions and transition rates.

Gp () = 0.0 ΔGATP = −23.0RT μPreF = −8.0RT μF = −16.0RT ρmin = 0.00005 ρmax = 10,000 βPreF = 0.54RT/nm2 βF = 0.54RT/nm2 θP→PreF = 10.0 θPreF→F = 10.0 GPreF() = μPreF + βPreF( − 7.0 nm)2 GF() = μF +βF2 RP→PreF()/RPreF→P() = exp [GP() − GPreF()] RPreF→F()/RF→PreF() = exp [GPreF() − GF()] RF→P()/RP→F() = exp [GF()−ΔGATP] RP→PreF() = max {ρmin, θP→PreF · [RP→PreF()/RPreF→P()]½} s−1 RPreF→P() = min {ρmax, θP→PreF · [RP→PreF()/RPreF→P()]½} s−1 RPreF→F() = min (ρmax, max {ρmin, θPreF→F · [RPreF→F()/RF→PreF()]½}) s−1 RF→PreF() = min (ρmax, max {ρmin, θPreF→F · [RPreF→F()/RF→PreF()]−1/2}) s−1 RF→P() = min (ρmax, max {ρmin, [RF→P()/RP→F()] · RP→F()}) s−1 RP→F() = 0.000051 s−1

The states PreF and F have parabolic profiles because attached crossbridges act like linear springs. The minima of the free-energy profiles of these states are given by μ_PreF and μ_F, respectively. The spring constants for the strongly bound states are given by β_PreF and β_F, respectively. Thus, for this model, the attached crossbridge stiffness is independent of whether it is in state PreF or F. However, there is a 7-nm offset in distortion for the PreF state that results from extension of the myosin head in the pre-rotation state.

The energy profiles defined above are used to set the ratio of transition rates between adjacent states but not the absolute values. The terms θ_P→PreF and θ_PreF→F are scale factors that set the absolute transition rates for the crossbridge cycle. The ρ_min and ρ_max terms are constraints on minimum and maximum values of the transition rates. The constraints prevent extremely large variations in rates that can lead to inaccuracies or poor performance in the numerical methods in the model execution. The arrows in the formulas indicate transitions between states.

*Trademark, service mark, or registered trademark of International Business Machines Corporation.

¹The term PreF refers to a strongly bound but pre-rotation state (essentially the pre-force generating state), while F indicates a strongly bound post-rotation state that induces a strain in the extensible neck region (essentially the force generating state).

Received November 18, 2005; accepted for publication January 4, 2006; Published online September 26, 2006.