Improved Hill-type musculotendon models with activation-force-length coupling

Abstract

BACKGROUND:

Hill-type musculotendon models are most commonly used in biomechanical simulations for their computational efficacy and efficiency. But these models are generally built for maximally-activated muscles and linearly scale muscle properties when applied to submaximal conditions. However, the precondition of this scaling, which is muscle activation and properties are independent each other, has been proven unreal in many studies. Actually, the maximal activation condition is not ubiquitous for muscles in vivo, so it is necessary to adapt the linear scaling approach to improve the model practicability.

OBJECTIVE:

This paper aimed at proposing two improved Hill-type musculotendon models that are better suited for submaximal conditions.

METHOD:

These two models were built by including the activation-force-length coupling and their biological accuracy and computation speed were evaluated by a series of benchmark simulations.

RESULTS:

Compared to experimental measurements, the percent root mean square errors of forces calculated by the two AFLC models were less than 13.98% and 13.81% respectively. However, the average running time of the second AFLC model was nearly 17 times that of the first one with only a little improvement in accuracy.

CONCLUSION:

The two AFLC models were validated more accurate than the common Hill-type model in submaximally activated conditions and the first one was recommended in the construction of upper-layer musculoskeletal models.

Keywords

Hill-type model musculotendon model activation force-length coupling submaximal activation

1. Introduction

The individual muscle force is a reflection of neural control and tissue loading and its accurate estimation can provide better diagnosis and treatment for neurological and orthopaedic diseases [1]. Unfortunately, direct measurement of muscle forces in vivo is not feasible in clinical setting. Therefore a mathematical musculotendon model, which takes anatomical and physiological characteristics of muscles into account is required. Generally, there are two categories of musculotendon models: one is the Huxley-type molecular model and the other is the Hill-type phenomenological model. In muscle-driven simulations of human and animal movements, Hill-type models are most widely used because they are simple in computation and make good use of commonly measured experimental variables [2].

The Hill-type muscle models were named after English physiologist Hill [3] who measured the heat produced by muscles during contractions in quick-release experiments on frog’s sartorius in 1938 and modeled muscles with three elements, which was a contractile element (CE) and two non-linear spring element, one in series (SEE) and the other in parallel (PE). A detailed introduction to the Hill-type model could be found in Zahalak [4] and Winter [5].

However, Hill [3] did not delve into how the tendon elasticity affected muscle forces. Later, Zajac [6] assessed the role of tendon in movement and proposed a one-parameter mathematical model, which took muscle and tendon as a whole, i.e., the musculotendon actuator. This model has become the foundation of almost all Hill-type musculotendon models used in biomechanical simulations. However, it will suffer a singularity problem in some situations, such as when function of the force-velocity relationship is not invertible or muscle activation is too small, which are common during normal movement conditions.

To obtain a stable musculotendon model, it is requisite to solve the singularity problem. Millard et al. [7] developed three non-singularity models: an equilibrium model, a damped equilibrium model and a rigid tendon model by applying boundary constraints to some variables or modifying the mathematical equation with some assumptions. However, in submaximal conditions, they revealed larger errors as they adopted the linear scaling method, which scaled the muscle properties linearly with activations and didn’t take their coupling into account.

Actually, this linear scaling approach should have the prerequisite that muscle activation and properties are independent each other, which has been proven unreal in many studies. Experimental validations of these models by Perreault et al. [2] revealed that the model errors during movements were largest at low motor unit firing rates, which were relevant to normal movement conditions. And Lloyd and Besier [8] employed the coupling between activation and force-length property and achieved an improvement in model accuracies in the simulations of an EMG-driven musculoskeletal model, which estimated muscle forces and knee joint moments in vivo. Actually, muscles in vivo are often unlikely to be maximally activated during daily life [9], so it is necessary to adapt the scaling approach to improve the practicability of Hill-type musculotendon models.

The main purpose of this paper is to propose two Hill-type musculotendon models with activation-force-length coupling (AFLC), which is significant in submaximal activation situations and validate their biological accuracy and computation speed in submaximal activation situations using benchmark experiment data. In Section 2, base component of the two models will be described first, then submaximal activation conditions will be analyzed and finally two methods of incorporating activation-force-length coupling will be introduced. In Section 3, benchmark simulations will be described and experimental results will be analyzed and discussed.

