Model simulation of the SPOC wave in a bundle of striated myofibrils.

SPOC (spontaneous oscillatory contraction) is a phenomenon observed in striated muscle under intermediate activation conditions. Recently, we constructed a theoretical model of SPOC for a sarcomere, a unit sarcomere model, which explains the behavior of SPOC at each sarcomere level. We also constructed a single myofibril model, which visco-elastically connects the unit model in series, and explains the behaviors of SPOC at the myofibril level. In the present study, to understand the SPOC properties in a bundle of myofibrils, we extended the single myofibril model to a two-dimensional (2D) model and a three-dimensional (3D) model, in which myofibrils were elastically connected side-by-side through cross-linkers between the Z-lines and M-lines. These 2D and 3D myofibril models could reproduce various patterns of SPOC waves experimentally observed in a 2D sheet and a 3D bundle of myofibrils only by choosing different values of elastic constants of the cross-linkers and the external spring. The results of these 2D and 3D myofibril models provide insight into the SPOC properties of the higher-ordered assembly of myofibrils.

Muscle contraction is explained based on a sliding mechanism of thick (myosin) and thin (actin) filaments [1-3]. The sliding of the two filaments occurs through the coordinated interaction between molecular motors (myosin II) and actin filaments. Contraction (On-state) and relaxation (Off-state) are regulated by the concentration of free Ca2+ in vivo [4]. On the other hand, it is known that in the contractile system of a muscle (a skinned muscle model prepared by the removal of membrane systems), spontaneous oscillatory contraction (termed SPOC) occurs under the intermediate activation conditions between the On- and Off-states, which can be controlled by ionic conditions in vitro [5,6]. Thus, the SPOC phenomenon is observed in both skeletal and cardiac myofibrils and in muscle fibers at intermediate Ca2+ concentrations for activation, of which the ionic condition is physiological [7] and also in the presence of high concentrations of MgADP and inorganic phosphate (Pi) with MgATP in Ca2+ free conditions, which are not evident physiologically [8]. We termed the former, Ca-SPOC and the latter, ADP-SPOC [5,6]. The 3D state diagram constructed against the concentrations of Ca2+, Pi and MgADP, in the presence of a fixed concentration of MgATP (covering not only physiological but also non-physiological conditions), showed that SPOC is considered as the third state of the contractile system of muscle; the SPOC region in the state diagram is sandwiched in between contraction and relaxation regions. The Ca-SPOC and ADP-SPOC regions are located at opposite ends of the SPOC region, implying that the Ca- and ADP-SPOCs occur at extreme conditions within the intermediate activation conditions. In the SPOC state, each sarcomere length (SL) changes with a saw tooth waveform, composed of a slow shortening phase and a rapid lengthening phase. The lengthening phase of sarcomeres sometimes travels to adjacent sarcomeres one by one along a myofibril, forming a traveling wave. There is a report showing that the auto-oscillation occurs even in an in vitro pure actomyosin motility system [9].

We recently constructed a unit sarcomere model [10] and a single myofibril model [11] that can capture the main characteristics of SPOC, including the phase diagram consisting of relaxation, contraction, and SPOC states, which shows various patterns of oscillation, i.e., in-phase oscillation, traveling wave, disrupted wave and out-of-phase oscillation, and contraction without oscillation. SPOC studies have two main purposes: one is to clarify how SPOC is self-controlled, i.e., the autonomous control mechanism of SPOC (note that there are several theoretical works in addition to our work on the mechanism of SPOC [12-15]), and the other is to determine whether SPOC has physiological significance, especially, in heartbeat. In the present study, we intended to construct a SPOC model to understand the 2D and 3D patterns of SPOC observed in a bundle of striated myofibrils, i.e., the contractile system of a muscle with the membrane system removed. It is to be noted here that myofibrils are three-dimensionally bundled in muscle fibers; in some cases, however, such as in cultured embryonic muscle cells, myofibrils are arranged relatively two-dimensionally, like a sheet of myofibrils. The present model is an extended version of a single myofibril model previously proposed [11]. That is, the single myofibril model was connected through simple elastic cross-linkers between the adjacent Z-lines and M-lines in the myofibrils, which are aligned side-by-side. To introduce a freedom of motion to the bundle of myofibrils, the end compliance was assumed to exist at one end of each myofibril. Here, the elastic constants for the elastic cross-linkers and the end compliance were represented by kz and ka, respectively. Based on the results of the simulation, we discuss the roles in SPOC of cross-linkers that are introduced in the present 2D and 3D myofibril models.

