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

BACKGROUND: The current excitation-contraction coupling model of fast-twitch skeletal muscle fibers cannot completely simulate the excitation-contraction process.
OBJECTIVE: To solve this problem, this study proposes an excitation-contraction model of fast-twitch skeletal muscle fibers based on the physiological structure and contractile properties of half-sarcomeres.
METHODS: The model includes the action potential model of fast-twitch fiber membranes and transverse tubule membranes, the cycle model of 𝐶𝑎2+ in myofibril, the cross-bridge cycle model, and the fatigue model of metabolism.
RESULTS: Finally, detailed analyses of the results from the simulation are conducted using the Simulink toolbox in MATLAB. Two conditions, non-coincidence and coincidence, are analyzed for both the thick and thin myofilaments.
CONCLUSIONS: The simulation results of two groups of models are the same as the previous research results, which validates the accuracy of models.

Introduction

At present, the action potential of the sarcolemma and transverse tubule membranes are studied using the Hodgkin-Huxley model and the Goldman-Hodgkin-Katz equation. The inwardly rectifying potassium (Kir) channel, together with the chloride channel (ClC-1), plays a vital role in the skeletal muscle physiology. Current research primarily focuses on the fibrous surface and transverse tubular system (TTS) of amphibians [1], and the measurement of Kir of mammalian skeletal muscle fibers [2, 3, 4]. The primary physiological function of mammalian skeletal muscle is to maintain the stability of resting membrane potential and to enable the excitation of electrogenic cells. It is essential to study the properties of Kir in mammalian skeletal muscle fibers to further understand the related muscular diseasesin humans [5, 6, 7].

In vertebrate skeletal muscle fibers, action potential controls contractile activity by inducing rapid changes of the free concentration in the sarcoplasm. , thus, plays a role in triggering and regulating skeletal muscle contraction, and its concentration affects contractile force and speed. In order to regulate in a myofibril accurately, according to the physiological structure properties of the sarcomere, researchers have used half-sarcomeres and divide its space [8] to analyze the change of concentration in the corresponding regions [9, 10, 11]. The indicator, with low affinity, is microinjected into the sarcoplasm to evaluate the amplitude and time historyof regional spatial variation of mean concentration [12]. This model can be used to evaluate the movement of in myofibrils, the release of in sarcoplasmic reticulum, the combination of in sarcoplasm, and primary buffers (such as troponin, ATP, parvalbumin, and sarcoplasmic reticulum pump). The in sarcoplasm is recycled to the sarcoplasmic reticulum by the pump.

Skeletal muscle fibers exhibit contractility through the relative sliding of the sarcomere’s thick and thin myofilaments in the myofibril. The sliding process depends on the hydrolysis coupling of ATP and the circulation of the fine filaments, i.e., actin filaments, along the molecules on the head of myosin, forming cross-bridges. In striated muscle fibers, it is believed that tropomyosin and troponin prevent the binding of myosin and actin, which form the cross-bridge circulation and act as doormen. Tropomyosin is the basis of many of actin’s biological activities, and the movement of troponin on the surface of actin is considered critical to the cooperative allosteric regulation of actin. The different positions of the thick and thin myofilaments will affect the speed at which and troponin combine and separate [13]. Skeletal muscle fibers become fatigued during the excitation-contraction process, and many factors affect fatigue. Metabolic fatigue is a multi-factored regulation process that includes the accumulation of phosphate, cross-bridge circulation, and reduction of calcium in the sarcoplasmic reticulum [14].

The currently existing excitation-contraction model of fast-twitch muscle fibers cannot simulate the excitation-contraction process completely. To solve this problem, this study proposes an excitation-contraction model for fast-twitch skeletal muscle fibers based on the physiological structure and contractile properties of the half-sarcomere. The study is structured as follows: (1) establishment of an action potential model for fast-twitch muscle fiber membranes and transverse tubule membranes; (2) establishment of the cycle model of in myofibril; (3) establishment of the cross-bridge cycle model and the fatigue model of metabolism in half-sarcomeres; and (4) analysis of the results from the excitation-contraction model simulation.

Methods

Modeling 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 (), delayed rectifier current (), inward rectifier current (), current (), and current (), 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 () on the T tube per unit area is then shown by Eq. (2) as:

where is the channel conductance, is the channel conductance, is the inward rectifier channel conductance, and is the channel conductance.

The changes inconcentration of in cells, in T tube, and in intercellular space are shown in Eqs (3)–(5), using time as the basis. Similarly, the changes inconcentration (with respect to time) of in cells, in T tube, and in intercellular space can be obtained.

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

Modeling of calcium cyclingin skeletal fast muscle fiber

We divided the sarcomere into six geometric regions, which are , , , , , and , as shown in Fig. 1. , , , , , and .

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 () and five open states). Based on the division of the half-sarcomere regions and the cycle, the change in concentration in regions are analyzed, as shown in Eqs (6)–(11), and the model’s parameters are shown in Table 1:

Table 1

Model parameters