Figure 1.

Base component of AFLC models, including an active contraction element (CE), a passive elastic element (PE), a parallel damp element and a serial elastic tendon.

2. Activation-force-length coupling (AFLC) models

2.1 Base component of the two models

The base component of our models consists of an active contractile element (CE), a passive elastic element (PE), a parallel damping element and a serial elastic tendon, as shown in Fig. 1. Forces generated by the four elements are represented by an activation-varied force-length-velocity relationship $af_{L}(a,\tilde{l}_{m})f_{V}(\tilde{v}_{m})$ , a passive force-relationship $f_{P}(\tilde{l}_{m})$ , a damping force $\beta\tilde{v}_{m}$ and a tendon force-length relationship $f_{T}(\tilde{l}_{t})$ , respectively. Assuming muscle mass is negligible, the forces of muscle and tendon should be in equilibrium, then the musculotendon model can be represented by Eq. (1).

$\displaystyle F^{\max}\left[af_{L}(a,\tilde{l}_{m})f_{V}(\tilde{v}_{m})+f_{P}(% \tilde{l}_{m})+\beta\tilde{v}_{m}\right]\cos\left(\alpha(\tilde{l}_{m})\right)% -F^{\max}f_{T}(\tilde{l}_{t})=0$ (1)

where the tilde means normalized dimensionless parameters, all the forces are normalized by maximum isometric force $F^{max}$ , $\tilde{l}_{m}$ , $\tilde{v}_{m}$ and $\tilde{l}_{t}$ represent muscle length, muscle velocity and tendon length normalized by optimal fiber length, maximum velocity and tendon slack length, respectively. The four musculotendon property curves, $f_{L}(1,\tilde{l}_{m})$ , $f_{V}(\tilde{v}_{m})$ , $f_{P}(\tilde{l}_{m})$ , $f_{T}(\tilde{l}_{t})$ are represented using quantic Bezier splines [10] and can be adjusted to different types of muscle and tendon [7]. We adopt the default force curves developed by Millard et al. [7] who curve fitted experimental data about the force-length, force-velocity, passive force-length and force-tendon relationship. $\beta$ is damping coefficient, whose default value is 0.1, $a$ is muscle activation ranging from $a_{\min}$ to 1, and $\alpha(\tilde{l}_{m})$ is length-varied pennation angle.

We assume that the muscle belly has a constant thickness and volume when contracting, so the pennation angle $\alpha(\tilde{l}_{m})$ will change with fiber length to maintain a constant height [7]. Equation (2) is used to calculate pennation angle at different fiber lengths.

$\displaystyle\alpha(\tilde{l}_{m})=\sin^{-1}(\sin\alpha_{o}/\tilde{l}_{m}),% \quad\alpha\in[0,\pi/2)$ (2)

where $\alpha_{o}$ is pennation angle at optimal fiber length.

Equation (1) is actually a nonlinear ordinary differential equation (ODE) taking muscle fiber length $\tilde{l}_{m}$ as an argument and can be solved by a Runge-Kutta-Merson [11] algorithm. Besides, a unilateral constraint is employed to avoid unrealistically short muscle length, i.e.,