Results and Discussion

Theoretical framework of the 2D and 3D myofibril models

The present mathematical model is the extension of a single myofibril model previously proposed [11]. The extended parts are twofold: (i) Two to three myofibrils are bundled at the sides of myofibrils to arrange the myofibrils in parallel. As illustrated in Figure 1A, one end of the bundle of myofibrils is fixed to a stiff needle and another end is attached to a flexible needle, X, of which the elastic constant is K. Here, the end compliance of myofibrils is introduced through a simple spring having an elastic constant ka. In practice, the end compliance should exist at both ends of myofibrils, but it was introduced only at one end for simplicity in our model, because the results of simulation were essentially the same. (ii) The Z-lines (and the M-lines) of adjacent myofibrils are connected through an elastic spring having an elastic constant k, in which the elastic force is assumed to work only in a longitudinal direction. The natural length of both springs with the elastic constants, ka and k, is assumed to be 0. The present model is represented by the following set of ordinary differential equations:

Figure 1

Schematic representation of a model on a bundle of myofibrils. (A) Myofibrils are connected side-by-side through an elastic spring having an elastic constant k between adjacent sarcomeres and to a flexible needle X, through an elastic spring ka. Half-sarcomeres are connected alternatively in the series across the Z- and M-lines. (B) Cross-sectional view of a bundle of myofibrils when they are three-dimensionally connected (3D connection). Three myofibrils are equivalent to each other. (C) Cross-sectional view of a bundle of myofibrils when they are two-dimensionally connected (2D connection). The three myofibrils are not equivalent to each other. (D) Schematic illustration showing a unit model composed of a pair of thick and thin filaments, which is taken from Figure 1(b) in [11]. (E) Classification of patterns of SPOC wave obtained in the present model simulation, which is shown from the point of view of phase difference in sarcomere oscillation with laterally and longitudinally adjacent sarcomeres. For more details, see the text.

The suffix, i, means the number of sarcomeres in each myofibril, and j means the number of myofibrils. The total number of sarcomeres in each myofibril is N, and the total number of myofibrils, M. Thus this model is expressed by 3N·M+1 differential equations. The meaning of Eqs. (1)–(4) is explained below in order.

Equation (1) is the time evolution equation for the fraction of attached myosin heads (cross-bridges) to the total number of myosin heads in the i th sarcomere in the j th myofibril, P. Here, we have assumed that myosin heads take only two states, attached and detached states, i.e., the so-called two state model (see Fig. 1D). The transition rates between the two states are given by α (attachment rate) and β (detachment rate). One of the characteristics of our model is that α and β depend on the sarcomere lattice spacing d, which differs for each sarcomere, i.e., d=d [10,11,16], namely, these are expressed by functions of d as α(d) and β(d). Their explicit forms are expressed as the following Eqs. (5) and (6):

where α1, d0 and β0 are some positive constants. Θ(x) is the Heaviside step function defined by Θ(x)=1 for x>0 and Θ(x)=0 otherwise. α(d) represents the attachment rate constant of myosin heads. It is assumed that α linearly decreases depending on its lattice spacing d and reaches 0 at the limit distance above which cross-bridge formation is prohibited because d exceeds the maximum length of myosin heads. We consider that β(d) is also a function of d, but it is assumed as a constant for simplicity. The value of α1 represents the activation level of myofibrils. In muscle, the values of α1 and β0 should depend on the environmental conditions such as the composition of solutions and temperature, and so on. β0 may also depend on how much external force is applied to the cross-bridges. Two extremes are a contraction condition realized at larger α1 and smaller β0, and a relaxation condition realized at smaller α1 and larger β0.

Equation (2) shows the force balance along the long axis. Fex j represents an external force imposed on the j th myofibril, which takes a positive value when it works in the direction for extension. Fex j is expressed as Eq. (7).