Parameter	Unit	Value	Parameter	Unit	Value
𝑃𝑎𝑟_𝑡𝑜𝑡	μM	1500	k𝐶𝑎,𝑃𝑎𝑟𝑜𝑓𝑓	ms^-1	0.0005
𝐶𝑠_𝑡𝑜𝑡	μM	31000	k𝐶𝑎,𝐴𝑇𝑃𝑜𝑛	μM^-1ms^-1	0.15
K_𝑆𝑅	μM	1	k𝐶𝑎,𝐴𝑇𝑃𝑜𝑓𝑓	ms^-1	30
υ_𝑆𝑅	μMmsμ-1m^-3	4.875	k𝑀𝑔,𝑃𝑎𝑟𝑜𝑛	μM^-1ms^-1	3.3 × 10^-5
L_e	μm^3ms^-1	2 × 10^-5	k𝑀𝑔,𝑃𝑎𝑟𝑜𝑓𝑓	ms^-1	0.003
τ_R=τR𝑆𝑅	μm^3ms^-1	0.75	k𝑀𝑔,𝐴𝑇𝑃𝑜𝑛	μM^-1ms^-1	1.5 × 10^-3
τ_𝐴𝑇𝑃	μm^3ms^-1	0.375	k𝑀𝑔,𝐴𝑇𝑃𝑜𝑓𝑓	ms^-1	0.15
τ_𝑀𝑔	μm^3ms^-1	1.5	k𝐶𝑎,𝐶𝑠𝑜𝑛	μM^-1ms^-1	4 × 10^-6
k𝐶𝑎,𝑃𝑎𝑟𝑜𝑛	μM^-1ms^-1	0.0417	k𝐶𝑎,𝐶𝑠𝑜𝑓𝑓	ms^-1	0.005
i	μm^3ms^-1	300

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

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

In the non-coincidence region () the function of combining and troponin, , is calculated according to Eq. (13):

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

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

Table 2

Model parameters

Parameter	Unit	Value	Parameter	Unit	Value
T_𝑡𝑜𝑡	μM	140	k2𝑜𝑛	ms^-1	0.15
k𝐶𝑎,T𝑜𝑛	μM^-1ms^-1	0.04425	k2𝑜𝑓𝑓	ms^-1	0.05
k𝐶𝑎,T𝑜𝑓𝑓	ms^-1	0.115	f_p	ms^-1	15
k0𝑜𝑛	ms^-1	0	f_0	ms^-1	1.5
k0𝑜𝑓𝑓	ms^-1	0.15	h_p	ms^-1	0.18
k1𝑜𝑛	ms^-1	0	h_0	ms^-1	0.24
k1𝑜𝑓𝑓	ms^-1	0.12	g_0	ms^-1	0.12

The 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 . When the product of solubilities of and in sarcoplasmic reticulum exceeds 6 mM, it is considered that the 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 in the sarcoplasmic reticulum to generate precipitation. A of 20 mM can reduce 29% of released by sarcoplasmic reticulum. When the muscle fiber contraction ends, in sarcoplasm can slowly be removed, and the sarcoplasm would then recover to a resting state of 3 mM.

Results and discussion

Analysis 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.

Table 3

Ion concentration in the myofibril at resting

Ions	Concentration	Unit
Calsequestrin (sarcoplasmic reticulum)	31000	μM
𝐶𝑎^2+ (sarcoplasmic reticulum)	1500	μM
𝐶𝑎^2+ (sarcoplasm)	0.05	μM
𝑀𝑔^2+ (sarcoplasm)	1000	μM
ATP (sarcoplasm)	8000	μM
Parvalbumin (sarcoplasm)	1500	μM
Troponin (sarcoplasm)	240	μM

Table 4

The percentage of ions in the binding site at resting

[height=0.8cm,width=3.4cm]IonsConjugates	𝐶𝑎^2+	𝑀𝑔^2+
Calsequestrin (sarcoplasmic reticulum)	54.55%	–
Troponin (sarcoplasm)	7.1%	–
Parvalbumin (sarcoplasm)	41%	54.1%
ATP (sarcoplasm)	–	90.9%

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 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 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 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 fluoresce in of 11 muscle fibers. When 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 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.

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 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, 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 and troponin. The simulation results of the binding of and troponin in region V4 are shown in Fig. 6. With the combination of 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 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 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 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.

Figure 7.

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

Simulation 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 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 , a 6-state model was used when and troponin combines. In region , due to the coincidence of thick and thin myofilaments, when troponin was combined with two ions, the cross-bridge was activated, so an 8-state model was used in this region.

When 1.1 m, the thick and thin myofilaments coincided completely. The conformations of troponin and tropomyosin regulated units are changed by the binding of 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 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 and troponin. This promoted the conformation change of tropomyosin and decreased the rate of separation of 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 8.

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

Figure 9.

Degree of cross-bridge activation.

Figure 10.

ATP hydrolyzed to produce phosphate.

When the length of the half-sarcomere is , 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 produced by ATP hydrolysis in the conditions of different degrees of coincidence of thick and thin myofilaments.