$\displaystyle\tilde{v}_{m}=\left\{\begin{array}[]{ll}0&\text{if $\tilde{l}_{m}% \leqslant\tilde{l}_{m}^{\min}$ and $\tilde{v}_{m^{\ast}}<0$}\\ \tilde{v}_{m^{\ast}}&\text{otherwise}\\ \end{array}\right.$ (3)

where $\tilde{v}_{m^{\ast}}$ is a candidate value for $\tilde{v}_{m}$ computed by Eq. (1).

It is worth mentioning that AFLC models are general versions and can be scaled to different types of muscles by maximum isometric force, optimal fiber length, tendon slack length, pennation angle at optimal fiber length, which can be measured or found in literature for different muscles.

2.2 Submaximal activation analysis

Huijing [12] indicated that the most limiting factor of the phenomenological models was that they were designed to describe effects in maximally activated muscles only. And the reason why submaximally activated situations are not paid enough attentions is that most experiments describing muscle properties were performed under conditions of supramaximal stimulation of the peripheral nerve. But muscles are not likely to be maximally active very often during daily activities, so the model has to solve this problem before being applied to the human musculoskeletal modelling.

The general way of estimating muscle forces at submaximally activated conditions is to scale those of the maximally activated ones linearly, so normalized force-length curves of maximally and submaximally activated muscles are identical. This linear scaling approach is under two hypotheses: (a) The relationship between the normalized force and activation is linear. (b) The muscle length is independent of activation. But unfortunately neither of these two hypotheses is true. As for a fully recruited muscle, submaximal activation means manipulating the firing rate, then the first hypothesis is equivalent to a linear force-firing rate relationship. Skinned muscles can be directly activated by the external application of Ca ${}^{2+}$ buffered solutions and changes in intracellular Ca ${}^{2+}$ concentration can be monitored by suitable Ca ${}^{2+}$ indicators. Therefore, the intracellular equivalent relationship of force-firing rate can be achieved by experiments on skinned muscle fibers as the force-pCa ( $-\log_{10}$ [Ca ${}^{2+}$ ]) curve. Stephenson and Wendt [13] presented a few isometric force-pCa curves obtained with different types of mechanically skinned skeletal muscle fibers. These curves are clearly sigmoidal in shape, which breaks the first hypothesis. And as the apparent sensitivity of the contractile proteins to Ca ${}^{2+}$ depends on muscle length, the shape of force-pCa curve varies with muscle lengths, then the second hypothesis fails.

When coming back to the force-length curve, the direct observation is the plateau shifting toward longer lengths as the level of activation decreases. Experiments on both sarcomeres and intact fibers of different types of muscles showed the similar results. In addition, Huijing [12] indicated that both active slack length and optimum fiber length shifted substantially to higher muscle length at lower constant frequencies and the shifting distance of the former is smaller than the latter. And the maximal length of force exertion was considered a fixed property of muscle regardless of its degree of activation.

Figure 2.

Force-length curves at different activations in the first AFLC.

2.3 Two methods of incorporating activation-force-length coupling

Actually, Huijing [12] only qualitatively described the coupling between muscle activation and force-length relationship and appealed to muscle modellers to take this property into account. Hatze [14] incorporated effects of firing frequency in his model by altering the slope of the force-frequency curve as a function of muscle length, but did not fully succeed for ignoring shifts of active slack length. Van Zandwijk et al. [15] modified Hatze’s model and excelled in predicting muscle forces by optimizing the model for an isometric twitch. However, both the models took Ca ${}^{2+}$ concentration as one input, which was only measured in special experiments meaning that the models maybe not appropriate to integrate to the musculoskeletal models for human movement analysis.

When constructing an EMG-driven musculoskeletal model of the human knee, Lloyd and Besier [8] built a muscle contraction dynamics model, which incorporated the coupling between activation and optimal fiber length using Eq. (4).

$\displaystyle l_{m}^{o}(t)=l_{m}^{o}(\gamma(1-a(t))+1)$ (4)

where $\gamma$ is the percentage change in optimal fiber length, $l_{m}^{o}$ is optimal fiber length at maximum activation, $a(t)$ and $l_{m}^{o}(t)$ are activation and optimal fiber length at time $t$ , respectively. This equation changes the force-length curve indirectly and shifts the whole curve to larger lengths. As shown in Fig. 2, Lloyd and Besier [8] examined the effect of the physiological parameter in a contrast test with the knee model, which predicted joint moments in vivo and revealed that the joint moment $R^{2}$ significantly increased 0.6 by including $\gamma$ , whose default value was 0.15.

Figure 3.

Activations at constant stimulating frequencies of 10, 20, 30 Hz and activations at random stimulating frequencies, of which mean values are 10, 20, 30 Hz.

In consideration of Lloyd’s satisfactory results, two possible methods have been designed, leading to two activation-force-length coupling (AFLC) models. As the variables in force-length relationship are both normalized by optimal values, the first method is to adjust the optimal fiber length to activations. In that way, the relative muscle length $\tilde{l}_{m}$ at time $t$ can be obtained by Eq. (5),

$\displaystyle\tilde{l}_{m}(t)=\frac{l_{m}(t)}{l_{m}^{o}(\gamma(1-a(t))+1)}$ (5)

where $l_{m}(t)$ is muscle length at time $t$ and appropriate value for $\gamma$ should be determined by experiments for different muscles.

The second method is to directly adjust the force-length curve, which shifts the curve directly by adjusting key track points to different activations. As AFLC models use the Bezier spline to represent the relationship, there are generally three controllable points representing the active slack length, transition length and maximum length of force exertion respectively and one fixed point representing the optimal fiber length. In fact, the three controllable points are constant during one simulation. Two points representing the optimal fiber length and active slack length are adjusted to activations using Eqs (4) and (6), respectively.

$\displaystyle l_{a}(t)=l_{a}(\gamma(1-a(t))^{2}+1)$ (6)

where $\gamma$ is the percentage change in optimal fiber length, $l_{a}$ is active slack length at maximum activation, $a(t)$ and $l_{a}(t)$ are activation and active slack length at time $t$ , respectively. Then the curve is redrawn with the varying activation at each time step. And the results of the two equations shift the active slack length less than optimal fiber length and ensure the fixation of maximal length of force exertion.

Figure 4.

Comparison of experimental and simulated force profiles estimated by the first AFLC model ( $\gamma=$ 0.30) in the two sets of submaximal benchmarking experiments. Results in the left columns were for maximum amplitudes of 1 mm and those in the right ones were for 8 mm. (A) Constant-frequency stimulation rates of 10, 20, 30 Hz were applied. (B) Random signals with mean stimulation rates of 10, 20, 30 Hz were applied. Note that experimental data are from Perreault et al. [2].

3. Biological benchmark simulations, results and discussion

3.1 Biological benchmark simulations

Millard et al. [7] provided a series of submaximal-activation biological benchmark data obtained from Perreault et al. [2], which contained forces measured from submaximally-activated cat soleus muscle in vivo. There were three groups of experimental trials and the first group included six trials, in which the musculotendon actuator was held at a constant length of $l_{s}^{T}+l_{o}^{M}\cos\alpha_{o}-4$ mm and stimulated at the mean rates of 10, 20 and 30 Hz constantly and randomly. The activations calculated by inverting the common Hill-type muscle model are shown in Fig. 3. The other two groups both had six trials and shared the same mode of stimulations with the first one, but applied length changes with maximum displacements of 1 mm and 8 mm to the free end of the tendon, respectively. The length changes were generated by low-pass filtering a normally distributed random waveform and then scaling it to corresponding amplitudes. Data of the first isometric group was used to calculate the activation signals elicited by three stimulation frequencies and these activations could be applied to musculotendon models in simulations of the following two non-isometric groups. Experimental protocols and muscle parameters were detailed by Perreault et al. [2] and Millard et al. [7]. All benchmark simulations were performed on a personal computer with a 3.6 GHz eight-core Intel processor.

It’s difficult to define activations in Hill-type models and to avoid effects of activation processes on model errors, they were achieved by calculating isometric experimental data inversely. In the process of inverse calculation, musculotendon length was constant and fiber velocity was zero, which greatly simplified the model and achieved the corresponding muscle activations during different kinds of simulations. As AFLC models incorporate muscle activation and force-length coupling, it was complicated to calculate activations from the measured forces, so activations calculated from the Millard damped model, which we took as the common Hill-type model, were adopted. Then the same activations and time-varying musculotendon lengths were input to the Millard model and two AFLC models. Errors between the experimentally measured muscle forces and those predicted by the three kinds of models were quantified using percent root mean square errors (RMSE), i.e.,

$\displaystyle\textit{RMSE}=\sqrt{\sum\nolimits_{N}(F_{\exp}-F_{\text{est}})^{2% }/\sum\nolimits_{N}F_{\exp}^{2}}\ast 100\%$ (7)

where $F_{\exp}$ is the experimentally-measured force and $F_{\text{est}}$ is the model-estimated force.

Figure 5.

Comparison of experimental and simulated force profiles estimated by the second AFLC model ( $\gamma=$ 0.30) in the two sets of submaximal benchmarking experiments. Results in the left columns were for maximum amplitudes of 1 mm and those in the right ones were for 8 mm. (A) Constant-frequency stimulation rates of 10, 20, 30 Hz were applied. (B) Random signals with mean stimulation rates of 10, 20, 30 Hz were applied. Note that experimental data are from Perreault et al. [2].

3.2 Results

Comparisons of experimentally-measured and two AFLC model-estimated force profiles from the submaximally-activated muscle undergoing length changes of maximum displacements of 1 mm and 8 mm are shown in Figs 4 and 5, separately. Throughout all the trials, model errors increased with increasing length change amplitudes and the relation between two curves in each trial had the trend that estimated forces were smaller at first and then kept larger than the measured ones. In constant-frequency stimulation trials, waveforms of the force change were similar in each column except for the force range size and the similarity was especially obvious between groups of 20 and 30 Hz. Moreover, there was no significant difference for curves between constant and random stimulation trials when the frequency reached 30 Hz.

The RMSEs of each model were compared in Fig. 6. Generally, two AFLC models performed better than the Millard model and exceeded the Millard model in all trials except the 10 and 20 Hz random trials of 1 mm displacement. To summarize, errors of the Millard model were less than 3.35% and 16.59% corresponding to displacement of 1 mm (top row in Fig. 6) and 8 mm (bottom row), respectively. Meanwhile, those of the first AFLC model were less than 2.47% and 13.98%, and the second 2.82% and 13.81% correspondingly.

However, the second AFLC model, which was a curve-modified version, cost much more time than the first one, as shown in Fig. 7. And the average running time of the second model is nearly 17 times that of the first one in each of the corresponding trial.

Figure 6.

Comparison of the root mean square errors (RMSE) of three musculotendon models for submaximally-activated muscle undergoing length changes of maximum displacement 1 mm (top) and 8 mm (bottom).

Figure 7.

Comparison of the average running time of each trial using two AFLC models. (A) Maximum amplitude is 1 mm. (B) Maximum displacement is 8 mm.

Figure 8.

Changing average RMS errors at different values of $\gamma$ factor ranging from 0 to 0.6.

In addition, value of $\gamma$ factor can also affect the simulation results, the optimal value for the first AFLC model in this simulation was about 0.32, which was achieved by calculating the average RMSE of each trial at different $\gamma$ values, as shown in Fig. 8. However, the default value adopted by Lloyd and Besier was 0.15, which was calibrated using the musculoskeletal model in their pilot work.

3.3 Discussion

This paper constructed two practicable AFLC musculotendon models, evaluated their biological accuracy and computational speed using submaximal benchmark data provided by Millard et al. [7], and further validated the activation-force-length coupling from the perspective of applications. The AFLC models incorporated the advantages of Millard damped equilibrium model, which solved the numerical singularities effectively, and the coupling between activation and force-length property, which affected the applicability of common Hill-type models. The model parameter $\gamma$ we’ve calculated out for cat soleus muscle is 0.32 and the first AFLC model is recommended in the construction of upper-layer musculoskeletal models.

Throughout all the trials, model errors increased with increasing length change amplitude, but not linearly. This makes sense as larger displacements brought a wider range of length changes in the model and greater challenges to the model robustness. As the musculotendon models are second-order models, the model errors changed nonlinearly with lengths. Meanwhile, muscle velocities also increased as length change amplitude increased, which might also be one cause of larger model errors. And there was a similar pattern in the change of relationship between experimental and simulation curves, in which simulation results underestimated experimental data at the early stage and overestimated that thereafter. This phenomenon indicated that there might be some time-variant factor that affected the model accuracy. For example, muscle fatigue will reduce muscle force-generating capabilities with time and is also quite common in daily activities [16].

As the stimulation frequency increased, muscle forces increased because more muscle fibers were activated and when all the fibers were activated, higher frequency could change muscle contraction mode from single twitch to tetanic contraction and generate more contraction forces. Therefore total muscle force magnitude did not increase linearly with the stimulation frequency. In constant stimulation trials, the magnitude of activations calculated at 20 Hz was closer to that at 30 Hz and accordingly their muscle force profiles were also similar. And the similarity between constant and random stimulation trials at the frequency of 30 Hz indicated that muscles almost came into the stage of tetanic contraction at this frequency. However, larger errors also came up with higher stimulations, it might also because of a wider range of forces.

In the 10 and 20 Hz random trials of 1 mm displacement, AFLC models increased the error a little, one reason might be the inaccurate activation. In fact, the activation, one of the inputs to the simulations, was inversely calculated using the Millard damped model, results should be better if it could be calculated using AFLC models Another reason could be the random conditions, which were not controllable. This also reflected the lacking of experimental data in some way and more experimental data could also be instructive in the selection of $\gamma$ value, which represented the coupling relationship between activation and force-length property. It should be noted that the value of $\gamma$ as fixed to 0.15 was an optimal choice in the EMG-driven musculoskeletal model of Lloyd and Besier [8], but our result was 0.32. In view of the differences between our experiments, we hypothesized that $\gamma$ value might vary with different muscles and application levels.

The first model was easier to implement and moved the whole force-length curve to right indirectly, but ignored the change of distance between the active slack length and the optimal fiber length, also neglected the property of fixed maximal fiber length. The second one took these properties into account but obtained no obvious improvement at the cost of much more time for redrawing the curve with new parameters at each time step. Therefore, the first modified model is recommended to adopt in the construction of upper-layer musculoskeletal models.

Hill-type musculotendon models are easy to implement and represent the major characteristics of muscle contraction activities. They are widely integrated into musculoskeletal models to simulate human movements. But they will cause inaccuracy problems if applied to the simulation of human daily activities. Except for activation force-length coupling, which we tried to combine into current models in this paper, there are still other issues that need further considerations. One study suggested scaling the maximum shortening velocity with activation when using Hill-type models in musculoskeletal modelling as it depended on activations [17]. And another study testified an obvious improvement in accuracy and confidence of their forward dynamic simulation results by varying the maximum shortening velocity with activations [18]. It seemed that there were more couplings than we expected, which highlighted in submaximal situations. Meanwhile, some time-related factors also need to be considered in real-life scenarios. For example, muscle contraction history could have an effect on the isometric force, to be specific, at the same length, muscles would generate lower isometric force after active shortening and higher after active stretching [19]. Another example, muscle fatigue that we mentioned before when explaining the overestimation, is also a common time-dependent physiological phenomenon in daily activities. However, the underlying mechanism of muscle fatigue is complicated and not well understood yet. Many studies have taken it as a research topic, which will promote the construction of better applicable muscle models [16, 20]. Furthermore, experimental studies have found that muscle contraction dynamics also depended on the muscle fiber-tendon length ratio, which varies with different muscles [21]. Gareis et al. [22] studied muscles in the hind legs of the cat and built the force-length relation for nine different skeletal muscles experimentally. It turned out that the patterns of force-length relations varied widely among different muscles. This could be dependent on the muscle architecture, such as ratio of fiber to tendon length, pennation pattern and angle, and the muscle function determined by the locations of attaching points and geometry paths [23]. To sum up, there’s still a lot of work to do when integrating Hill-type musculotendon models into musculoskletal models. In the meantime, more benchmark experiments and data are needed to evaluate these models.

In consideration of the limits of Hill-type musculotendon models and what we have achieved, future work can be divided in the following steps. Firstly, validating our AFLC models using natural activations, experiments like crossed-extension reflex (CXR), by which normal patterns of motor unit recruitment and rate modulation can be obtained. Secondly, as muscle properties vary among different muscles, more tests should be performed on muscles other than cat soleus muscles. Thirdly, integrating our first model into human musculoskeletal models and simulate daily movement tasks for further validations. Finally, we could consider other coupling relationships [18, 24] and submaximal-related factors [23] that should be incorporated into AFLC models to create more accurate submaximal models. In conclusion, more efforts should be devoted to the study of Hill-type musculotendon models for further applications.

Footnotes

Acknowledgments

The authors wish to acknowledge the support of the National Natural Science Foundation of China (61431017) and the assistance of Thomas Uchida in implementing the benchmark simulations.

Conflict of interest

None to report.

References

Erdemir

McLean

Herzog

van den Bogert

. Model-based estimation of muscle forces exerted during movements. Clin Biomech 2007; 22: 131-54. doi: 10.1016/j.clinbiomech.2006.09.005.

Perreault

Heckman

Sandercock

. Hill muscle model errors during movement are greatest within the physiologically relevant range of motor unit firing rates. J Biomech 2003; 36: 211-8. doi: 10.1016/S0021-9290(02)00332-9.

Hill

. The heat of shortening and the dynamic constants of muscle. Proc R Soc B Biol Sci 1938; 126: 136-95. doi: 10.1098/rspb.1938.0050.

Zahalak

. Modeling muscle mechanics (and energetics). in: Winters

Woo

S-Y

, editors. Mult Muscle Syst SE-1, Springer New York; 1990, p. 1-23. doi: 10.1007/978-1-4613-9030-5_1.

Winter

. Biomechanics and motor control of human movement. D.A. Winter. Biomech Mot Control Hum Mov 2009; 74: 94.

Zajac

. Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor control. Crit Rev Biomed Eng 1988; 17: 359-411.

Millard

Uchida

Seth

Delp

. Flexing computational muscle: Modeling and simulation of musculotendon dynamics. J Biomech Eng 2013; 135: 21005. doi: 10.1115/1.4023390.

Lloyd

Besier

. An EMG-driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo . J Biomech 2003; 36: 765-76.

Mau-Moeller

Jacksteit

Jackszis

Feldhege

Weippert

Mittelmeier

, et al. Neuromuscular function of the quadriceps muscle during isometric maximal, submaximal and submaximal fatiguing voluntary contractions in knee osteoarthrosis patients. PLoS One 2017; 12: e0176976. doi: 10.1371/journal.pone.0176976.

10.

Andersson

Kvernes

. Bezier and B-spline Technology. Master Sci Thesis 2003; 58.

11.

Wanner

Hairer

. Solving ordinary differential equations II. Stiff Differ Probl 1991. doi: 10.1007/978-3-642-05221-7.

12.

Huijing

. Muscle, the motor of movement: Properties in function, experiment and modelling. J Electromyogr Kinesiol 1998; 8: 61-77. doi: 10.1016/S1050-6411(97)00023-0.

13.

Stephenson

Wendt

. Length dependence of changes in sarcoplasmic calcium concentration and myofibrillar calcium sensitivity in striated muscle fibres. J Muscle Res Cell Motil 1984; 5: 243-72. doi: 10.1007/BF00713107.

14.

Hatze

. A myocybernetic control model of skeletal muscle. Biol Cybern 1977; 25: 103-19. doi: 10.1007/BF00337268.

15.

van Zandwijk

Bobbert

Baan

Huijing

. From twitch to tetanus: Performance of excitation dynamics optimized for a twitch in predicting tetanic muscle forces. Biol Cybern 1996; 75: 409-17.

16.

Tang

Stojanovic

Tsui

Kojic

. Modeling of muscle fatigue using Hill’s model. Biomed Mater Eng 2005; 15: 341-8.

17.

Chow

Darling

. The maximum shortening velocity of muscle should be scaled with activation. J Appl Physiol 1999; 86: 1025-31.

18.

Camilleri

Hull

. Are the maximum shortening velocity and the shape parameter in a Hill-type model of whole muscle related to activation? J Biomech 2005; 38: 2172-80. doi: 10.1016/j.jbiomech.2004.09.037.

19.

McGowan

Neptune

Herzog

. A phenomenological muscle model to assess history dependent effects in human movement. J Biomech 2013; 46: 151-7. doi: 10.1016/j.jbiomech.2012.10.034.

20.

Green

Stefanovic

Warman

Askew

. Empirical modelling of the dynamic response of fatigue during intermittent submaximal contractions of human forearm and calf muscles. J Electromyogr Kinesiol 2015; 25: 20-7.

21.

Mörl

Siebert

Häufle

. Contraction dynamics and function of the muscle-tendon complex depend on the muscle fibre-tendon length ratio: A simulation study. Biomech Model Mechanobiol 2015. doi: 10.1007/s10237-015-0688-7.

22.

Gareis

Moshe

Baratta

Best

D’Ambrosia

. The isometric length-force models of nine different skeletal muscles. J Biomech 1992; 25: 903-16. doi: https://doi.org/10.1016/0021-9290(92)90230-X.

23.

Huijing

. Parameter interdependence and success of skeletal muscle modelling. Hum Mov Sci 1995; 14: 443-86.

24.

Biewener

Wakeling

Lee

Arnold

. Validation of hill-type muscle models in relation to neuromuscular recruitment and force-velocity properties: predicting patterns of in vivo muscle force. Integr Comp Biol 2014; 54: 1072-83. doi: 10.1093/icb/icu070.