where x is the position of X and SL is the length of the i th sarcomere in the j th myofibril. As the value of ka becomes larger, the length of every myofibril approaches x, implying that the length of every myofibril becomes equal. The second and the third terms of Eq. (2) are the active and frictional forces generated by cross-bridges, respectively, where a is the strength of mean active force generated by a single cross-bridge and s0 is a spatial interval between adjacent myosin heads along the thick filament. ξ is the length of overlap between the thick and thin filaments, which is related with the sarcomere length SL as SL=SL(0)–ξ, where SL(0) is the sarcomere length at no overlap. Note in Eq. (2) that the quantity ξP/s0 represents the total number of cross-bridges in the i th sarcomere in the j th myofibril. Although the molecular friction constant, η, is reported to be asymmetrical against lengthening and shortening of sarcomeres, at least under a rigor condition [17], we assumed η as a constant. Finally, the last term of Eq. (2) represents the interaction between the laterally adjacent sarcomeres through the position of the Z-line (or the M-line). Here Z is the position of the Z-line (or the M-line) of i th sarcomere in the j th myofibril, which is expressed as . k in Eq. (2) is the spring constant that laterally connects sarcomeres in adjacent myofibrils. Note that the force derived from an elastic protein such as connectin/titin is not considered in the present model. In skeletal muscle, in practice, it has been shown that connectin/titin is not essential for SPOC to occur [18-20].

Equation (3) shows the force balance along the short axis. Here four types of force are considered. The first term shows the elastic force generated by cross-bridges, where k and l are, respectively, the elastic constant of cross-bridges and the natural length of myosin heads, both of which are some positive constants. The second term is attributable to various interaction effects between the thick and thin filaments, and to the elastic force working to maintain the lattice structures of the Z- and M-lines [16], and osmotic pressure [21], and so on, of which the elastic constant is expressed as k. lr is the natural length of sarcomere lattice spacing in the relaxing condition. Based on the experimental observations showing that the lattice spacing, lr, decreases as the SL increases [16], lr is considered as an increasing function of overlap length ξ between the two myofilaments. The first part in the second term can be expressed as a linear function,

where l0 represents the lattice spacing where there is no overlap between the two myofilaments in the relaxed state, and l1 means the slope of the lattice spacing vs. overlap relationship. The third term in Eq. (3) is the friction force that the filament lattice bears along the short axis, where η is assumed to be a small positive constant. The last term in Eq. (3) represents the interaction between longitudinally adjacent sarcomeres through the lattice spacing. The value of kMZ is the activity coefficient that couples between adjacent sarcomeres along the long axis of myofibrils.

Equation (4) shows the equation of motion of the flexible needle X, of which the position is denoted by x. The first term on the right side in Eq. (4) represents the elastic force which is the result of the elasticity of the flexible needle, and the second term represents the total force coming from springs with ka. Y0 is the specific value of x when the flexible needle is not bent, that is, when no force is imposed on X. By changing Y0, we can control the external force acting on the bundle of myofibrils. For a given Y0, there is a stationary state where all sarcomeres have the same length and the same active force , where ξ̄ and P̄ are stationary solutions of ξ and P obtained by solving the equations (1)–(4) for the given Y0. These quantities Y0, ξ̄ and F0 satisfy the following relation:

where the first term of the right side is the total length of myofibril, and the second and third terms are extensions of the springs ka and K, respectively. M is the number of myofibrils. Thus, the three variables, Y0, ξ̄ and F0 in Eq. (9) are related to each other, and therefore we can specify the values of other variables if one of three values can be determined. Here, the value of F0 will be used for the specification of Y0. The term on the left side in Eq. (4) is the friction force the flexible needle bears, where η is estimated to be a small positive constant.

The meaning and the values of the parameters contained in this model are summarized in Table 1. The initial conditions used here are given by adding a random number δ ∈ [−0.0001, 0.0001] to a stationary solution (P̄, ξ̄, d̄ ) for a given F0. The numerical simulations of this model were performed with sufficiently high precision.

Table 1

Summary of parameters used in the present model simulation

