Modeling and simulation of excitation- contraction coupling of fast-twitch skeletal muscle fibers

2.1Modeling of the skeletal fast muscle fiber membrane potential

A two-compartment model was used to simulate the action potential of skeletal fast muscle fibers. The total ionic current on the surface of the muscle fibers is calculated by adding 𝑁𝑎+ current (I𝑁𝑎,s), K+ delayed rectifier current (I𝐷𝑅,s), K+ inward rectifier current (I𝐼𝑅,s), 𝐶𝑙- current (I𝐶𝑙,s), and 𝑁𝑎+-K+ current (I𝑁𝑎𝐾,s), as shown by Eq. (1):

Due to the difference between ion channel densities in the T tube and myolemma, the ratio of T tube membrane channel density to myolemma channel density can be represented by η. The calculation of ionic current (I𝑖𝑜𝑛𝑖𝑐,t) on the T tube per unit area is then shown by Eq. (2) as:

where g𝑁𝑎,t is the 𝑁𝑎+ channel conductance, g𝐷𝑅,t is the K+ channel conductance, g𝐼𝑅,t is the K+ inward rectifier channel conductance, and g𝐶𝑙,t is the 𝐶𝑙- channel conductance.

The changes inconcentration of 𝑁𝑎+⁢(𝑁𝑎i) in cells, 𝑁𝑎+⁢(𝑁𝑎t) in T tube, and 𝑁𝑎+⁢(𝑁𝑎e) in intercellular space are shown in Eqs (3)–(5), using time as the basis. Similarly, the changes inconcentration (with respect to time) of K+⁢(Ki) in cells, K+⁢(Kt) in T tube, and K+⁢(Ke) in intercellular space can be obtained.

In Eqs (2)–(5), fT represents the proportion of fiber occupied by the T tube, ξ1 represents the ratio of volume to surface area of the T tube, ξ2 represents the ratio of volume to surface area of cells, ξ3 represents the ratio of volume to surface area of intercellular space, τ𝑁𝑎1 represents diffusion time constant of the T tube, and τ𝑁𝑎2 represents the diffusion time constant of intercellular space.

2.2Modeling of calcium cyclingin skeletal fast muscle fiber

We divided the sarcomere into six geometric regions, which are V1, V2, V3, V4, V5, and V6, as shown in Fig. 1. V1=5.5%⁢V, V2=3.5%⁢V, V3=6%⁢V, V4=85%⁢(1-0.75lx)⁢V, V5=85%⁢(1.75lx-1)⁢V, and V6=85%⁢(1-1lx)⁢V.

Figure 1.

Geometric division diagram of half sarcomere.

A 10-state model is used to describe the release process from the T tube voltage to the Reynolds channel. The ten states consist of five states of four voltage sensor molecules and two states of a Reynolds channel (including five closed states (C0-C4) and five open states). Based on the division of the half-sarcomere regions and the 𝐶𝑎2+ cycle, the change in 𝐶𝑎2+ concentration in V1-V6 regions are analyzed, as shown in Eqs (6)–(11), and the model’s parameters are shown in Table 1:

2.3Modeling of the cross-bridge dynamic in skeletal fast muscle fiber

A 6-state model is used when 𝐶𝑎2+ and troponin combine in the non-coincidence region (V4), and the binding site of free troponin (T0) is shown in Eq. (12):

In the non-coincidence region (V4) the function of combining 𝐶𝑎2+ and troponin, F1⁢(𝐶𝑎4,Tn), is calculated according to Eq. (13):

In the coincidence region (V5), an 8-state model is used. During the process of cross-bridge cycling, the binding site of free troponin (T0) is shown in Eq. (14):

In the coincidence region (V5), during the process of cross-bridge cycling, the function of combining 𝐶𝑎2+ and troponin, F2⁢(𝐶𝑎5,Tn), is calculated according to Eq. (15), and the model’s parameters are shown in Table 2.

2.4The fatigue modeling of metabolismin skeletal fast muscle fiber

