1. Introduction
Running has gained popularity as one of the most common and important forms of physical exercise. Currently, popular running methods can be divided into two types, running on the ground and running on a treadmill. Ground running is one of the simplest and earliest forms of running, requiring little extra equipment. But ground running can be affected by many environmental factors, such as road conditions and weather. Therefore, more and more people are choosing to run on treadmills that can allow them to run indoors.
Treadmills are frequently utilized in a variety of settings, such as physical therapy, exercise testing, rehabilitation, and research. However, as one of the most common running surfaces, the existing standard only defines the structural safety aspects of the treadmill while neglecting to specify and study the mechanical properties of the surface [
1,
2]. Similarly, there is no clarity on the mechanical features of treadmills in the numerous scientific studies comparing treadmills to running on the ground [
3,
4,
5]. In fact, the mechanical properties of the treadmill deck have an important effect on physiological response. It has been reported that athletes adjust leg stiffness and dynamics when running on surfaces with different mechanical properties [
6,
7]. During running, the body is continuously subjected to repetitive impact forces that are caused by GRF during the stance phase. Repeated impact force during exercise is the main cause of chronic injury and overuse injury [
8,
9]. Therefore, the characteristics of the sports surface itself are one of the important factors in the running process, because they are related to the performance and safety of the runner.
The potential effects of GRF and injuries during running have been widely discussed. A number of multi-degree-of-freedom models have been developed to study the biomechanics of human gait, given the intricate nature of the human body, which consists of multiple interconnected systems. Among them, the MSD model is the most commonly used method to study human gait characteristics. Kim [
10] and Derick [
11] have proposed some major MSD models that assess the GRF during the stance phase of human running. Günther et al. [
12] developed a multi-body model of the anthropomorphic lower limb that also takes into account the movement of soft tissue blocks relative to the bone to obtain more accurate estimates of bone kinematics and the distribution of internal forces in the leg. Liu and Nigg [
13] proposed a 4-degree-of-freedom MSD model called the LN model. This model is the first to comprehensively consider the rigid mass and wobbling mass of the upper and lower limbs of the human body. In addition, the shoe–foot nonlinear model is used to study the GRF in the stance phase of human running, so the LN model has received wide attention. In fact, GRF is a result of model calculations and is consistent with experimental results. Subsequently, Zadpoor et al. [
14] modified some of the parameters of the LN model to minimize the difference between the numerical and experimentally obtained GRF peaks.
While some researchers believed that stiffness and damping properties of shoes and sports surfaces alter GRF characteristics [
15,
16], nevertheless, research with expert runners has demonstrated that changing the hardness parameters of footwear and active surfaces has little impact on GRF [
17]. The GRF’s dominant frequency is relatively close to the natural frequency of LLST [
18], so during running the calf muscle may experience a resonance phenomenon, increasing soft tissue vibration amplitude. However, running does not exhibit such phenomena according to experimental studies [
19]. Experiments have shown that the stiffness of leg muscles varies with the environment, such as the hardness of shoes and sports surfaces. As a result of these physiological phenomena, some researchers have proposed a theory of muscle activity based on central nervous system (CNS) control [
20,
21,
22,
23]. It is also claimed that the CNS will keep the first peak (happening at heel strike) and the second peak (happening at the toe-off) of the GRF stable and will be able to modulate the viscoelasticity of the LLST to reduce the amplitude of vibration in the soft tissues of the lower limbs when the heel hits the ground.
However, the LN model, as a passive model, is unable to observe changes in soft tissue stiffness and damping in the leg. Therefore, for the purpose of predicting these biomechanical actions of the LLST during running, several active motion models have been developed. Zadpoor and Nikooyan [
24] proposed a new active model, the LNZN model, which is improved from the LN model, and the stiffness and damping coefficient of lower-limb muscles are adjustable. Based on the theory of muscle activity controlled by the CNS, they developed an optimal control strategy to modify the coefficients, a model of central nervous system function designed to achieve goals determined by two physiological assumptions: “constant force” and “constant vibration”.
In the past, studies of biomechanical models, including the LNZN model, have not explicitly considered the effect of the treadmill deck, whereas sports surface properties have a significant impact on the running process [
25]. Therefore, it is crucial to tie the GRF to the mechanical characteristics of the treadmill deck using a mathematical model. For this reason, in this paper we improve the LNZN model by considering the stiffness and damping parameters of the treadmill deck. The newly proposed model was used to explore the impact of treadmill running on muscle activity and GRF. Based on the theory of muscle activity regulation, the objective function is established. Using the particle swarm optimization algorithm, the optimal solution is determined as the safe range of the hardness parameters of the treadmill deck, and the muscle activity in this range is considered reasonable. This range should be considered as one of the factors considered in the safety area when designing the mechanical parameters of the treadmill deck. The overall architecture flow diagram of our work is shown in
Figure 1.
2. Materials and Methods
2.1. MSD Model
Among the crucial roles of the musculoskeletal system, movement stands out as a paramount function. There have been extensive experimental studies and modeling techniques discussed by scholars. While experimental methods hold value in elucidating multiple facets of locomotion, they are subject to certain constraints. For example, not all quantities can be easily measured, experiments cannot study isolated effects of parameters, and experiments might require access to numerous participants and specific testing conditions.
Zadpoor and Nikooyan improved the 4-degree-of-freedom active MSD model based on the LN model and provided a detailed description of the model parameters [
26]. The model presented in this paper is based on the LNZN model. We added a layer to the model that represents the subsystem of the treadmill deck. The changed model is shown in
Figure 2. There are five masses in this model. The upper and lower limbs’ rigid mass units are shown as
m1 and
m3, respectively, which are connected by a set of linear springs and dampers, denoted by
k1 and
c1, respectively.
m2 is the soft tissue mass unit of the lower limb, connected to the rigid mass
m3 of the upper limb by a spring with a stiffness coefficient of
k3, and connected to the rigid mass
m1 of the lower limb by a set of variable stiffness
k2 and damping
c2.
m4 is the soft tissue mass of the upper limb connected to the rigid mass of the upper limb by a spring with a stiffness coefficient of
k5 and a set of parallel linear stiffness (
k4) and damping (
c4), where rigid mass represents the mass of the human skeleton and wobbling mass represents the mass of the human body muscle and blood, etc.
m5,
kt, and
ζt represent the mass, stiffness, and damping ratio of the treadmill deck, respectively.
The default model parameter values utilized in this paper’s numerical simulation are consistent with those provided in the LNZN model, summarized and proposed by Zadpoor and Nikooyan [
24], and the specific values are shown in
Table 1. The initial displacement of all mass blocks is 0.
According to the modified model, the equation governing motion is derived as follows:
where the coefficients
m,
c, and
k, respectively, reflect mass, damping, and stiffness in
Figure 2, and
g denotes the gravitational acceleration. Furthermore, following the proposal of Nikooyan and Zadpoor [
26], the magnitude of GRF is represented by
Fg and we modify the vertical contact force
Fg according to the modified model as:
In effect, this is a nonlinear viscoelastic model of
Fg, where the parameters
a,
b,
c,
d, and
e are shoe-specific parameters. In the LNZN model, Zadpoor and Nikooyan [
24] proposed the parameter values of soft shoes and hard shoes, as shown in
Table 2. The same values are used in the simulation in this article.
2.2. Determine the Objective Function
According to biomechanical hypotheses, the CNS is said to be able to maintain constant impact peak and activity peak by altering muscle activity in the LLST [
27]. On the other hand, under typical working conditions, the peak value of GRF is affected by the treadmill deck parameters
kt and
ζt. Thus, the CNS may be viewed as a controller that, in order to satisfy the biomechanical assumptions, minimizes the objective function for each introduced pair of treadmill deck hardness parameters to determine the best tuning parameters (
k2 and
c2). The objective function described below, proposed by Zadpoor and Nikooyan, seeks to reduce the disparity between peak values of the GRF and the amplitude of leg soft tissue vibrations from their respective experimental constants during running.
P1 and
P2 refer to the values of the first and second peaks in the GRF waveform, correspondingly. The default settings for GRF peaks are
P10 and
P20 (
P10 = 1436.8 N,
P20 = 2026.4 N), which are determined according to specific parameters in the LNZN model. Λ represents the magnitude of oscillation in the lower-limb mass, which can be determined by fitting a sine curve to the
x2 graph. Nonetheless, research indicates that the oscillations of the wobbling mass in the lower limb do not conform to a sine wave approximation [
24]. Thus, the maximum
x2 value is regarded as the amplitude, and the default amplitude is Λ
0 = 10.1 cm.
In order to normalize the objective function, we need to determine the P1, P2, and Λ values determined after optimisation for the different treadmill deck parameters and the difference between their corresponding default values, and then divide by the default values. As a result, different boundaries for the regulation parameters k2 and c2 are explored to obtain the outcomes presented in this paper.
2.3. Determine the Boundaries and Constraints
Proficient runners can regulate the peaks of GRF in their initial running gait step, and it can be assumed that muscle activity prepares the musculoskeletal system for ground impact prior to landing. In addition, it has been shown in the literature that the stiffness of leg muscles can change by up to 3.6 times and the damping properties can change by up to 4 times compared to rest [
20,
28], depending on the physiological capacity of the subject. Based on this fact, it is assumed that the interval of muscle activity adjustment will change and contract with muscle fatigue. Therefore, on the basis of the experiment, muscle fatigue is considered to be the weakening of the human body’s ability to regulate the soft tissues of the lower extremities, and Zadpoor proposed Formula (4) to narrow the interval [
29].
where
t is the time that has elapsed since the beginning.
ξ (
t) is the parameter that determines how much the interval shrinks.
ξmin and
ξmax are the minimum and maximum values of
ξ, respectively, set to 1/4 and 4. At each moment, the optimization boundaries for the maximum (
UB) and minimum (
LB) of the stiffness (
k2) and damping (
c2) parameters are determined by:
where
k20 and
c20 are the default values of
k2 and
c2, respectively, as shown in
Table 1 for the typical overdamped system of a treadmill. In this study, we considered the definition interval of the stiffness and damping ratio parameters
kt and
ζt of suitable treadmill decks according to the common treadmill performance in the market, and investigated them within the definition domain given below. For each pair of
kti and
ζti increasing within a given range, we search for the optimal value of the tuning parameter.
Therefore, for each pair of kti and ζti, we need to solve an optimization problem, which will lead to a large number of iterative solving optimization processes. Broadly, two categories of optimization methods exist. The gradient-based strategy, for example, begins the search with an initial estimate, and the reliance on the original guess may lead to an accuracy decrease by slipping into local minima. In this paper, we use another well-known optimisation method, particle swarm.
2.4. Simulation
In this section, we use numerical methods to solve the optimization problem defined in the previous section and simulate two fatigue conditions to simulate the function of muscle to regulate stiffness and damping.
The k2 and c2 values that minimize the objective function (Formula (3)) are determined for each set of suggested treadmill deck parameters (Formula (6)).
According to the given initial conditions, the state variables are determined by solving the equations of motion. The value of the objective is evaluated by finding the most suitable muscle stiffness and damping values to calculate the ground reaction force and its peaks, along with the amplitude of the LLST oscillation. For the different boundary limits defined for k2 and c2, the iteration is repeated using optimization techniques to look for the best k2 and c2 values. Similarly, these iterative optimization processes are replicated for each pair of treadmill deck parameters kti and ζti within the appropriate range.
Ultimately, the optimal values of the objective function, the optimal values of k2 and c2, as well as the optimal values of the components of the objective function, including the first and second peaks of the GRF and the amplitude of the vibrations, are plotted separately in the space of the treadmill deck hardness parameter kti–ζti. The region with the smallest optimal value for each pair of kti and ζti identified in the graph is defined as the safe region. Minimal muscle activity is expected if the model’s conditions fall within a safe range.
3. Results
Similar to the results of Zadpoor and Nikooyan [
29], there is a specific region in the
kti–ζti plane where the optimal value of the objective function is close to zero. The region where the optimal value of the objective function is less than 0.07 is called the safe region, as shown in
Figure 3a,b, and is identified by two dotted lines in each subgraph of
Figure 3. It can be seen that with the generation of fatigue, the size of the safe area is obviously reduced. This occurred due to the reduction in muscle capability while modifying its viscoelastic properties to align with the two underlying physiological hypotheses.
In
Figure 3c–f, the adjustment parameters
k2 and
c2 of lower-limb stiffness and damping are shown, respectively. It can be observed that higher damping and stiffness coefficients of the treadmill deck increase the lower-limb stiffness parameter
k2 and decrease the lower-limb damping parameter
c2. For the pre-fatigue state or post-fatigue state, the value of
k2 does not reach the upper bound, which means that the leg muscles do not need to be very stiff in order to minimize the objective function. We note that the stiffness and damping values of the lower-limb muscles also increase as fatigue arises within the safe area. This finding is consistent with the findings of Khassetarash et al. [
30], who observed that an increase in damping coefficient is observed when runners tire during running.
In
Figure 3g–j, the GRF’s first and second peak values acquired through each optimization are depicted. These outcomes show that by reducing the hardness of the treadmill deck, the first peak is decreasing, and conversely the second peak is increasing. However, within the safe area, there was no significant change in the values of the first and second peaks. This shows that when the hardness parameters of the treadmill deck are within a safe range, the CNS can keep the GRF peak stable, reducing the impact of shoe and movement surface changes. This is also consistent with the findings of Clarke et al. [
17] and Ferris and Farley [
28], who determined that GRF peaks remain constant over the appropriate range of parameters.
The vibration amplitudes that may be noticed on the calf muscles during the landing support phase are depicted in
Figure 3k,l. Within the hardness range of the treadmill deck, the vibration amplitude changes only very slightly (less than 5%) and is almost constant. This is because leg muscles have high stiffness and damping. This fact should also be attributed to the CNS, although it contributes less to CNS function than keeping GRF spikes constant.
4. Discussion
In order to evaluate the effectiveness of the objective optimization function, we will compare the simulation results with the experimental data of previous scholars. First, we will compare the simulated GRF calculated in this paper with the GRF values measured by Clarke et al. [
17] in actual running and the results of Ferris and Farley’s research [
28]. In both the simulation results and the experimental results, there is a range of hardness (safe area) within which the GRF remains constant. Outside the safe area, the GRF will deviate from the default value. Dixon [
8] and Kerdok [
31] have also observed that by varying the hardness of the sports surface or shoe, the ground reaction force does not change significantly within a certain range. Second, we compared the model predictions (
Figure 3k,l) with Nigg and Wakeling’s experimental observations of vibration levels [
20]. We can see that consistent with experimental observations, the simulated vibration amplitude changes only insignificantly, even if the hardness of the sports surface changes considerably. Thirdly, we compared the simulated stiffness and damping coefficients of lower-limb soft tissue with the previous experimental results. Previous experimental tests by Wakeling and Wilson [
32] have shown that damping plays a greater role than stiffness in muscle tuning. According to these studies, the change of stiffness coefficient is much smaller than that of damping coefficient. We also observed the same trend in our simulation results (see
Figure 3e,f).
Compared with the results of the LNZN model, our model shows that the maximum impact peak value is about 1500 N, which is only slightly higher than the expected impact peak value. Our results show that the GRF peak value does not change much over the whole range of hardness parameters. It is also worth noting that our results show that the size of the safe area is significantly reduced after exercise fatigue. This suggests that the design of the active shock absorption system of the treadmill needs to consider this point. Especially for runners with exercise fatigue, our results suggest increasing the hardness of the treadmill deck a little to provide more impact energy recovery, thus providing better protection and comfort for runners.
When compared to 3D human motion models [
33], the model utilized here has certain advantages and some drawbacks. The 3D model contains critical information such as joint angle and joint torque, which are required for both GRF and fatigue simulation. However, in the 3D model, the human body is simplified as a system of rigid masses, disregarding the influence of wobbling masses and thus not accounting for the dynamic stiffness of muscles. In addition, both the 3D model and the previous MSD human motion model did not take into account the influence of treadmill deck parameters, which is a component of the model employed in our investigation. Also, the MSD modelling approach gives us better physical insight into the phenomena.
The MSD model used in this paper is just a simplified model. The model mainly considers the vertical GRF of a standard human; however, the horizontal GRF may also have an important effect on human movement. Furthermore, within this model, the variation of stiffness and damping characteristics among individual muscles is overlooked, and the established safe area is derived from simulations. Each person’s central nervous system may exhibit different behaviors. We must also recognize that human muscle movement is an extremely complex process that is influenced by many factors, such as physical activity level, age, gender, etc. These individualized factors can lead to differences in muscle activity between individuals, and it is clear that the models considered in this paper have not been able to incorporate these factors. Therefore, a more realistic exercise model can be established in future work to consider these diversity factors more comprehensively, and these factors can be incorporated into the mechanical parameter optimization of the treadmill.
5. Conclusions
In this paper, by improving the original MSD human motion model and introducing the stiffness and damping of the treadmill deck, the optimum range of the running table parameters is discussed from the angle of the theory of muscle activity regulation. The ability of the muscle to move is considered to be the boundary of the range of muscle viscoelastic parameters used to simulate pre- and post-fatigue states. The results show that: (1) the mechanical parameters of the treadmill deck play a key role in the process of muscle activity regulation, there exists such a safe area provided by a theory based on muscle activity; (2) there is little variation in the peak values of GRF in the safe area, however, it is worth noting that the safe area post-fatigue is less than 1/3 of the pre-fatigue area; (3) the damping ratio is recommended to have an inverse relationship with stiffness when designing a treadmill and its stiffness is not recommended to exceed 650 kn/m. Considering the cost and time of experimental testing, manufacturers can refer to the results of this study when designing treadmills, either by adopting a specific hardness or by introducing treadmills with adjustable deck hardness to meet the physiological needs of different exercisers and reduce the risk of sports injuries.
Through this study, we have a deeper understanding of the design of treadmill deck and the optimization of human motion model. Later, this research could be extended to consider more factors and parameters to explore the treadmill deck parameters’ accuracy further. In general, this paper provides insights for understanding the influence of muscle activity regulation on treadmill deck parameters and provides a useful reference for the optimal design of treadmill deck and the improvement of human movement model.