Parameter [dimensiona]	Meaning	Value
a [Fx]	Average active force per each cross-bridge	1 [pN]
s0 [Lx]	Interval between adjacent cross-bridges along the thick filament	0.01 [μm]
η [Fx·T/Lx]	Molecular frictional constant for cross-bridges	1.5 [pN·sec/μm]
km [Fy/Ly]	Lateral elastic constant of cross-bridges	1 [(a.u.)b/nm]
lm [Ly]	Average natural length of cross-bridges	23 [nm]
ηd [Fy·T/Ly]	Frictional constant of the sarcomeric structure	0.0005 [(a.u.)b·sec/nm]
kr [Fy/Ly]	Lateral elastic constant of the sarcomeric structure	60 [(a.u.)b/nm]
α1 [1/T·Ly]	Slope of the attachment rate vs. lattice spacing relationship (the indicator of the activation level in model simulations)	19 [1/sec·nm]
β0 [1/T]	Detachment rate of cross-bridges	20 [1/sec]
d0 [Ly]	Lattice spacing beyond which the activation does not occur	25.3 [nm]
SL(0) [Lx]	SL at no overlap between the thick and thin filaments	3.6 [μm]
lr0 [Ly]	Lattice spacing at no overlap in the relaxed state	23 [nm]
lr1 [Ly/Lx]	Slope of the lattice spacing vs. SL relationship	2.53 [nm/μm]
F0	Constant part of the auxotonic force	10 [pN]
K [Fx/Lx]	Spring constant of the flexible needle X	0.01–10.0d [pN/μm]
kMZ [Fy/Ly]	Elastic modulus of the Z- and M-lines	0.5 [(a.u.)b/nm]
N	Number of sarcomeres for each myofibril	24
M	Number of myofibrils	1–3d
kz [Fx/Lx]	Elastic constant of the cross-linkers connecting the Z- and M-lines	0.0–10.0d [pN/μm]
ka [Fx/Lx]	Elastic constant of the spring that connects the end of each myofibril to a flexible needle X	20–100d [pN/μm]
ηn [Fx·T/Lx]	Frictional constant of the flexible needle X	0.0005 [pN·sec/μm]

Dimensions are represented by [Lx] (for the length along the long axis; μm), [Ly] (for the length along the short axis; nm), [T] (for time; sec), [Fx] (for force along the long axis; pN) and [Fy] (for force along the short axis; arbitrary units).

The unit of force along the short axis is arbitrary (a.u.), because precise measurements of the force have not been performed yet.

In our model, only one pair of the thick and thin filaments is considered. Hence this value corresponds to the contractile force generated by one pair of the filaments at a stationary state. To apply this model to an actual myofibril, this value needs to be multiplied by several thousands.

Variable in this range.

With this model, we considered two types of boundary conditions, a circular connection (periodic boundary condition) (Fig. 1B) and a sheet-type connection (Neumann boundary condition) (Fig. 1C). In the present study, we investigated the characteristics derived from various values of parameters, e.g., a cross-linker stiffness k, a needle’s stiffness K, and boundary conditions. As a result, several SPOC patterns experimentally observed have been revealed in the present model simulation as summarized in Figure 1E.

Model simulation of SPOC patterns in a single myofibril

First, to confirm that the present model contains the single myofibril model previously presented [11], we performed a model simulation for the case of a single myofibril. The difference between the two models is that the present model (Fig. 2A) contains end compliance, whereas the previous one does not (Fig. 2B).

Figure 2

Examples of model simulation showing typical patterns of sarcomere length (SL) oscillations for a single myofibril in which 24 sarcomeres are connected in series. (A) Schematic representation of a myofibril in the present model. (B) Schematic representation of a myofibril of a previous model [11]. (C) Phase diagram of the SL oscillation patterns of a myofibril obtained from the present model simulation. Ordinate and abscissa, the spring constant of the external spring K and the external force F0, respectively. (D) A correspondence table of colors in (C). (E)–(H) Typical examples showing a phase relationship for various SPOC patterns obtained in the present model simulation: an in-phase (E: F0=20, K=0.01), a traveling wave (F: F0=20, K=0.2), a disrupted traveling wave (G: F0=20, K=1.0), an out-of-phase (H: F0=20, K=3.8). The abscissa is the number of sarcomeres, and the ordinate is an oscillation phase in each sarcomere.

The previous model should be equivalent to the case when the value of ka in the present model is extremely large. Then, the comparison was made with large ka, and F0=10 and 20 (previously used). As a result, we could confirm that the present model is consistent with the previous one. As the value of K increases, all of the SPOC patterns presented in the phase diagram (Fig. 3 of [11]), i.e., in-phase synchronization, traveling wave, disrupted traveling wave and out-of-phase synchronization, appeared in this order. That is, the result shown in Figure 2C corresponds to that at N=24 in Figure 3 of [11], and the results shown in Figure 2E, F, G and H, respectively, correspond to those in Figure 2(a), (b), (c) and (d) of [11].

Figure 3