During the cross-bridge cycle, when the cross-bridge attains a strong binding state from a weak binding state, the ATP will be hydrolyzedto generate ADP and Pi. When the product of solubilities of P𝑖𝑆𝑅 and 𝐶𝑎1 in sarcoplasmic reticulum exceeds 6 mM2, it is considered that the Pi in sarcoplasmgoes through the passive channel and is transported to the sarcoplasmic reticulum,at a speed ranging from 30 to 170 μm/s, and combines with 𝐶𝑎2+ in the sarcoplasmic reticulum to generate precipitation. A Pi of 20 mM can reduce 29% of 𝐶𝑎2+ released by sarcoplasmic reticulum. When the muscle fiber contraction ends, Pi in sarcoplasm can slowly be removed, and the sarcoplasm would then recover to a resting state of 3 mM.

3.1Analysis of the simulated excitation-contraction model of skeletal fast muscle fiber in conditions of non-coincidence of the thick and thin myofilaments

When skeletal fast muscle fibers are in a resting state, the number of ions in the myofibrilremains constant, and their concentrations are shown in Table 3. Calsequestrin in the sarcoplasmic reticulum combines a large amount of calcium ions to be stored. Magnesium ions in sarcoplasm combine at most of the ATP binding sites. The percentages of calcium and magnesium ions in binding sites, respectively, are shown in Table 4.

Figure 2.

Release rate of calcium ion in terminal cistern.

Figure 3.

The simulation module of the binding of calcium ion and troponin.

Figure 4.

The simulation system of the excitation-contraction model of fast-twitch fiber of skeletal muscle when non-coincidence of the thick and thin myofilaments occurs.

When a single action potential is transmitted to the transverse tubule membranes, the dihydropyridine voltage receptor (DHPR) on the transverse tubule membranes is coupled to the Reynolds channel on the terminal cistern to control the release of 𝐶𝑎2+ on the terminal cistern. According to the mathematical model established in this study, a simulation module is built in Simulink to simulate the release of 𝐶𝑎2+ on the terminal cistern membrane under the control of a single action potential. The simulation results are shown in Fig. 2. The terminal cistern began to release the 𝐶𝑎2+ 1.4 ms after the action potential was generated, the peak value of the release flow was 204 μM/ms, and the time corresponding to the peak value was 2.8 ms. The simulation results were the same asthe mean value detected by Baylor et al., upon calculation of the 𝐶𝑎2+ fluoresce in of 11 muscle fibers. When 𝐶𝑎2+ was released into the sarcoplasm, it bound to the binding sites of the buffer in the sarcoplasm, such as ATP and parvalbumin. In the V4 region, besides ATP and parvalbumin, troponin was also included. Considering the V4 region as an example, the simulation process carried out in Simulink can be described as follows: (1) in the V4 region, the calcium ion and troponin are bound with a 6-state model; (2) a blank model editing window is opened in Simulink; (3) the corresponding modules are created; (4) the parameters of the module are set and the values are shown in Table 4; and (5) the modules are connected, as shown in Fig. 3. The simulation processes for the other regions are similar to those employed for the V4 region. The modules built in Simulink are connected to form the excitation-contraction model of fast-twitch skeletal muscle fibers when non-coincidence of the thick and thin myofilaments occur. The module system is shown in Fig. 4. The system parameters of the model are set, the start time of the simulation is set to 0, the stop time of the simulation is set to 30, and the type is set to Variable-step [15].

Figure 5.

The average concentration of free calcium ions in sarcoplasm.

Figure 6.

Changes in the concentration of the free troponin binding sites and the binding of troponin and calcium ions in V4.

