1. Introduction
Skeletal muscle, the impetus of human movement, is an important part of the human motor system. It generates active contractions, leading to a variety of complex and subtle movements. The biomechanics of the skeletal muscle have attracted increased attentions due to their role in regulating human movements, improving competitive sports performance, preventing sports injury, and designing effective training programs for rehabilitation.
Muscle fiber, composed of paralleled myofibrils with average diameters of approximately 1 μm, is the main component of individual skeletal muscle [
1]. Myofibrils possess distinct patterns alternating between dark A-bands and light I-bands. The center line of light bands is defined as the Z-line, and the section between neighboring Z-lines is referred to as sarcomere, the fundamental unit of skeletal muscle. From the microscopic point of view, sarcomere contains thick and thin filaments. While the thin filament is composed of actin, tropomyosin-Tm, and troponin-Tn, the thick filament is formed by myosin molecular motors, which are key components of skeletal muscle contractions. Each thick filament is surrounded by six thin filaments and connected to the Z-line through an elastic element made of titin, whose passive elasticity plays an important role in muscle contraction.
Biomechanics research on skeletal muscles has been focused on the macro and micro levels, separately [
2,
3,
4,
5]. At the macro level, the Hill’s model, based on
in vitro muscle experiments, was widely adopted. However, this macro phenomenal model neglected the linkage to the microscopic characteristics of the skeletal muscle. At the micro level, Huxley proposed the sliding-filament theory which described that the muscle contraction resulted from of the relative slip between the thick and thin filaments, induced by myosin molecular motors. However, in his research, the description for the cycling working process of myosin motors is insufficient. Myosin motor is one kind of protein that transfers chemical energy from catalytic hydrolysis ATP (Adenosine Triphosphate) into mechanical work,
i.e., directional movement [
6,
7,
8,
9]. The mechanochemical coupling of myosin motors is then used to understand the micro-mechanical behavior of skeletal muscle.
Lymn and Taylor first described myosin motor’s work cycle, considering four steps, in 1971 [
8,
10]. Julicher
et al. [
11] then proposed a two-state brown ratchet model for a single motor and obtained the collective properties of molecular motors. This model has been extended to different families of molecular motors, such as kinesin, dynein, and myosin. Chin
et al. [
12] divided the biochemical cycle of myosin motor into seven states, and proposed a simulation model that allows altering the rates of transitions, mimicking the role of chemicals, including ATP, ADP,
etc., between any of the states in a stochastic cycle. Wang [
13] developed a mathematical model for investigating the mechanochemical coupling mechanism of molecular motors. However, these models lack the quantitative linkage to the bulk behaviors of muscle contractions.
Lecarpentier
et al. [
14] related the macro properties of skeletal muscle, such as maximum unloaded shortening velocity, peak force, the cross-sectional area of a muscle strip with the micro characteristics, including the force per single cross-bridge, probability distribution between each working state, and peak muscle efficiency. The statistical mechanics combined with Huxley’s formula was used in their study.
The goal of this study is to develop a multiscale mathematical model for quantifying the mechanical behaviors of the active contraction skeletal muscle in relation with the mechanochemical coupling of molecular motors. First, the four states of the working mechanisms of myosin molecular motors in muscle contraction are analyzed. The evolution process of myosin molecular motors, regulated by load, chemicals, and Brownian force, was described. Second, the Fokker-Planck Equation was used to describe the collective behavior of myosin molecular motors. Third, the active contraction force and velocity of muscle fiber are calculated using the non-equilibrium statistical mechanics; and, finally, the biomechanical properties, at multiple scales, include probability density distribution of myosin motors and muscle fiber contraction force and velocity, are quantitatively related. Our results were then compared with published experiments and predicted models.
2. Biomechanical Model of Skeletal Muscle Contraction
The biomechanical model of skeletal muscle is generally characterized by parameters that reflect the physical structure of fiber (e.g., cross-sectional area, fiber length) and the behavior of molecular motors (e.g., transition rates between different chemical state, the free energy barriers established in the chemical reaction). The magnitudes of these parameters correspond to various fiber behaviors, such as fast and slow fibers. Regardless of these alternations, the mechanochemical coupling induced contraction force is the common feature among skeletal muscle.
2.1. Working Mechanism of Myosin Molecular Motor
In general, the muscle-contraction process consists of three steps: the muscle cell is excited, filaments slip and muscle contracts, and muscle stretches. Relative slip between the thick and thin filament is the direct cause of muscle contractions. The slip driving force is generated by the interaction between myosin motor and actin. The Ca
2+ concentration in sarcomere increases suddenly under nerve impulse excitation, then, the binding sites on thin actin filaments are opened and the cross-bridge of myosin can blind with actin. Then, a power-stroke starts through the hydrolysis of ATP and causes the slipping between thin and thick filaments [
4]. A large number of myosin motors acting on the thin filaments leads to the contraction of sarcomere.
The cycle of each myosin motor, first described by Lymn and Taylor [
8,
10], was divided into four main steps, also referred to as strokes,
i.e., Detachment-stroke, when myosin motor detaches from actin with the new attachment of ATP; Recovery-stroke, when myosin moves forward to be ready for reattachment to the actin; Attachment-stroke, when myosin adheres to the actin; and Power stroke, when the myosin moves backward and forms a strong binding to the actin and the stored chemical energy released through ATP hydrolysis. One working cycle is completed by four chemical states, coupled with conformational changes, and the movement of a myosin motor, from one binding site to the next binding site, is then achieved.
2.2. Collective Dynamic Model of Myosin Molecular Motors in Sarcomere
A myosin motor has many degrees of freedom associated with the unidirectional motion of the motor. The one-dimensional motion along the filament’s slip direction is considered in this work. The distance between two binding sites on a thin filament is designated as L. A working cycle of a myosin motor is divided into N = 4 steps, each step corresponding to one chemical state. σi is used to represents the state i. Myosin motor transitions stochastically between different states, under the effect of existing external energy (for example, released from ATP hydrolysis), and a myosin motor can move directionally and push the thin filaments. Let x be the variable as the displacement of a myosin motor along the thin filament, due the conformational change of the cross-bridge, is the relative slipping distance between thin and thick filaments. Let ρ(x,t) be the probability density that the myosin motor, at the position of x in the occupancy state at time t, possesses.
As the size and mass of molecular motors are relatively small, their stochastic motion follows the law of Brownian motion, which was well described by the Langevin Equation, regardless of inertia. In addition, the change of chemical state occupancy is depicted by a discrete Markov process. The combination of the two aforementioned governing equations then regulate the stochastic evolution of molecular motors.
For the molecular motor alone, the mechanochemical coupling mechanism is the function of multi-force interactions, including the motor force, the external load, the Brownian force obeying the fluctuation dissipation theorem, and the active force produced by chemical reactions. Its stochastic evolution of the probability density is represented by the Fokker-Planck Equation [
14,
15,
16], as follows:
where,
D =
kBT/γ
m, is the diffusion coefficient of myosin motor in aqueous solution,
kB is the Boltzmann constant,
T is the absolute temperature, and γ
m is a drag coefficient. The value
Vi(
x) is the active potential energy produced by chemical reaction in state
i, and is a function of the geometrical coordinate
x in the chemical occupancy state
i.
FLoad is the external load force on the myosin motor and
kij(
x) is the transition rate from chemical occupancy state σ
i to state σ
j, the value depends on the potential energy in the two states. In addition, the probability density equation should satisfy the constraint defined by:
2.3. Solving Fokker-Planck Equation
The continuity equation (Equation (1)) is discretized as a spatially discrete Markov chain [
16] to calculate the probability density distribution of myosin motors. A jump progress can then be used to estimate the continuous motion of myosin molecular motors, as demonstrated in
Figure 1. The computational region [0,
L] is divided into
M sub-intervals, so the interval between two sites is Δ
x =
L/
M. Let
and
indicate the forward and backward jump rates between sites
xn–1 and
xn. The jump rates can be calculated as follows:
where,
and
.
The probability of the myosin motor between [
xn,
xn+1] is calculated as:
Therefore, Equation (1) can be approximated by the master equation for the jump process:
Figure 1.
The spatial discretization and jump rates between two sites.
Figure 1.
The spatial discretization and jump rates between two sites.
2.4. Model Simplification
It is challenging to solve the above equations due to the four steps in one complete working cycle of myosin motor. Considering that myosin motors only act in the case of a strong binding to the actin, thus, the working cycle of myosin motor can be simplified to two strokes: attachment-stroke and detachment-stroke, as depicted in
Figure 2. Correspondingly, Equation (6) can be simplified as:
where
p1 and
p2 are the probabilities in the two simplified states. The transition rates between the two states (
k12 and
k21) are calculated as:
where [ATP] represents the concentration of ATP in the chemical reaction. The transition rate is influenced by the load force in the power stroke state, while it depends on free energy in other states. The relationship between transition rate and load force can be expressed as:
where θ
+ and θ
− are the distribution coefficients that indicate the effect on the forward and backward transition rates, respectively [
17].
d is the moving distance of the molecular motor on the actin filament. The other transition rates
and
k−i can be calculated with the detailed balance requirements [
14].
Figure 2.
The simplification of working cycle of myosin motor.
Figure 2.
The simplification of working cycle of myosin motor.
Equation (7) could be written in a matrix form as:
The probability distribution subject to the constraints:
where,
The steady-state solution of the probability density distribution of molecular motors can be obtained by making the left side of Equation (12) equal to zero. In addition, according to the asymmetrical periodicity of potential, the two potentials are given by:
where Δ
G indicates the energy derived from one ATP hydrolysis. The potentials determine the values of jump rates
Fin and
Bin in Equation (12). The jump rates can be calculated with Equations (3) and (4). Equation (15) is used to represent the effect of the chemical reaction on the myosin motion, and link the chemical reaction and motor motion.
2.5. Active Contraction Force and Velocity of Muscle Fiber
Muscle fiber is composed of sarcomeres that are arranged in series and in parallel. According to mechanical principles, the muscle active contraction force is proportional to the number of sarcomeres in parallel. That is to say, the muscle active force is proportional to the cross-sectional area of muscle, once the force generated by each sarcomere is determined. In this work, we hypothesized that all sarcomeres’ working processes are synchronized and that the Ca
2+ concentration is large enough to detach all the binding sites on thin filaments during movement. Then, the probability distribution of myosin motors at any time can be used to calculate the magnitude of force acting on thin filament. The contraction force produced by all myosin motors in muscle fiber can be obtained as:
where,
is the first-order of displacement, which also means the average displacement along the thin filament of all myosin motors,
n0 is the number of myosin motors in one thick filament,
ke is the elasticity of cross-bridge,
A is the cross-sectional area of muscle fiber,
s is the area of hexagonal surrounded by 6 thin filaments around a thick filament, α is the degree of overlap between thin and thick filaments, and as the sliding-filament theory confirms, there exists an optimal overlap to generate a maximum force [
5], a Gaussian function is used to imitate the process in this paper:
where,
lopti indicates the optimal length of muscle fiber.
The muscle fiber contraction velocity of the sarcomere can be calculated by the evolution of the molecular motor as:
where:
where the probability flux and average drift velocity of molecular motor in state
i was defined through Equations (1) and (7) [
16]. For a muscle fiber that is composed of sarcomeres, the contraction velocity will be:
where,
ms is the number of sarcomere that one myofibril contains in series.
4. Conclusions
In this work, a collective dynamic model of myosin molecular motors is developed, based on both the Langevin Equation and Markov process. Its behaviors are then coupled with the active contraction force and contraction velocity of skeletal muscle at the macroscopic scale. The roles of motor potential profile, ATP concentration, and load force on myosin motors are investigated Results are summarized as: (i) The probability of the attachment state in the two-step working cycle increases with a larger load force of motors; (ii) the probability density in the attachment state increases with [ATP] but reaches a saturation state at a certain [ATP]; (iii) the obtained force-length and velocity-force properties of muscle are validated against experimental observations. This indicated that our multiscale model could capture the major features of muscle contractions; and (iv) the proposed multiscale model in this work could be used for characterizing any muscle contraction.
In the present model, the force-length relationship was assumed as symmetric regarding optimal muscle length. Its asymmetry could be implemented in the future work using a modified Gaussian function. The sarcomere was also assumed to be fully activated, i.e., the Ca2+ concentration is high enough for myosin motors to bind strongly with the actin. However, the Ca2+ concentration is regulated by the firing rate of nerve impulses. Despite these simplifications, the present work demonstrated the importance of multiscale modeling for a better understanding the mechanism of skeletal muscle active contractions. Quantitative relation between myosin molecular motors and muscle contractions could be used to provide guidance for designing effective rehabilitation exercises, improving performance in competitive sports, and exploiting the potential to prevent muscle injuries.