Phase diagram and typical examples of the SL oscillation patterns obtained from the present model simulation in the case of two myofibrils connecting in parallel. (A) Schematic representation of the model. Two myofibrils are equivalent to each other. (B)–(E) Phase diagram of the SL oscillation patterns from a point of view of elastic constants k and ka when K=3.0 (B), K=1.0 (C), K=0.2 (D), K=0.01 (E). (F) A correspondence table of colors in (B)–(E). (G-1)–(K-3) Typical examples of in-phase [(G), K=0.01, ka=100, k=5.0], Type-1 traveling [(H), K=0.5, ka=100, k=0.5], disrupted [(I), K=3.0, ka=20, k=5.0], in-phase traveling [(J), K=1.0, ka=100, k=10.0], and in-phase disrupted [(K), K=3.0, ka=20, k=10.0] patterns. (G-1)–(K-1) Typical examples showing a phase relationship of SL oscillation for various SPOC patterns obtained in the present model simulation, where the abscissa is the number of sarcomeres, and the ordinate is an oscillation phase of each sarcomere. These phase relationships were obtained from the snapshot at t=396 as shown by red arrows in (G-3)–(K-3). (G-2)–(K-2) Typical examples showing time course of the changes in x(t). (G-3)–(K-3) show enlarged views of (G-2)–(K-2).

Model simulation of SPOC patterns in a bundle composed of two myofibrils

Next, we examined how the SPOC patterns appear when two myofibrils are aligned through the elastic spring having k, connecting the adjacent Z- and M-lines as illustrated in Figure 3A. In the case of a bundle composed of two myofibrils, 3D and 2D connections are equivalent to each other, except that the number of springs in the 3D connection is twofold that in the 2D connection. Figure 3B–E are a type of two dimensional phase diagram showing how the SPOC patterns appear against different values of ka and k at fixed values of K (respectively, 3.0, 1.0, 0.2, 0.01). Different from the mode for a single myofibril (Fig. 2), the relationship of the oscillation phase between laterally adjacent sarcomeres, consisting of different myofibrils, can be obtained. Figure 3G–K show typical examples of (1) the relationship in the oscillation phase of sarcomeres in two myofibrils, (2) the time course of the change in amplitude of length oscillation in the bundle of myofibrils, and (3) enlarged views of (2).

The SPOC patterns observed in this case are in-phase SPOC (Fig. 3G), Type-1 traveling wave (Fig. 3H), disrupted wave (Fig. 3I), in-phase traveling wave (Fig. 3J), and inphase disrupted wave (Fig. 3K). The noticeable characteristics of these phase diagrams were: 1) the in-phase SPOC pattern tends to appear when both k and K are small, and 2) the disrupted pattern tends to appear as K becomes larger (and ka is larger), which is because the total length of myofibrils tends to be fixed, so that the change in the SL is restricted by this boundary condition. It was also noticeable that the traveling wave region is smaller when compared with the case of a single myofibril (Fig. 2C) at k=0. In the case of a disrupted wave observed at k=0 (Fig. 3B, C), a wave locally traveling in an opposite direction is sometimes observed, but, as k is introduced, the traveling wave consistently travels in one direction, suggesting that k has a function to match the direction of traveling waves in the two myofibrils.

As the value of k is further increased, the in-phase traveling wave becomes dominant (Fig. 3B–D). At larger k, the phase difference of oscillation between the adjacent sarcomeres tends to become smaller, so that the oscillation patterns tend to be in order. It is noted that the disrupted wave region is mixed here and there with the in-phase traveling wave region (Fig. 3C). This result suggests that some uncertainty exists in the present model simulation; one reason is that we have introduced a randomness, δ ∈ [−0.0001, 0.0001], for the initial values of parameters because there must be some fluctuation in the parameters, e.g., SL dispersion (in practice, the SL in heart muscle has a distribution of about 10% even in the relaxing conditions [22]). Therefore, the initial conditions of this model simulation are always slightly different from each other, which may yield different SPOC patterns. Another possibility to explain this uncertainty is that the present model simulation may essentially contain this uncertainty due to its non-linearity. This problem should be clarified in the future.

Regarding the effects of ka, it is noticeable that, on increasing the ka value, the SPOC patterns that appear at a region of smaller k values shift toward a larger k region, whereas the region of the in-phase traveling wave shifts toward a smaller k region. This means that the region of Type-1 traveling wave becomes narrow, and the sarcomere oscillation tends to be in phase between the adjacent myofibrils as ka becomes larger.

In the case of in-phase SPOC, the sarcomeres in each myofibril were consistently synchronized; it is however conceivable that the oscillation phase between the adjacent two myofibrils is unsynchronized. But an unsynchronized in-phase SPOC was not found in the present model simulation.