The simulation results showed that the concentrations of calcium ions in V3, V4, and V5 regions increased after 1.4 ms, and the concentrations of free calcium ions in V3 reached a peak value of 24 μM at 3.5 ms. With the diffusion of free calcium ions, the concentrations of calcium ions in V4 and V5 also reached peak values. When the action potential was stopped, the calcium pump began to transport the calcium ions from the sarcoplasm to the sarcoplasmic reticulum, which led to the decrease of the calcium ion concentration in the sarcoplasm, upon which it attained resting state. As shown in Fig. 5, Baylor and Hollingworth measured the change of calcium ion concentration by injecting the calcium ion indicator, furaptra, into the fast-twitch muscle fibers of mice. The peak value of free calcium ion concentration measured was 17 μM, and the time corresponding to the peak value was 4.2 ms. Baylor’s simulation by a multi-compartment model showed that the average free calcium ion concentration in sarcoplasm was 16.3 μM, and the corresponding time was 3.5 ms. The simulation results in this study show that the average free calcium ion concentration in sarcoplasm is 16 μM, and the time corresponding to the peak value is 3.2 ms. The free 𝐶𝑎2+ concentration amplitude in this study are similar to the experimental results and Baylor’s simulation results, and the time taken for free calcium ion concentration to reach the peak value is slightly smaller than that in Baylor’s study. The small size of the half-sarcomere used in this model may have led to the free calcium ion concentration in the sarcoplasm reaching the peak value more quickly. In the region of V4, 𝐶𝑎2+ binds to troponin. Since the thick and thin myofilaments do not overlap, there is no cross-bridge activation, and the 6-state model is used to combine 𝐶𝑎2+ and troponin. The simulation results of the binding of 𝐶𝑎2+ and troponin in region V4 are shown in Fig. 6. With the combination of 𝐶𝑎2+ and troponin, the concentration of tropomyosin-binding site of free troponin decreased, and its lowest value was approximately 170 μM, which was 71% of the concentration of free binding site of troponin in the resting state. The concentration of combining 𝐶𝑎2+ and troponin increased, and its maximum value was found to be 70 μM, which was 29% of the combination of calcium ions and troponin at resting state, and the corresponding time was 7 ms. Zot et al. simulated the combination of 𝐶𝑎2+ and troponin, and the peak value of free troponin binding site was 75% of that in its resting state. The peak value of the combination of calcium ion and troponin was 25% of that in its resting state, and the time corresponding to the peak value was 10 ms. In comparison to the simulation results of Zot et al. [16], the peak value obtained in this study is slightly higher, and the time corresponding to the peak value is a little shorter. However, the time corresponding to the peak value in this study is closer to that of Baylor et al. The simulation results of the relationship between troponin and tropomyosin regulated unit concentration is shown in Fig. 7. Only a small number of tropomyosin structures have conformational changes when a single action potential is stimulated.

Figure 7.

Changes in the concentration of troponin and tropomyosin regulated the units in V4.

Figure 8.

Changes in the concentration of troponin and tropomyosin regulated the units.

3.2Simulation analysis of the excitation-contraction model of skeletal fast muscle fiber in the condition of coincidence of the thick and thin myofilaments

When the half-sarcomere’s length was lx∈[1.1,1.75)μm, the thick and thin myofilaments coincided. With the shortening of the half-sarcomere, the number of the coincident thick and thin myofilaments increased; at that time, the structure of the half-sarcomere was divided into 6 regions. In region V4, a 6-state model was used when 𝐶𝑎2+ and troponin combines. In region V5, due to the coincidence of thick and thin myofilaments, when troponin was combined with two 𝐶𝑎2+ ions, the cross-bridge was activated, so an 8-state model was used in this region.

When lx= 1.1 μm, the thick and thin myofilaments coincided completely. The conformations of troponin and tropomyosin regulated units are changed by the binding of 𝐶𝑎2+ and troponin during a single contraction. The head of myosin (cross bridge) on the thick myofilament binds to the myosin-binding site on the thin myofilament, actin, by hydrolyzing ATP. The simulation results of the combination of 𝐶𝑎2+ and troponin, and cross-bridge activation are shown in Fig. 8. In comparison to Fig. 7, the concentration of troponin and tropomyosin regulated units increased, and the combination of myosin to actin changed the rate of binding and separation of 𝐶𝑎2+ and troponin. This promoted the conformation change of tropomyosin and decreased the rate of separation of 𝐶𝑎2+ and troponin. The result was consistent with that of Zot’s. The simulation results of the degree of cross-bridge activation are shown in Fig. 9. The degree A2 of cross-bridge activation represents the intensity of the excitation-contraction of the skeletal fast muscle fibers. The simulation results of ATP hydrolysis are shown in Fig. 10, with the circulation of the cross-bridge, and the phosphoric acid produced by hydrolysis of ATP accumulating continuously.

Figure 9.

Degree of cross-bridge activation.

Figure 10.

ATP hydrolyzed to produce phosphate.

When the length of the half-sarcomere is lx∈(1.1,1.75), the model can be used to simulate the change of free calcium ion concentration in the sarcoplasm, the degree of cross-bridge activation, and the Pi produced by ATP hydrolysis in the conditions of different degrees of coincidence of thick and thin myofilaments.