Journal of Biomechanics 47 (2014) 3448–3453
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
A bioenergetic model for simulating athletic performance of intermediate duration Gilbert Gede n, Mont Hubbard University of California Davis, Mechanical and Aerospace Engineering, 1 Shields Ave, Davis, CA 95616, United States
art ic l e i nf o
a b s t r a c t
Article history: Accepted 17 September 2014
Simulating factors affecting human athletic performance, including fatigue, requires a dynamic model of the bioenergetic capabilities of the athlete. To address general cases, the model needs inputs, outputs, and states with a set of differential equations describing how the inputs affect the states and outputs as functions of time. We improve an existing phenomenological muscle model, removing unnecessarily fast dynamic behavior, adding force–velocity dependence, and generalizing it to task level activities. This makes it more suitable for simulating and calculating optimal strategies of athletic events of medium duration (longer than a sprint but shorter than a marathon). To examine the validity and limitations of the model, parameters have been identified from numerical fits to published experimental data. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Bioenergetics Model Fatigue Simulation Muscle Model Athletic Performance
1. Introduction When considering athletic events that involve racing against the clock such as a 5 km road race or cycling time trial, the following question inevitably arises: what is the optimal pacing strategy? Commonly practiced strategies range from maintaining a constant power output to finishing with a sprint. The question is especially difficult to answer in situations where the athlete's speed and power output vary over the duration of the event (e.g., a course with hills or wind). The first step in answering it consists of determining a model to represent the athlete's energetic dynamics in the event. Numerous authors have examined competitive athletic pacing. A runner's optimal strategy for different distances was calculated and the results were compared to the then current world records (Keller, 1974). This model used the concept of an oxygen debt. Anaerobic metabolism effects on energy production have also been investigated (Ward-Smith, 1999), but in a manner which is not applicable to longer events or events with a variable power output. The minimum time problem for running has been studied using the critical power model by Morton (2009) who found that the optimal solution is all-out for the entire race. All these efforts have limitations, either in fidelity or extensibility to more complex events. Simulating an athlete's performance in a general situation requires a dynamic model with a control input which governs
n
Corresponding author. E-mail addresses: [email protected] (G. Gede), [email protected] (M. Hubbard). http://dx.doi.org/10.1016/j.jbiomech.2014.09.017 0021-9290/& 2014 Elsevier Ltd. All rights reserved.
the athlete's output. While such dynamics models have existed for decades, a more recent approach (Liu et al., 2002) divides the muscle into three compartments: resting, active, and fatigued. This model has also been extended (Xia and Frey Law, 2008) to include multiple muscle fiber types and the control logic required to activate the fibers in a realistic fashion. More recently, a 4 compartment, single fiber type model was proposed which appears to faithfully replicate experimental behavior (Sih et al., 2012). However, if the goal is to perform optimal control calculations which can tell an athlete what power output, as a function of time, they should produce to achieve the minimum elapsed time in a cycling time trial, none of these models are appropriate. The models above all have significant weaknesses which prevent them from filling this role. To best answer the question, “what is the optimal strategy?”, we present a new phenomenological model which represents as much behavior as needed and no more. In this paper, we present a new dynamic model, designed for answering pacing strategy questions. It is an improvement of the existing phenomenological model of Xia and Frey Law (2008); it is more appropriate for general simulation and use in optimal control calculations. It does not attempt to replace current knowledge on muscle physiology, nervous system control, or decades of clinical experiments. Instead, it attempts to represent what an athlete can do, without considering the details of the physiological mechanisms of how this happens. Compared to current models, it removes unnecessary detail and adds important features. A limited validation using published experimental data is shown. Finally, interpretations of the identified physiological parameters are discussed.
G. Gede, M. Hubbard / Journal of Biomechanics 47 (2014) 3448–3453
3449
2. Methods A qualitative overview of the model is first presented, followed by a mathematical description, and the methods used for its validation.
2.1. Background Compartment theory has been used for modeling a single muscle (Liu et al., 2002; Xia and Frey Law, 2008). Each motor unit (MU) pool within the muscle is a collection of a finite number of similar fibers, but for modeling purposes it is convenient to treat it instead as continuous. The continuous pool is apportioned into three compartments: resting fraction Mr, active fraction Ma, and fatigued fraction Mf. The dynamics of the model describe the transfer between compartments: active muscle becoming fatigued, fatigued muscle recovering, and resting muscle becoming active. Fig. 1 shows a schematic of this process during and after a maximal effort exercise. We have taken the approach of Xia and Frey Law (2008) a step further by applying it to groups of muscles, and use a 3 compartment model with multiple MU pools to represent task-level activities such as cycling, as the model of Sih et al. (2012) did with success. The fundamentals of the 3 compartment model are unchanged: the amount of force produced is proportional to the amount of active muscle, the rate at which fatigue accrues is proportional to how much muscle is active, and the rate of recovery is proportional to the current level of fatigue. There are differences in the model we present however, in the rate of activation, elimination of redundant quantities, and the inclusion of the force–velocity curve, among other changes. Additionally, the original model was developed for examining the behavior of motor unit pools within a single muscle. The present model was formulated with the assumption it could also represent a group of muscles. Fig. 2 shows how the compartments are allocated within a single motor unit pool.
Fig. 1. A schematic representation of muscle transfer between compartments before, during, and after a maximal effort exercise. Initially the muscle is completely at rest. Once maximum effort starts, no muscle is at rest, and fatigue builds up. At the end of the effort period, recovery begins with active muscle being released and fatigued muscle recovering.
2.2. Model dynamics Unlike the Xia and Frey Law (2008) model, we assume there are only slow and fast fibers (but no intermediate ones). The quantities involved with the slow oxidative fiber MU pool are identified by a superscript o and the fast glycolytic fiber MU pool quantities by a superscript g. The differential equations for the single state variable for each pool, the fatigued muscle fraction Mf, are g
dM f dt
¼ M ga F g M gf Rg
ð1aÞ
o
dM f ¼ M oa F o M of Ro dt g
g
o
ð1bÞ o
where F , R and F , R are the fatigue and recovery rates for the fast and slow fiber types, respectively. As will be shown, only one differential equation per pool is required. This has been accomplished by removing mathematically redundant equations and instead partitioning the motor unit pools into fatigued and nonfatigued compartments (where it is assumed the non-fatigued compartment can be instantaneously apportioned between resting and active). This model, which has 3 compartments per pool, can be represented with only one state variable for each MU pool for two reasons. First, there is clearly a redundant state, due to the total muscle size constraint 1 ¼ Ma þ Mf þ Mr
ð2Þ
making one state dependent on the other two. Second, by not considering the rate of muscle activation, Ma no longer needs to be considered a state variable. The purpose of the present model is to simulate athletic events of medium duration (longer than a sprint but shorter than a marathon). We claim that the sub-second (much shorter than a second) time dynamics of muscle fiber recruitment are not relevant to athletic events of this and longer durations. If these sub-second dynamics are not considered, the nonfatigued muscle in a MU pool can be instantaneously split between the active and resting compartments. This allows for the control input to the system to simply be the amount of active muscle Ma, rather than the recruitment rate of active muscle. These new control inputs (Ma for each fiber type) have constraints. Within each pool the amount of active muscle must be non-negative, and it cannot be greater than the amount of non-fatigued (i.e. available for recruitment) muscle 0 r M ga r 1 M gf
ð3aÞ
0 r M oa r 1 M of
ð3bÞ
There is also a constraint imposed by the prioritization of fiber recruitment. Previous compartment models have used a logic sequence to enforce the correct order of fiber recruitment (Henneman et al., 1965). The following constraint
Fig. 2. A visualization of the compartments within a MU pool. The bar dividing the fatigued and non-fatigued compartments moves down at a rate which is the difference of the recovery ðM f RÞ and fatigue ðM a FÞ. The system control Ma is bounded between zero and the amount of non-fatigued muscle. The resting compartment Mr is simply the non-active, non-fatigued fraction.
equation enforces the same behavior: M ga M oa þ M of 1 ¼ 0
ð4Þ
The allowable space for the control inputs to the system is shown as a bold line in Fig. 3. This constraint enforces that the amount of active fast fiber can be nonzero only when all of the slow fibers are active. Using a constraint equation is necessary when applying optimal control techniques, because logic sequences introduce nonsmoothness, possibly preventing an optimal solution from being found. Outside optimal control applications, an equivalent logic sequence is acceptable. The final change from other compartment models is the incorporation of the force–velocity curve within the equations. For medium duration events an athlete is unlikely to spend time at either extreme of the force–velocity curve, and we claim that an exponential force–velocity approximation is sufficient. This allows the
3450
G. Gede, M. Hubbard / Journal of Biomechanics 47 (2014) 3448–3453
τ ¼ 1=ðF þ RÞ, is the time constant of the exponential decay to the steady-state force output, when the muscle is at maximum effort (see Fig. 1). The fiber type with the smaller time constant (shorter period) is the faster fiber, as it will take less time to become fatigued. We now have an expression (7) for the torque output of each fiber type as a function of time and crank speed assuming that the system parameters are known. A least squares algorithm will be used to compute parameters to fit the experimental data, scipy:optimize:leastsq (Jones et al., 2001). The eight model parameters (S0, c, F, and R for the two fiber types) are found by comparing the predicted torque based on the exact solutions of the differential equations to the experimental torque. The error is simply the difference between the measured and the predicted torque at each time. Because our model does not consider muscle activation time, we consider the experimental data only after the peak torque has been reached. 2.4. Experiments
Fig. 3. Allowable configuration space for the two inputs, Moa and Mga, is shown by the bold line. This configuration space enforces the prioritized recruitment behavior required by the Henneman et al. (1965) size principle. final muscle model output force to be expressed as F ¼ Sg0 M ga e c v þ So0 M oa e c g
ov
ð5Þ
where Sg0, cg and So0, co are the muscle size and force–velocity curve coefficients, for the fast g and slow o fibers respectively, v is the muscle speed (or an analog thereof, such as crank speed when cycling or ground speed when running, which is always positive), an exogenous input, and F is the muscle force output. 2.3. Validation methods The validation presented below consists of selecting appropriate experimental data, solving equations, and numerical parameter fitting. Because it is difficult to measure the individual activation amount of both slow and fast motor unit pools, we have limited the types of experiments used for validation. The selected experiments are maximal effort tests, in which the amount of active muscle will always be the maximum possible for both fiber types: M a ¼ 1 M f , i.e. Mr ¼0. By selecting maximal effort tests with recorded force or power outputs and recorded velocities (or analogs thereof), the parameter fitting process will identify Fo, Fg, Ro, Rg, co, cg, So0, and Sg0 in (1) and (5). However, the model in general is not limited to the maximal effort case. The experiments are task level (cycling) activities, involving collections of joint groups which are themselves collections of muscles. Thus far, the model has been intentionally described without a scale of application. For the task level activity of cycling, we attempt to show that this level of model abstraction can accurately represent the combined behavior of all the muscles involved. At this point, two substitutions in the equations will be introduced, in order to use them with cycling. The muscle force output F is instead substituted with crank torque T, and muscle velocity v with crank angular velocity (speed) ω. This represents a transformation from linear to rotational variables; aside from the change in symbols and units Eq. (5) is unaffected. Because the differential equations presented here are relatively simple, they can be solved exactly both during maximal effort tests and during rest. After solving the differential equations for Mof and Mgf , the amount of active muscle is found from M a ¼ 1 M f . The total torque output is the sum of the torque output of each fiber type. During a maximal effort test the solution to the differential equations for a single MU pool (one fiber type) is F F þ M f ðt 0 Þ M f ðtÞ ¼ ð6Þ e ðt t0 ÞðF þ RÞ F þR F þR
Two experiments suitable for validation have been published, both using a cycling ergometer where power was a quadratic function of crank speed: P ¼ r ω2 . The power dissipated P is the product of a resistance r and the speed ω squared. Unlike Coulomb friction-braked ergometers, this type of ergometer is associated with a gradual decrease in crank torque produced by the athlete over a maximum effort test. In contrast, friction-braked ergometers show an initial peak in torque (due to flywheel acceleration), followed by a constant torque (maintaining equilibrium against friction), then an abrupt end (inability to maintain the torque required to balance friction). One experiment (variable tempo) consisted of maximum effort tests at three resistive loads, with fully rested subjects (Vanhatalo et al., 2008). Power and speed were recorded while subjects performed with low, medium, and high resistances, giving high, medium, and low peak speeds and end-test speeds (see Fig. 4). A similar experiment (rest interval) used different subjects (see Fig. 6). It involved periods of maximum effort expenditure separated by rest intervals varying in duration (Vanhatalo and Jones, 2009). In each test fully rested subjects exercised at maximal effort for 30 s, followed by a rest period of either 2 or 15 min, followed by another maximal effort exercise for 3 min. For a fully rested athlete, it was assumed that the 30 s expenditure was identical during each test, and was identical to the first 30 s of a 3 min control test. Data from the figures in each publication were digitized using software (Tummers, 2006). In both experiments the figures displayed results averaged across many subjects. Individual subject data would be preferable for these validation analyses, but we believe there are still worthwhile insights to be gained from examining the averaged data.
3. Results For both experiments crank torque was computed by dividing the reported power by the reported crank speed. A least squares fit was used to compute a single set of eight model parameters (four for each MU pool) in each experiment. To show the ability of the model to work over a range of muscle speeds and recovery periods, parameters were found using the data from all tests. For
The output torque for each MU pool is then TðtÞ ¼ S0 e cωðtÞ ð1 M f ðtÞÞ
ð7Þ
When starting from a fully rested state there is no fatigued muscle and M f ðt 0 Þ ¼ 0. To find Mf at time t during a rest period, after exertion ending at t1, the differential equations are again solved to yield M f ðtÞ ¼ M f ðt 1 Þe Rðt t 1 Þ
ð8Þ
We also show two quantities which are calculated from the model parameters: the steady-state force fraction and the time constant. The steady-state force output for each fiber type is defined as G ¼ 1 ðF=ðF þ RÞÞ. This is the force that can be still produced by each fiber type when the athlete has reached the maximum sustainable level of fatigue. The time constant for each fiber type, defined as
Fig. 4. Crank speed from three tests with different ergometer resistances, leading to different tempos (Vanhatalo et al., 2008).
G. Gede, M. Hubbard / Journal of Biomechanics 47 (2014) 3448–3453
both experiments, the identified parameters which were either very large or very small in magnitude were rounded to infinity or 0, respectively. The difference between the exact output of the least squares solver and the rounded parameters was less than machine precision, justifying the use of the rounded parameters. The variable tempo experiment included power and speed from maximal effort cycling ergometer tests (Vanhatalo et al., 2008). Fig. 4 displays the speed ω in the tests, the exogenous input term used in parameter identification. Fig. 5 displays torque for the three tests, found by dividing the reported power by the speed. Also plotted are the fit lines showing the predicted model torque resulting from parameter identification. The RMS torque prediction error is 0.49 Nm. The variable tempo parameters are in Table 1. The rest interval experiment also measured power and speed from maximal effort cycling ergometer tests (Vanhatalo and Jones, 2009). Fig. 6a and b shows the calculated speed, derived from the reported power and maximum speed, using the equation defining the ergometer's power/speed relationship. Fig. 6c and d portrays the torque, calculated by dividing the recorded power by calculated speed. Also plotted in Fig. 6c and d are the fit lines for predicted torque from the identified parameters. Fig. 6a and c shows speed and torque from the 30 s effort; Fig. 6b and d plots speed and torque from the 3 min efforts. The RMS torque prediction error is 0.37 Nm. The rest interval parameters are also in Table 1.
4. Discussion In formulating the model, we have assumed that equations which govern the behavior of MU pools within a single muscle
Fig. 5. Experimentally calculated (data from Vanhatalo et al., 2008) and model fit torque for three maximal effort tests with different ergometer resistances.
Table 1 Parameters generated by least squares fit of experimental data for both the variable tempo experiment and the rest interval experiment. Parameter
S0 (N m) F ðs 1 Þ R ðs 1 Þ c ðs 1 Þ G τ (s)
Variable tempo test
Rest interval test
Slow
Fast
Slow
Fast
9.58 0.0062 0 0.0186 0 162
371 0.042 0.014 0.146 0.25 17.8
10.2 1.7E þ 5 Inf 0 1 0
42.3 0.017 0.012 0 0.407 34.8
3451
could also be used to describe the behavior of a muscle group or task level activity involving many muscle groups (which each have their own MU pools). Even though the original model (Xia and Frey Law, 2008) referred to individual muscles, we cannot claim our model represents individual muscles, because it has not been experimentally validated on this scale. Instead, we only examine the model's ability to replicate the behavior of task level activities. The success of this model in fitting the examined experimental behavior shows that it is apparently able to represent such an activity (cycling). For each experiment, the distinction between the slow and fast MU pools is made by examining the time constants. The pool which reaches full fatigue sooner is the fast pool, generally identified by a smaller time constant τ. Examining the parameters found from the variable tempo experiment (Table 1), we see that a relatively large isometric muscle size S0 has been found for the fast fibers, compared to that of the slower fiber MU pool. The faster muscle pool has a time constant of 17.8 s and a static torque capability of 371 Nm. The slower pool's time constant is much longer, 162 s, with a much smaller static torque capacity of 9.58 Nm. Admittedly, the variation in MU pool size does not seem to be reasonable: the value for S0 for the fast fiber MU pool is approximately 40 times that of the slower MU pool. While these static torque values appear to be significantly different than what might be expected physiologically, they must be viewed in the context of the force–velocity curve approximation. For the average crank speed in the standard target speed run, the actual torque, Se cω , is lower for each MU pool: 82.2 Nm for the fast pool, and 7.91 Nm for the slow pool. This is shown more clearly in Fig. 7, where the force velocity curve for each fiber type is plotted against the speed range for the experiments. Also indicated are the average crank speed and the time distribution of the crank speed. This indicates that the operational difference between S0 for the slow and fast fiber groups is not as large as the parameters would suggest. Discrepancies such as this are to be expected, because the model contains an approximation to the force velocity curve. Interestingly, when fitting parameters to a single one of the three variable tempo tests, the identified force–velocity curve parameter c and recovery coefficients R are very small. This shows that these model components, recovery and the force–velocity curve, come into play only when used with test data which exhibit recovery behavior or significantly different speed ranges. The rest interval experiment was conducted with different subjects from those in the variable tempo ones, and thus identical parameters should not be expected. In this experiment very small force–velocity coefficients were identified for both MU pools. Additionally, the slow MU pool showed very large fatigue and recovery rate coefficients. It is important to recognize though that these coefficient only appear as parts of terms in (6). Therefore the magnitude of the coefficients is less important than how they interact within that equation. The previously defined terms G and τ (steady-state force output and time constant) provide an interpretation of these ratios. In this experiment the fast MU pool has a time constant of 34.8 s, the slow MU pool has a time constant of essentially 0 s. At first glance, this seems to be inconsistent with the previous definition of fast and slow fibers. However, there is a steady-state force output of 100% for the slow pool. This indicates that the slow pool immediately reaches its steady-state value, which is equal to its unfatigued value. Its maximum force output does not change. Between the two experiments (variable tempo and rest interval) there is a similar size ratio between the slow and fast fibers. Over the operational speed ranges, the rest interval experiment shows a factor of 4 between fast and slow pool torque outputs, while the variable tempo experiment shows a factor of 10 between fast and slow pool torque outputs. The fast MU pool fatigue rate
3452
G. Gede, M. Hubbard / Journal of Biomechanics 47 (2014) 3448–3453
Fig. 6. Experimentally calculated speeds (a and b) and experimentally calculated and numerically fit torques (c and d) for 30 s and subsequent 3 min maximal effort expenditures, after different rest intervals (experimental data from Vanhatalo and Jones, 2009). Note that the 30 s data is assumed to be identical in both cases, as the athletes were completely rested before each test.
coefficients identified from the experiments are of the same order of magnitude: 0.042 vs. 0.017. Also, in both experiments the slow fibers are found to have little or no decrease in force production over the test duration; only the fast fibers fatigue appreciably. Considering all these points, there is some reassurance in the model's self-consistency across various test conditions. 4.1. Examination of reduced models
Fig. 7. Fast and slow MU pool force–velocity curves identified from the variable tempo experiment. Also indicated are the average speed (dashed line) and the distribution of speeds (shaded grey) over the experiment. There is a large difference between the fast and slow static torque parameters S0, 370 N m and 10 N m. As the static torque parameter appears in S0 M a e cω , their output is reduced by the force–velocity curve. In this case, over the main operating range (9–11 rad/s) the difference between product of the static torque parameter and the force velocity curve approximation of the slow and pools is 80 N m and 8 N m, a much smaller difference than the base values of S0 suggest.
Examining the results in Table 1 should raise concerns over the accuracy of a human energetics model which involves coefficients of 0 or infinity. In addition to the preceding discussion of the fit parameters, we have examined how the results of the fitting process are influenced by examining the maximal effort case; specifically how some model components combine, which leads to these unrealistic parameters. By examining the parameters when considering only the fullyrested case (i.e. M f ðt 0 Þ ¼ 0), the analytic solutions to the differential equations are greatly simplified. When there is no initial fatigue, the torque output equation can be rewritten in terms of G and τ; for a single MU pool, the torque output is ð9Þ TðtÞ ¼ S0 e cωðtÞ G þ ð1 GÞe t=τ Re-examining the parameters in Table 1, we can see how the identification of slow and fast fibers can be made. For the variable
G. Gede, M. Hubbard / Journal of Biomechanics 47 (2014) 3448–3453
tempo experiments, the slow fibers' torque output is TðtÞ ¼ 9:58e0:0186ωðtÞ e t=162 and the fast fibers' torque output is TðtÞ ¼ 371e0:146ωðtÞ 0:25 0:75e t=17:8 . Examining these equations shows that the fast fibers will clearly reach their steady-state output sooner. However, if these model parameters were used in a longer simulation, the model would become inaccurate, as it would suggest that all the slow fibers would be fatigued and only fast fibers would generate torque – something that is not consistent with physiology. One clear conclusion to draw from this is to generate the parameters from longer test sessions. For the rest interval experiments, we obtain the following equations: the slow MU pool torque is TðtÞ ¼ 10:2 and the fast MU pool torque is TðtÞ ¼ 42:3 0:407 þ 0:593e t=34:8 . The slow MU pool is identified as the pool which shows increased resistance to fatigue. In this case it shows no change at all! The fast MU pool does show signs of fatigue however. This raises the question: is there a need for the complete model that has been presented in this paper? The answer is unequivocally yes. It is important to remember that these submodels are the result of the numerical parameters being substituted into the solutions to the differential equations solved for the maximal effort case. There are no simpler forms of the differential equations which can be solved and manipulated into all the forms above. This reduced model analysis does reveal some model weaknesses. In many dynamic systems when the inputs are saturated the solutions to the differential equations become simplified. This is what has happened here. This also shows that the model needs more thorough parameter identification tests; the parameters shown for the maximal effort case are not completely translatable to the general case (i.e. non-maximal effort). Parameter identification studies are more effective when the system inputs are known. Because we have not considered the mechanisms for MU pool activation, the inputs in this case cannot be known exactly. This might be a potential weakness of the present model when attempting to find parameters in more general applications. 5. Conclusions and future work We have presented a model which reproduces the experimental dynamic energetic behavior of human subjects while cycling, displaying both fatigue and recovery behavior. Although further validation of the model is required, we believe it may be as capable as other more complex models, that it is more appropriate for dynamic simulations of intermediate duration events, and that it is better suited for optimal control usage. The model has been validated against data from only two experiments, and some deficiencies are apparent. Neither experiment completely explores all possible behaviors of the model: one experiment exercised more completely the force–velocity behavior and the other focused on recovery periods. In each case there are some model deficiencies regarding the unexercised behaviors. An obvious next step would be to conduct experiments which explore
3453
a wider range of recovery periods, target muscle speeds, and nonmaximal exertion. Successfully predicting new behaviors, and not just demonstrating the ability to generate a close fit, would provide the strongest validation. This could be done in a number of ways, but determining the model parameters for a test subject would always be the starting point. Validation could then be accomplished by using different recovery periods for a rest interval style test, at multiple tempos. A more complex experiment might involve a period of steady-state non-maximal effort inducing partial fatigue, followed by a period of maximal effort. If the subsequent maximal effort exercise was accurately predicted then the model would be further validated.
Conflicts of interest statement The authors have no conflicts of interest to report.
Acknowledgments The authors would like to thank Prof. George A. Brooks for discussions on human physiology, and Prof. David Hawkins for valuable advice on modeling. References Henneman, E., Somjen, G., Carpenter, D., 1965. Functional significance of cell size in spinal motoneurons. J. Neurophysiol. 19 (5), 560–580, URL: 〈http://www. neurosciences.us/courses/systems/SensoryMotor/henneman.pdf〉. Jones, E., Oliphant, T., Peterson, P., et al., 2001. SciPy: Open Source Scientific Tools for Python. URL: 〈http://www.scipy.org〉. Keller, J., 1974. Optimal velocity in a race. Am. Math. Mon. 81 (5), 474–480 URL: 〈http://www.jstor.org/stable/2318584〉. Liu, J.Z., Brown, R.W., Yue, G.H., 2002. A dynamical model of muscle activation, fatigue, and recovery. Biophys. J. 82 (May (5)), 2344–2359, URL: 〈http://www. pubmedcentral.nih.gov/articlerender.fcgi?artid=1302027&tool=pmcentrez& rendertype=abstract〉. Morton, R.H., 2009. A new modelling approach demonstrating the inability to make up for lost time in endurance running events. IMA J. Manag. Math. 20 (August (2)), 109–120, URL 〈http://imaman.oxfordjournals.org/cgi/doi/10.1093/imaman/ dpn022〉. Sih, B., Ng, L., Stuhmiller, J., 2012. Generalization of a model based on biophysical concepts of muscle activation, fatigue and recovery that explains exercise performance. Int. J. Sports Med. 33 (4), 258–267, URL 〈http://www.ncbi.nlm. nih.gov/pubmed/22403006〉. Tummers, B., 2006. DataThief III. URL 〈http://datathief.org〉. Vanhatalo, A., Doust, J.H., Burnley, M., 2008. Robustness of a 3 min all-out cycling test to manipulations of power profile and cadence in humans. Exp. Physiol. 93 (March (3)), 383–390, URL 〈http://www.ncbi.nlm.nih.gov/pubmed/17951327〉. Vanhatalo, A., Jones, A.M., 2009. Influence of prior sprint exercise on the parameters of the 'all-out critical power test' in men. Exp. Physiol. 94 (February (2)), 255–263, URL 〈http://www.ncbi.nlm.nih.gov/pubmed/18996948〉. Ward-Smith, A., 1999. The kinetics of anaerobic metabolism following the initiation of high-intensity exercise. Math. Biosci. 159 (June (1)), 33–45, URL 〈http:// linkinghub.elsevier.com/retrieve/pii/S0025556499000152〉. Xia, T., Frey Law, L.a., 2008. A theoretical approach for modeling peripheral muscle fatigue and recovery. J. Biomech. 41 (October (14)), 3046–3052, URL 〈http:// www.ncbi.nlm.nih.gov/pubmed/18789445〉.
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
A bioenergetic model for simulating athletic performance of intermediate duration Gilbert Gede n, Mont Hubbard University of California Davis, Mechanical and Aerospace Engineering, 1 Shields Ave, Davis, CA 95616, United States
art ic l e i nf o
a b s t r a c t
Article history: Accepted 17 September 2014
Simulating factors affecting human athletic performance, including fatigue, requires a dynamic model of the bioenergetic capabilities of the athlete. To address general cases, the model needs inputs, outputs, and states with a set of differential equations describing how the inputs affect the states and outputs as functions of time. We improve an existing phenomenological muscle model, removing unnecessarily fast dynamic behavior, adding force–velocity dependence, and generalizing it to task level activities. This makes it more suitable for simulating and calculating optimal strategies of athletic events of medium duration (longer than a sprint but shorter than a marathon). To examine the validity and limitations of the model, parameters have been identified from numerical fits to published experimental data. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Bioenergetics Model Fatigue Simulation Muscle Model Athletic Performance
1. Introduction When considering athletic events that involve racing against the clock such as a 5 km road race or cycling time trial, the following question inevitably arises: what is the optimal pacing strategy? Commonly practiced strategies range from maintaining a constant power output to finishing with a sprint. The question is especially difficult to answer in situations where the athlete's speed and power output vary over the duration of the event (e.g., a course with hills or wind). The first step in answering it consists of determining a model to represent the athlete's energetic dynamics in the event. Numerous authors have examined competitive athletic pacing. A runner's optimal strategy for different distances was calculated and the results were compared to the then current world records (Keller, 1974). This model used the concept of an oxygen debt. Anaerobic metabolism effects on energy production have also been investigated (Ward-Smith, 1999), but in a manner which is not applicable to longer events or events with a variable power output. The minimum time problem for running has been studied using the critical power model by Morton (2009) who found that the optimal solution is all-out for the entire race. All these efforts have limitations, either in fidelity or extensibility to more complex events. Simulating an athlete's performance in a general situation requires a dynamic model with a control input which governs
n
Corresponding author. E-mail addresses: [email protected] (G. Gede), [email protected] (M. Hubbard). http://dx.doi.org/10.1016/j.jbiomech.2014.09.017 0021-9290/& 2014 Elsevier Ltd. All rights reserved.
the athlete's output. While such dynamics models have existed for decades, a more recent approach (Liu et al., 2002) divides the muscle into three compartments: resting, active, and fatigued. This model has also been extended (Xia and Frey Law, 2008) to include multiple muscle fiber types and the control logic required to activate the fibers in a realistic fashion. More recently, a 4 compartment, single fiber type model was proposed which appears to faithfully replicate experimental behavior (Sih et al., 2012). However, if the goal is to perform optimal control calculations which can tell an athlete what power output, as a function of time, they should produce to achieve the minimum elapsed time in a cycling time trial, none of these models are appropriate. The models above all have significant weaknesses which prevent them from filling this role. To best answer the question, “what is the optimal strategy?”, we present a new phenomenological model which represents as much behavior as needed and no more. In this paper, we present a new dynamic model, designed for answering pacing strategy questions. It is an improvement of the existing phenomenological model of Xia and Frey Law (2008); it is more appropriate for general simulation and use in optimal control calculations. It does not attempt to replace current knowledge on muscle physiology, nervous system control, or decades of clinical experiments. Instead, it attempts to represent what an athlete can do, without considering the details of the physiological mechanisms of how this happens. Compared to current models, it removes unnecessary detail and adds important features. A limited validation using published experimental data is shown. Finally, interpretations of the identified physiological parameters are discussed.
G. Gede, M. Hubbard / Journal of Biomechanics 47 (2014) 3448–3453
3449
2. Methods A qualitative overview of the model is first presented, followed by a mathematical description, and the methods used for its validation.
2.1. Background Compartment theory has been used for modeling a single muscle (Liu et al., 2002; Xia and Frey Law, 2008). Each motor unit (MU) pool within the muscle is a collection of a finite number of similar fibers, but for modeling purposes it is convenient to treat it instead as continuous. The continuous pool is apportioned into three compartments: resting fraction Mr, active fraction Ma, and fatigued fraction Mf. The dynamics of the model describe the transfer between compartments: active muscle becoming fatigued, fatigued muscle recovering, and resting muscle becoming active. Fig. 1 shows a schematic of this process during and after a maximal effort exercise. We have taken the approach of Xia and Frey Law (2008) a step further by applying it to groups of muscles, and use a 3 compartment model with multiple MU pools to represent task-level activities such as cycling, as the model of Sih et al. (2012) did with success. The fundamentals of the 3 compartment model are unchanged: the amount of force produced is proportional to the amount of active muscle, the rate at which fatigue accrues is proportional to how much muscle is active, and the rate of recovery is proportional to the current level of fatigue. There are differences in the model we present however, in the rate of activation, elimination of redundant quantities, and the inclusion of the force–velocity curve, among other changes. Additionally, the original model was developed for examining the behavior of motor unit pools within a single muscle. The present model was formulated with the assumption it could also represent a group of muscles. Fig. 2 shows how the compartments are allocated within a single motor unit pool.
Fig. 1. A schematic representation of muscle transfer between compartments before, during, and after a maximal effort exercise. Initially the muscle is completely at rest. Once maximum effort starts, no muscle is at rest, and fatigue builds up. At the end of the effort period, recovery begins with active muscle being released and fatigued muscle recovering.
2.2. Model dynamics Unlike the Xia and Frey Law (2008) model, we assume there are only slow and fast fibers (but no intermediate ones). The quantities involved with the slow oxidative fiber MU pool are identified by a superscript o and the fast glycolytic fiber MU pool quantities by a superscript g. The differential equations for the single state variable for each pool, the fatigued muscle fraction Mf, are g
dM f dt
¼ M ga F g M gf Rg
ð1aÞ
o
dM f ¼ M oa F o M of Ro dt g
g
o
ð1bÞ o
where F , R and F , R are the fatigue and recovery rates for the fast and slow fiber types, respectively. As will be shown, only one differential equation per pool is required. This has been accomplished by removing mathematically redundant equations and instead partitioning the motor unit pools into fatigued and nonfatigued compartments (where it is assumed the non-fatigued compartment can be instantaneously apportioned between resting and active). This model, which has 3 compartments per pool, can be represented with only one state variable for each MU pool for two reasons. First, there is clearly a redundant state, due to the total muscle size constraint 1 ¼ Ma þ Mf þ Mr
ð2Þ
making one state dependent on the other two. Second, by not considering the rate of muscle activation, Ma no longer needs to be considered a state variable. The purpose of the present model is to simulate athletic events of medium duration (longer than a sprint but shorter than a marathon). We claim that the sub-second (much shorter than a second) time dynamics of muscle fiber recruitment are not relevant to athletic events of this and longer durations. If these sub-second dynamics are not considered, the nonfatigued muscle in a MU pool can be instantaneously split between the active and resting compartments. This allows for the control input to the system to simply be the amount of active muscle Ma, rather than the recruitment rate of active muscle. These new control inputs (Ma for each fiber type) have constraints. Within each pool the amount of active muscle must be non-negative, and it cannot be greater than the amount of non-fatigued (i.e. available for recruitment) muscle 0 r M ga r 1 M gf
ð3aÞ
0 r M oa r 1 M of
ð3bÞ
There is also a constraint imposed by the prioritization of fiber recruitment. Previous compartment models have used a logic sequence to enforce the correct order of fiber recruitment (Henneman et al., 1965). The following constraint
Fig. 2. A visualization of the compartments within a MU pool. The bar dividing the fatigued and non-fatigued compartments moves down at a rate which is the difference of the recovery ðM f RÞ and fatigue ðM a FÞ. The system control Ma is bounded between zero and the amount of non-fatigued muscle. The resting compartment Mr is simply the non-active, non-fatigued fraction.
equation enforces the same behavior: M ga M oa þ M of 1 ¼ 0
ð4Þ
The allowable space for the control inputs to the system is shown as a bold line in Fig. 3. This constraint enforces that the amount of active fast fiber can be nonzero only when all of the slow fibers are active. Using a constraint equation is necessary when applying optimal control techniques, because logic sequences introduce nonsmoothness, possibly preventing an optimal solution from being found. Outside optimal control applications, an equivalent logic sequence is acceptable. The final change from other compartment models is the incorporation of the force–velocity curve within the equations. For medium duration events an athlete is unlikely to spend time at either extreme of the force–velocity curve, and we claim that an exponential force–velocity approximation is sufficient. This allows the
3450
G. Gede, M. Hubbard / Journal of Biomechanics 47 (2014) 3448–3453
τ ¼ 1=ðF þ RÞ, is the time constant of the exponential decay to the steady-state force output, when the muscle is at maximum effort (see Fig. 1). The fiber type with the smaller time constant (shorter period) is the faster fiber, as it will take less time to become fatigued. We now have an expression (7) for the torque output of each fiber type as a function of time and crank speed assuming that the system parameters are known. A least squares algorithm will be used to compute parameters to fit the experimental data, scipy:optimize:leastsq (Jones et al., 2001). The eight model parameters (S0, c, F, and R for the two fiber types) are found by comparing the predicted torque based on the exact solutions of the differential equations to the experimental torque. The error is simply the difference between the measured and the predicted torque at each time. Because our model does not consider muscle activation time, we consider the experimental data only after the peak torque has been reached. 2.4. Experiments
Fig. 3. Allowable configuration space for the two inputs, Moa and Mga, is shown by the bold line. This configuration space enforces the prioritized recruitment behavior required by the Henneman et al. (1965) size principle. final muscle model output force to be expressed as F ¼ Sg0 M ga e c v þ So0 M oa e c g
ov
ð5Þ
where Sg0, cg and So0, co are the muscle size and force–velocity curve coefficients, for the fast g and slow o fibers respectively, v is the muscle speed (or an analog thereof, such as crank speed when cycling or ground speed when running, which is always positive), an exogenous input, and F is the muscle force output. 2.3. Validation methods The validation presented below consists of selecting appropriate experimental data, solving equations, and numerical parameter fitting. Because it is difficult to measure the individual activation amount of both slow and fast motor unit pools, we have limited the types of experiments used for validation. The selected experiments are maximal effort tests, in which the amount of active muscle will always be the maximum possible for both fiber types: M a ¼ 1 M f , i.e. Mr ¼0. By selecting maximal effort tests with recorded force or power outputs and recorded velocities (or analogs thereof), the parameter fitting process will identify Fo, Fg, Ro, Rg, co, cg, So0, and Sg0 in (1) and (5). However, the model in general is not limited to the maximal effort case. The experiments are task level (cycling) activities, involving collections of joint groups which are themselves collections of muscles. Thus far, the model has been intentionally described without a scale of application. For the task level activity of cycling, we attempt to show that this level of model abstraction can accurately represent the combined behavior of all the muscles involved. At this point, two substitutions in the equations will be introduced, in order to use them with cycling. The muscle force output F is instead substituted with crank torque T, and muscle velocity v with crank angular velocity (speed) ω. This represents a transformation from linear to rotational variables; aside from the change in symbols and units Eq. (5) is unaffected. Because the differential equations presented here are relatively simple, they can be solved exactly both during maximal effort tests and during rest. After solving the differential equations for Mof and Mgf , the amount of active muscle is found from M a ¼ 1 M f . The total torque output is the sum of the torque output of each fiber type. During a maximal effort test the solution to the differential equations for a single MU pool (one fiber type) is F F þ M f ðt 0 Þ M f ðtÞ ¼ ð6Þ e ðt t0 ÞðF þ RÞ F þR F þR
Two experiments suitable for validation have been published, both using a cycling ergometer where power was a quadratic function of crank speed: P ¼ r ω2 . The power dissipated P is the product of a resistance r and the speed ω squared. Unlike Coulomb friction-braked ergometers, this type of ergometer is associated with a gradual decrease in crank torque produced by the athlete over a maximum effort test. In contrast, friction-braked ergometers show an initial peak in torque (due to flywheel acceleration), followed by a constant torque (maintaining equilibrium against friction), then an abrupt end (inability to maintain the torque required to balance friction). One experiment (variable tempo) consisted of maximum effort tests at three resistive loads, with fully rested subjects (Vanhatalo et al., 2008). Power and speed were recorded while subjects performed with low, medium, and high resistances, giving high, medium, and low peak speeds and end-test speeds (see Fig. 4). A similar experiment (rest interval) used different subjects (see Fig. 6). It involved periods of maximum effort expenditure separated by rest intervals varying in duration (Vanhatalo and Jones, 2009). In each test fully rested subjects exercised at maximal effort for 30 s, followed by a rest period of either 2 or 15 min, followed by another maximal effort exercise for 3 min. For a fully rested athlete, it was assumed that the 30 s expenditure was identical during each test, and was identical to the first 30 s of a 3 min control test. Data from the figures in each publication were digitized using software (Tummers, 2006). In both experiments the figures displayed results averaged across many subjects. Individual subject data would be preferable for these validation analyses, but we believe there are still worthwhile insights to be gained from examining the averaged data.
3. Results For both experiments crank torque was computed by dividing the reported power by the reported crank speed. A least squares fit was used to compute a single set of eight model parameters (four for each MU pool) in each experiment. To show the ability of the model to work over a range of muscle speeds and recovery periods, parameters were found using the data from all tests. For
The output torque for each MU pool is then TðtÞ ¼ S0 e cωðtÞ ð1 M f ðtÞÞ
ð7Þ
When starting from a fully rested state there is no fatigued muscle and M f ðt 0 Þ ¼ 0. To find Mf at time t during a rest period, after exertion ending at t1, the differential equations are again solved to yield M f ðtÞ ¼ M f ðt 1 Þe Rðt t 1 Þ
ð8Þ
We also show two quantities which are calculated from the model parameters: the steady-state force fraction and the time constant. The steady-state force output for each fiber type is defined as G ¼ 1 ðF=ðF þ RÞÞ. This is the force that can be still produced by each fiber type when the athlete has reached the maximum sustainable level of fatigue. The time constant for each fiber type, defined as
Fig. 4. Crank speed from three tests with different ergometer resistances, leading to different tempos (Vanhatalo et al., 2008).
G. Gede, M. Hubbard / Journal of Biomechanics 47 (2014) 3448–3453
both experiments, the identified parameters which were either very large or very small in magnitude were rounded to infinity or 0, respectively. The difference between the exact output of the least squares solver and the rounded parameters was less than machine precision, justifying the use of the rounded parameters. The variable tempo experiment included power and speed from maximal effort cycling ergometer tests (Vanhatalo et al., 2008). Fig. 4 displays the speed ω in the tests, the exogenous input term used in parameter identification. Fig. 5 displays torque for the three tests, found by dividing the reported power by the speed. Also plotted are the fit lines showing the predicted model torque resulting from parameter identification. The RMS torque prediction error is 0.49 Nm. The variable tempo parameters are in Table 1. The rest interval experiment also measured power and speed from maximal effort cycling ergometer tests (Vanhatalo and Jones, 2009). Fig. 6a and b shows the calculated speed, derived from the reported power and maximum speed, using the equation defining the ergometer's power/speed relationship. Fig. 6c and d portrays the torque, calculated by dividing the recorded power by calculated speed. Also plotted in Fig. 6c and d are the fit lines for predicted torque from the identified parameters. Fig. 6a and c shows speed and torque from the 30 s effort; Fig. 6b and d plots speed and torque from the 3 min efforts. The RMS torque prediction error is 0.37 Nm. The rest interval parameters are also in Table 1.
4. Discussion In formulating the model, we have assumed that equations which govern the behavior of MU pools within a single muscle
Fig. 5. Experimentally calculated (data from Vanhatalo et al., 2008) and model fit torque for three maximal effort tests with different ergometer resistances.
Table 1 Parameters generated by least squares fit of experimental data for both the variable tempo experiment and the rest interval experiment. Parameter
S0 (N m) F ðs 1 Þ R ðs 1 Þ c ðs 1 Þ G τ (s)
Variable tempo test
Rest interval test
Slow
Fast
Slow
Fast
9.58 0.0062 0 0.0186 0 162
371 0.042 0.014 0.146 0.25 17.8
10.2 1.7E þ 5 Inf 0 1 0
42.3 0.017 0.012 0 0.407 34.8
3451
could also be used to describe the behavior of a muscle group or task level activity involving many muscle groups (which each have their own MU pools). Even though the original model (Xia and Frey Law, 2008) referred to individual muscles, we cannot claim our model represents individual muscles, because it has not been experimentally validated on this scale. Instead, we only examine the model's ability to replicate the behavior of task level activities. The success of this model in fitting the examined experimental behavior shows that it is apparently able to represent such an activity (cycling). For each experiment, the distinction between the slow and fast MU pools is made by examining the time constants. The pool which reaches full fatigue sooner is the fast pool, generally identified by a smaller time constant τ. Examining the parameters found from the variable tempo experiment (Table 1), we see that a relatively large isometric muscle size S0 has been found for the fast fibers, compared to that of the slower fiber MU pool. The faster muscle pool has a time constant of 17.8 s and a static torque capability of 371 Nm. The slower pool's time constant is much longer, 162 s, with a much smaller static torque capacity of 9.58 Nm. Admittedly, the variation in MU pool size does not seem to be reasonable: the value for S0 for the fast fiber MU pool is approximately 40 times that of the slower MU pool. While these static torque values appear to be significantly different than what might be expected physiologically, they must be viewed in the context of the force–velocity curve approximation. For the average crank speed in the standard target speed run, the actual torque, Se cω , is lower for each MU pool: 82.2 Nm for the fast pool, and 7.91 Nm for the slow pool. This is shown more clearly in Fig. 7, where the force velocity curve for each fiber type is plotted against the speed range for the experiments. Also indicated are the average crank speed and the time distribution of the crank speed. This indicates that the operational difference between S0 for the slow and fast fiber groups is not as large as the parameters would suggest. Discrepancies such as this are to be expected, because the model contains an approximation to the force velocity curve. Interestingly, when fitting parameters to a single one of the three variable tempo tests, the identified force–velocity curve parameter c and recovery coefficients R are very small. This shows that these model components, recovery and the force–velocity curve, come into play only when used with test data which exhibit recovery behavior or significantly different speed ranges. The rest interval experiment was conducted with different subjects from those in the variable tempo ones, and thus identical parameters should not be expected. In this experiment very small force–velocity coefficients were identified for both MU pools. Additionally, the slow MU pool showed very large fatigue and recovery rate coefficients. It is important to recognize though that these coefficient only appear as parts of terms in (6). Therefore the magnitude of the coefficients is less important than how they interact within that equation. The previously defined terms G and τ (steady-state force output and time constant) provide an interpretation of these ratios. In this experiment the fast MU pool has a time constant of 34.8 s, the slow MU pool has a time constant of essentially 0 s. At first glance, this seems to be inconsistent with the previous definition of fast and slow fibers. However, there is a steady-state force output of 100% for the slow pool. This indicates that the slow pool immediately reaches its steady-state value, which is equal to its unfatigued value. Its maximum force output does not change. Between the two experiments (variable tempo and rest interval) there is a similar size ratio between the slow and fast fibers. Over the operational speed ranges, the rest interval experiment shows a factor of 4 between fast and slow pool torque outputs, while the variable tempo experiment shows a factor of 10 between fast and slow pool torque outputs. The fast MU pool fatigue rate
3452
G. Gede, M. Hubbard / Journal of Biomechanics 47 (2014) 3448–3453
Fig. 6. Experimentally calculated speeds (a and b) and experimentally calculated and numerically fit torques (c and d) for 30 s and subsequent 3 min maximal effort expenditures, after different rest intervals (experimental data from Vanhatalo and Jones, 2009). Note that the 30 s data is assumed to be identical in both cases, as the athletes were completely rested before each test.
coefficients identified from the experiments are of the same order of magnitude: 0.042 vs. 0.017. Also, in both experiments the slow fibers are found to have little or no decrease in force production over the test duration; only the fast fibers fatigue appreciably. Considering all these points, there is some reassurance in the model's self-consistency across various test conditions. 4.1. Examination of reduced models
Fig. 7. Fast and slow MU pool force–velocity curves identified from the variable tempo experiment. Also indicated are the average speed (dashed line) and the distribution of speeds (shaded grey) over the experiment. There is a large difference between the fast and slow static torque parameters S0, 370 N m and 10 N m. As the static torque parameter appears in S0 M a e cω , their output is reduced by the force–velocity curve. In this case, over the main operating range (9–11 rad/s) the difference between product of the static torque parameter and the force velocity curve approximation of the slow and pools is 80 N m and 8 N m, a much smaller difference than the base values of S0 suggest.
Examining the results in Table 1 should raise concerns over the accuracy of a human energetics model which involves coefficients of 0 or infinity. In addition to the preceding discussion of the fit parameters, we have examined how the results of the fitting process are influenced by examining the maximal effort case; specifically how some model components combine, which leads to these unrealistic parameters. By examining the parameters when considering only the fullyrested case (i.e. M f ðt 0 Þ ¼ 0), the analytic solutions to the differential equations are greatly simplified. When there is no initial fatigue, the torque output equation can be rewritten in terms of G and τ; for a single MU pool, the torque output is ð9Þ TðtÞ ¼ S0 e cωðtÞ G þ ð1 GÞe t=τ Re-examining the parameters in Table 1, we can see how the identification of slow and fast fibers can be made. For the variable
G. Gede, M. Hubbard / Journal of Biomechanics 47 (2014) 3448–3453
tempo experiments, the slow fibers' torque output is TðtÞ ¼ 9:58e0:0186ωðtÞ e t=162 and the fast fibers' torque output is TðtÞ ¼ 371e0:146ωðtÞ 0:25 0:75e t=17:8 . Examining these equations shows that the fast fibers will clearly reach their steady-state output sooner. However, if these model parameters were used in a longer simulation, the model would become inaccurate, as it would suggest that all the slow fibers would be fatigued and only fast fibers would generate torque – something that is not consistent with physiology. One clear conclusion to draw from this is to generate the parameters from longer test sessions. For the rest interval experiments, we obtain the following equations: the slow MU pool torque is TðtÞ ¼ 10:2 and the fast MU pool torque is TðtÞ ¼ 42:3 0:407 þ 0:593e t=34:8 . The slow MU pool is identified as the pool which shows increased resistance to fatigue. In this case it shows no change at all! The fast MU pool does show signs of fatigue however. This raises the question: is there a need for the complete model that has been presented in this paper? The answer is unequivocally yes. It is important to remember that these submodels are the result of the numerical parameters being substituted into the solutions to the differential equations solved for the maximal effort case. There are no simpler forms of the differential equations which can be solved and manipulated into all the forms above. This reduced model analysis does reveal some model weaknesses. In many dynamic systems when the inputs are saturated the solutions to the differential equations become simplified. This is what has happened here. This also shows that the model needs more thorough parameter identification tests; the parameters shown for the maximal effort case are not completely translatable to the general case (i.e. non-maximal effort). Parameter identification studies are more effective when the system inputs are known. Because we have not considered the mechanisms for MU pool activation, the inputs in this case cannot be known exactly. This might be a potential weakness of the present model when attempting to find parameters in more general applications. 5. Conclusions and future work We have presented a model which reproduces the experimental dynamic energetic behavior of human subjects while cycling, displaying both fatigue and recovery behavior. Although further validation of the model is required, we believe it may be as capable as other more complex models, that it is more appropriate for dynamic simulations of intermediate duration events, and that it is better suited for optimal control usage. The model has been validated against data from only two experiments, and some deficiencies are apparent. Neither experiment completely explores all possible behaviors of the model: one experiment exercised more completely the force–velocity behavior and the other focused on recovery periods. In each case there are some model deficiencies regarding the unexercised behaviors. An obvious next step would be to conduct experiments which explore
3453
a wider range of recovery periods, target muscle speeds, and nonmaximal exertion. Successfully predicting new behaviors, and not just demonstrating the ability to generate a close fit, would provide the strongest validation. This could be done in a number of ways, but determining the model parameters for a test subject would always be the starting point. Validation could then be accomplished by using different recovery periods for a rest interval style test, at multiple tempos. A more complex experiment might involve a period of steady-state non-maximal effort inducing partial fatigue, followed by a period of maximal effort. If the subsequent maximal effort exercise was accurately predicted then the model would be further validated.
Conflicts of interest statement The authors have no conflicts of interest to report.
Acknowledgments The authors would like to thank Prof. George A. Brooks for discussions on human physiology, and Prof. David Hawkins for valuable advice on modeling. References Henneman, E., Somjen, G., Carpenter, D., 1965. Functional significance of cell size in spinal motoneurons. J. Neurophysiol. 19 (5), 560–580, URL: 〈http://www. neurosciences.us/courses/systems/SensoryMotor/henneman.pdf〉. Jones, E., Oliphant, T., Peterson, P., et al., 2001. SciPy: Open Source Scientific Tools for Python. URL: 〈http://www.scipy.org〉. Keller, J., 1974. Optimal velocity in a race. Am. Math. Mon. 81 (5), 474–480 URL: 〈http://www.jstor.org/stable/2318584〉. Liu, J.Z., Brown, R.W., Yue, G.H., 2002. A dynamical model of muscle activation, fatigue, and recovery. Biophys. J. 82 (May (5)), 2344–2359, URL: 〈http://www. pubmedcentral.nih.gov/articlerender.fcgi?artid=1302027&tool=pmcentrez& rendertype=abstract〉. Morton, R.H., 2009. A new modelling approach demonstrating the inability to make up for lost time in endurance running events. IMA J. Manag. Math. 20 (August (2)), 109–120, URL 〈http://imaman.oxfordjournals.org/cgi/doi/10.1093/imaman/ dpn022〉. Sih, B., Ng, L., Stuhmiller, J., 2012. Generalization of a model based on biophysical concepts of muscle activation, fatigue and recovery that explains exercise performance. Int. J. Sports Med. 33 (4), 258–267, URL 〈http://www.ncbi.nlm. nih.gov/pubmed/22403006〉. Tummers, B., 2006. DataThief III. URL 〈http://datathief.org〉. Vanhatalo, A., Doust, J.H., Burnley, M., 2008. Robustness of a 3 min all-out cycling test to manipulations of power profile and cadence in humans. Exp. Physiol. 93 (March (3)), 383–390, URL 〈http://www.ncbi.nlm.nih.gov/pubmed/17951327〉. Vanhatalo, A., Jones, A.M., 2009. Influence of prior sprint exercise on the parameters of the 'all-out critical power test' in men. Exp. Physiol. 94 (February (2)), 255–263, URL 〈http://www.ncbi.nlm.nih.gov/pubmed/18996948〉. Ward-Smith, A., 1999. The kinetics of anaerobic metabolism following the initiation of high-intensity exercise. Math. Biosci. 159 (June (1)), 33–45, URL 〈http:// linkinghub.elsevier.com/retrieve/pii/S0025556499000152〉. Xia, T., Frey Law, L.a., 2008. A theoretical approach for modeling peripheral muscle fatigue and recovery. J. Biomech. 41 (October (14)), 3046–3052, URL 〈http:// www.ncbi.nlm.nih.gov/pubmed/18789445〉.