Another noticeable point is that the in-phase SPOC basically occurs when the values of both K and k are small, but when the value of K is small and the value of k is large, we found that the in-phase SPOC region again appears (Fig. 3E). When k is fixed at a large value, e.g., 10.0 in Fig. 3B–E, and as the value of K is increased, the SPOC patterns appear in the order of in-phase oscillation, in-phase traveling wave, and in-phase disrupted wave, which is similar to the case of single myofibrils (Fig. 2). Thus, this result is considered to be attributable to the fact that the two myofibrils behave just like a single myofibril because of the large mechanical coupling between the laterally adjacent sarcomeres.

Finally, it is to be noted that the oscillation amplitude of the total length of myofibrils largely depends on the SPOC patterns: in the case of in-phase SPOC (Fig. 3G-2) the amplitude of the total length of myofibrils becomes the sum of the oscillation amplitude of each sarcomere, such that the oscillation amplitude of a bundle of myofibrils becomes maximal (in the present simulation, the peak-to-peak amplitude reached about 15% of the total length, although this depends upon the activation conditions). In the other cases, the oscillation amplitude of the total length of a bundle quickly decays starting from the initial in-phase oscillation with a large amplitude. Only in the case of the in-phase traveling wave, the oscillation amplitude remains at about 3% of the total length of myofibrils at steady state. If the SPOC dynamics are involved in the heartbeat, the total length of myofibrils should change to some extent. In this sense, our interests are particularly focused on the conditions where the in-phase SPOC and the in-phase traveling wave of SPOC occur. Also, it is interesting that the total length of a bundle largely fluctuates in the case of disrupted SPOC, which may correspond to a pathological state.

Model simulation of SPOC patterns in a bundle composed of three myofibrils: 1) the case for a 3D (triangular) connection

Next, we examined how the SPOC patterns appear in the case of a bundle composed of three myofibrils for different combinations of the values of K, ka and k. First, the case of three myofibrils connected three dimensionally was examined. As a result, as shown in Figure 4, we found a new SPOC state in addition to those observed in the bundle of two myofibrils. i.e., the traveling velocity of the lengthening phase of SL in one myofibril is different from that in the other two myofibrils, so that the oscillation pattern does not maintain a steady state, but rather it changes with time. We termed this traveling wave as “Type-2”, and the traveling wave that keeps the mutual phase difference nearly constant as “Type-1”. The Type-2 traveling wave tends to appear at larger values of k and ka (Fig. 4B–E).

Figure 4

Phase diagram of the SL oscillation patterns obtained from the present model simulation in the case of sterically connecting three myofibrils and typical examples of each oscillation patterns. (A) A schematic representation of the present model, where three myofibrils are equivalent to each other. (B)–(E) Phase diagram of the SL oscillation patterns from the point of view of elastic constants k and ka when K=3.0 (B), K=1.0 (C), K=0.2 (D), K=0.01 (E). (F) A correspondence table of colors in (B)–(E). (G-1)–(L-3) Typical examples of in-phase [(G), K=0.01, ka=100, k=10.0; see Supplementary movie 1 “In-phase.avi”], Type-1 traveling [(H), K=0.2, ka=20, k=1.0; see Supplementary movie 2 “Type-1 traveling. avi”], Type-2 traveling [(I), K=1.0, ka=100, k=1.0; see Supplementary movie 3 “Type-2 traveling.avi”], disrupted [(J), K=3.0, ka=100, k=0; see Supplementary movie 4 “disrupted.avi”], In-phase traveling [(K), K=0.2, ka=20, k=10.0; see Supplementary movie 5 “In-phase-traveling.avi”], and In-phase disrupted [(L), K=3.0, ka=20, k=10.0; see Supplementary movie 6 “in-phase-disrupted.avi”] patterns. (G-1)–(L-1) Typical examples showing a phase relationship of SL oscillation for various SPOC patterns obtained in the present model simulation, where abscissa is the number of sarcomeres, and ordinate is an oscillation phase of each sarcomere. These phase relationships were obtained from the snapshot at t=396 as shown by red arrows in (G-3)–(L-3). (G-2)–(L-2) show typical examples of time course of the changes in x(t). (G-3)–(L-3) show enlarged views of (G-2)–(L-2).

The phase diagram obtained in the present case is basically similar to the case of two myofibrils: that is, at k=0, the SPOC pattern changes from the in-phase SPOC to the disrupted SPOC as the value of K is increased. When the k is introduced and increased, Types 1 and 2 traveling wave SPOCs and then the in-phase traveling wave appear in this order independent of the value of ka. And, similarly to the case of two myofibrils, the in-phase SPOC appears at large k and small K. This is likely because a bundle of three myofibrils behaves just as a single myofibril. Again, when the value of ka is small, the states of SPOC are shifted toward the side of larger k, which is a similar occurrence to that in the case of two myofibrils. One of the characteristics of the Type-2 traveling wave is that the total length of the bundle oscillates with some period (about 50 in the case of Fig. 4I-2) and a small amplitude. Such a traveling pattern in which the pattern dynamically changes with time is often observed experimentally.

Model simulation of SPOC patterns in a bundle composed of three myofibrils: 2) the case for a 2D (planar) connection

Finally, we examined how the SPOC patterns appear in a bundle composed of three myofibrils for a planar (2D) connection between sarcomeres in adjacent myofibrils (Fig. 5A), again under various combinations of the values of K, ka and k (Fig. 5B–E). In muscle tissue which is composed of a bundle of muscle cells, the myofibrils, the major components of striated muscle, are three-dimensionally arranged through cross-linkers connecting the Z- and M-lines. But, in some cases such as cultured embryonic muscle cells, myofibrils are rather arranged two-dimensionally as a sheet of myofibrils. So, it is worth examining the case of 2D connection [23,24].

Figure 5

Phase diagram of the SL oscillation patterns obtained from the present model simulation in the case of planar connecting three myofibrils. (A) A schematic representation of the present model, where three myofibrils are not equivalent to each other. (B)–(E) Phase diagram of the SL oscillation patterns from the point of view of elastic constants k and ka when K=3.0 (B), K=1.0 (C), K=0.2 (D), K=0.01 (E). (F) A correspondence table of colors in (B)–(E). (G-1)–(G-3) Typical example of Type-2 traveling wave (K=3.0, ka=20, k=3.0). (G-1) shows a phase relationship of SL oscillation obtained at t=376 in (G-3) as shown by a red arrow, where the abscissa is the number of sarcomeres, and the ordinate is an oscillation phase of each sarcomere. (G-2) shows a typical example of time course of the changes in x(t). (G-3) shows an enlarged view of (G-2).

As summarized in Figure 5, the results are essentially similar to those for the 3D connection (Fig. 4). One major difference is that in the phase diagram the region of Type-2 traveling wave is significantly larger than that in the 3D connection. This is understandable because the three myofibrils are not equivalent to each other, that is, the central one is different from the two peripheral ones, so that the in-phase traveling wave and the Type-1 traveling wave hardly occur. In the case of the Type-2 traveling wave, the central myofibril tends to show a higher traveling velocity, resulting in a more distinct periodicity of oscillation of the total length of a bundle as shown in Figure 5G. In the case of the 3D connection (Fig. 4), the in-phase traveling wave appears as the value of k is increased, but in the 2D connection case, only the disrupted pattern appears. And, the Type-1 traveling wave appears only in the central small region of the phase diagram (Fig. 5B–E). Further, the Type-1 traveling wave is divided into two types, i.e., the oscillation phase of the central myofibril which goes either ahead (a convex-type) or behind (a concave-type) the other two peripheral myofibrils.

An example showing the SPOC pattern observed in a myofibril bundle

An example of ADP-SPOC pattern, experimentally observed, is shown in Figure 6. A phase-contrast microscopy shows that the Type-1 traveling wave occurs in a bundle of glycerinated skeletal myofibrils prepared from rabbit psoas muscle (see Supplementary movie 7 “movie_1.avi”); this is a Type-1 traveling wave keeping the phase difference of SL oscillation nearly constant among adjacent myofibrils. The condition for this observation was as follows: KCl 120 mM, MgCl2 4 mM, MOPS 20 mM, EGTA 4 mM, ATP 0.1 mM, ADP 4 mM, Pi 4 mM, pH 7.0 at room temperature [5,8].

Figure 6

An example showing a Type-1 traveling wave in a bundle of skinned myofibrils. (A) Phase contrast image of rabbit psoas glycerinated myofibrils under an ADP-SPOC condition [8]. Numbers shown by arrows indicate the sarcomeres used for image analyses shown in (B) and (C). The yellow outlined rectangle on the upper side defines “Top”, and on the lower side, “Bottom”. (B) Time course of SL changes. The numbers of sarcomeres correspond to those indicated by the arrows in (A). The “Top” and “Bottom” also correspond to the rectangle regions in (A). (C) Time course of the SL changes in single sarcomeres in the “Top” and “Bottom” rectangle regions. The red (Top) and blue (Bottom) lines, respectively, correspond to the red and blue arrows in (A), and the red and blue lines in (B).

Characteristics of SPOC patterns of a bundle composed of three myofibrils obtained in the present model simulation

We can summarize the results as follows (Fig. 7):

Figure 7

Schematic representation of the summary on the phase diagram obtained by the present 2D and 3D myofibril model simulations. To construct this summary, some minor uncertainties discussed in the text were ignored.

In-phase SPOC occurs at the regions where the value of K is small, and the value of k is either small or large. As the value of ka is smaller, the in-phase SPOC region becomes larger. This characteristic is independent of the 2D and 3D connections (Figs. 4 and 5).

Disrupted SPOC occurs at the region where the value of K is large and the value of k is small. In the 3D connection (Fig. 4B–E), it is to be noted that the disrupted SPOC appears here and there at larger k (as schematically shown by several green circles in Fig. 7), although the appearance probability is low. This may be attributable to the uncertainty of this non-linear model, which includes some randomness on the initial values of parameters as ±δ. In the 2D connection, on the other hand, the disrupted SPOC region appears at both smaller and larger values of k, especially at larger K (Fig. 5B–E).

Type-1 (phase-locked) and Type-2 (phase-unlocked) traveling SPOCs occur at intermediate values of k (Fig. 4B–E). In-phase traveling SPOC occurs at larger k. Thus, as the value of k is increased, Type-1, Type-2 and in-phase traveling SPOCs appear in this order.

We have newly introduced the parameters, k and ka, in the present myofibril bundle model to extend a single myofibril model previously presented [11]. The spring, of which an elastic constant is k, connecting sarcomeres belonging to different myofibrils tends to reconcile the oscillation phase of connected sarcomeres in order. The spring having the elastic constant ka, introduced as an end compliance of each myofibril, contributes to the coordination of SPOCs in the bundle of myofibrils. That is, as elastic constant ka is increased, the total length of myofibrils tends to become equal, so that ka contributes to the coordinated oscillation of myofibrils. The parameter K, which is experimentally controllable, is important for global coupling in the mechanism of SPOC as previously described in detail [11].

A possible physiological significance of SPOC patterns obtained in the present model simulation

Among several SPOC patterns obtained in the present simulation, the pattern showing the largest amplitude of oscillation of the total length of myofibrils is the in-phase SPOC, and the second largest one is the in-phase traveling SPOC. If the SPOC properties have some functions in heartbeat, those patterns are candidates having a physiological role. In other words, the SPOC patterns in which the total length of myofibrils is not oscillated or randomly fluctuated because of the out-of-phase or disrupted oscillation may be non-physiological or pathological. Those discussions are still speculative, but it will be worth examining, based on the present model simulation, what is a condition for organized SPOC or disordered SPOC to occur, and how the transition from a normal (healthy) SPOC to an abnormal (pathological) SPOC occurs.

The present model only provides a small seed for such a line of research approach for the future. But, the present model simulation suggests that the disorganized SPOC tends to occur as 1) the value of K is increased, which may correspond to the stiffening of muscle tissue, and/or 2) the value of k comes out of an optimal region, which may be due to stiffening or breakage of the cross-linkers between myofibrils.

In the present model, we introduced mainly three parameters, i.e., K, ka and k. Among these parameters, the value of K is experimentally controllable. It may be possible to control the value of k by selective cleavage with proteolytic enzyme. Thus, there is a possibility that the present model simulation presents a guideline for future experimental studies on muscle physiology.

SPOC (SPontaneous Oscillatory Contraction) is a phenomenon observed in striated muscle under the intermediate activation conditions. To understand the SPOC properties in a bundle of myofibrils, we have extended the single myofibril model to a two-dimensional (2D) and a three-dimensional (3D) model, in which myofibrils are elastically connected side-by-side through cross-linkers between the Z-lines and the M-lines. These models could reproduce various patterns of SPOC waves experimentally observed in a bundle of myofibrils by choosing appropriate values of elastic constants of the cross-linkers and external spring.