Fatigue is a major limitation to the clinical application of functional electrical stimulation. The activation pattern used during electrical stimulation affects force and fatigue. Identifying the activation pattern that produces the greatest force and least fatigue for each patient is, therefore, of great importance. Mathematical models that predict muscle forces and fatigue produced by a wide range of stimulation patterns would facilitate the search for optimal patterns. Previously, we developed a mathematical isometric force model that successfully identified the stimulation patterns that produced the greatest forces from healthy subjects under nonfatigue and fatigue conditions. The present study introduces a four-parameter fatigue model, coupled with the force model that predicts the fatigue induced by different stimulation patterns on different days during isometric contractions. This fatigue model accounted for 90% of the variability in forces produced by different fatigue tests. The predicted forces at the end of fatigue testing differed from those observed by only 9%. This model demonstrates the potential for predicting muscle fatigue in response to a wide range of stimulation patterns.
- muscle fatigue
- functional electrical stimulation
- stimulation pattern
- stimulation frequency
functional electrical stimulation (FES) is the use of electrical stimulation to activate muscles to produce functional movements. FES has been used primarily to restore movement and function in patients with neurological disorders (1, 10, 22, 25). However, muscles are recruited differently when stimulated by FES from when they are activated by the central nervous system. During FES-induced contractions, motor units are activated synchronously, whereas during voluntary contraction, they are activated asynchronously (14,17). In addition, when stimulated by FES, the motor unit recruitment order differs from that seen during voluntary contractions (5). As a result, muscles fatigue more rapidly during FES than during contractions activated by the central nervous system. This rapid fatigue is a major limitation to the clinical application of FES (1, 22, 25).
It has been shown that the mean frequency and the pattern of pulses markedly affect force production and fatigue (3, 9). One approach to reducing fatigue is to identify the stimulation pattern that minimizes fatigue and maximizes force. Traditionally, during FES, skeletal muscles are activated with constant-frequency trains (CFTs), where the pulses within each train are separated by regular interpulse intervals (IPIs; Fig. 1). One stimulation pattern that has been shown to produce more force from a fatigued muscle than any CFT is a variable-frequency train (VFT) that begins with a doublet (i.e., 2 pulses separated by a brief IPI) (7). Recently, a stimulation pattern containing doublets throughout the train (DFT) was found to produce greater forces than CFTs and VFTs from fatigued human muscles (4, 15). However, the pattern that maximizes the force varies with the physiological conditions of the muscle, such as level of fatigue (15) or muscle length (24), and varies from person to person (15, 21). In addition, although VFTs have been shown to produce greater forces from fatigued muscles, they also produce greater fatigue than CFTs (8, 9). Moreover, muscle fatigue is a complex physiological phenomenon, and the underlying mechanisms are not well understood (2, 19). Numerous tests and testing sessions would be needed to identify the stimulation pattern that produces the most force and least fatigue for each patient. Mathematical models that can predict muscle forces and fatigue can speed up this process.
Several models have been developed to study muscle fatigue, including nonphysiological analytic models (11, 26, 28) and models that used physiological information such as myoelectric [electromyogram (EMG)] (12, 26) or metabolic (e.g., NMR) measurements (20, 27) as fatigue indexes added to muscle force models. All nonphysiological analytic models have only tested their ability to fit the forces produced during a fatigue test (11, 26, 28). Although studies have supported a good correlation between EMG amplitude and the forces during fatigue (12, 26), the predictive ability of an EMG-based fatigue model has never been demonstrated. Mizrahi and colleagues (26) developed a three-time-constant analytic function to curve fit the isometric forces produced by 20-Hz continuous stimulation (180 s). The model produced an accurate fit. They also found a strong exponential relationship (r = 0.986) between the force and the peak-to-peak amplitude of the EMG M wave during fatigue. They were not able, however, to predict forces using EMG data when the muscle was fatigued by a different fatiguing protocol (26). One of the limitations in modeling fatigue from an EMG is the difficulty in obtaining clean, artifact-free EMG data, because the stimulation currents are orders of magnitude larger than those produced by muscle action potentials (12). In addition, factors such as muscle fiber composition and stimulation history of the muscle complicate the force-EMG relationships, making it less practical to predict forces during fatigue from the EMG signals alone (26).
Mizrahi and colleagues (20, 27) used FES to develop a series of musculotendon models of fatigue and recovery of the quadriceps femoris muscles of paralyzed individuals. The 5-element musculotendon model needed 17 non-muscle-specific parameters taken from the literature and 5 muscle-specific parameters calculated from the femur length and circumference of the thigh of each subject. In addition, another seven parameters needed for the fatigue index were obtained by fitting the fatigue equations to pH data during fatigue and recovery taken from the literature. This model produced fairly good fits for isometric forces during 180-s fatigue tests at four different angles but did not include stimulation pattern or frequency as inputs; therefore, it is not capable of predicting the force responses to a variety of activation frequencies.
Previously, we developed a force model that accurately predicted the force responses to a wide range of stimulation patterns when muscles were nonfatigued or fatigued (15). We noticed that when the muscles were fatigued the parameter values of this force model differed significantly from those when the muscles were nonfatigued (15). Our initial approach to modeling fatigue is to model the changes in our force model parameter values during muscle fatigue. The purpose of this study is to present a physiologically based fatigue model that predicts changes in forces during repetitive activation of skeletal muscles. The modifications of our previous force model are presented first, and the theoretical development of the fatigue model is described immediately thereafter. Next the predictions of the model and comparison with the experimental data are evaluated.
DEVELOPMENT OF THE FATIGUE MODEL
Previously, we developed a mathematical model that successfully predicted isometric muscle forces to a wide range of stimulation patterns (15). This force model has two differential equations Equation 1with and Equation 2The symbols used in Eqs. 1 and 2 are defined in Table 1. Briefly, Eq. 1 represents the dynamics of the rate-limiting step leading to the formation of a Ca2+-troponin complex (CN, unitless), in which Ri (unitless) accounts for the nonlinear summation of the Ca2+ transient in single muscle fibers when stimulated with two closely spaced pulses (16). This summation enhancement decays with IPIti −ti − 1 (ms). Thus,Ri ≈ 1 for a pulse that occurs a long time after the preceding pulse, and Ri ≈R 0 for a pulse that occurs immediately after the preceding pulse. Equation 2 represents the development of mechanical force (F, in N), which is driven by force-producing cross-bridges, modeled by a Michaelis-Menten term: CN/(1 + CN).
The force model is governed by only five free parameters:R 0, A, τc, τ1, and τ2 (Table 1). In our previous study, τc was fixed at 20 ms for nonfatigued and fatigued muscles and τ1 was calculated using experimental forces (15). The other three parameter values were identified by fitting the model to the measured forces produced by stimulating the muscle with a combination of two trains (15). For nonfatigued muscles, the two trains that best parameterized the model were two VFTs: one with an initial 5-ms IPI followed by constant IPIs of 10 ms (VFT10) and the other with an initial 5-ms IPI followed by constant IPIs of 100 ms (VFT100). For fatigued muscles, the train combination was a CFT with IPIs of 100 ms (CFT100) and VFT100.
The major shortcoming of this force model is that it needed two different train combinations to parameterize the model for nonfatigued and fatigued muscles, making it impossible to model the fatigue process and to predict muscle forces under different levels of fatigue. Systematic search for the best train or train combination to parameterize the model under nonfatigue and fatigue conditions was conducted in a pilot study (unpublished observations). We found the combination of a VFT with an initial 5-ms IPI followed by four IPIs of 50 ms (VFT50) and CFT100 (i.e., a V50-C100 combination) to be the best. The V50-C100 combination improved correlation coefficients (R 2) between the predicted and experimental forces from 0.92 (old combination) to 0.95 (new combination) for nonfatigued muscles and from 0.89 to 0.92 for fatigued muscles.
As with any other mathematical representations of muscle force, this force model was inevitably based on some nonphysiological assumptions such as the same value of τc (20 ms) for nonfatigued and fatigued muscles. A prolonged Ca2+ transient with fatigued muscles has been observed by Westerblad and Allen (30), which might be due to depressed Ca2+ uptake (31), reduced rate of Ca2+ release (18), and elevated resting intracellular Ca2+concentration, however (13). This observation suggests that fixing τc at 20 ms in our previous study for fatigued muscles may not be appropriate physiologically, although doing so may not have affected the model's predictive ability because of the compensation effect from parameters τ1 and τ2. Conversely, the prolonged relaxation of the force (controlled by τ1 and τ2 in the force model) could also be modeled by the slower Ca2+ dynamics with the fatigued muscle. For simplicity, in this study for fatigued muscles, therefore, we fixed τ1 and τ2 at the values obtained in the nonfatigued condition and fit the values of the other three parameters (A, R 0, and τc) to the experimental forces. Further testing of the V50-C100 combination showed that it was still the best combination for parameterizing the force model under nonfatigue and fatigue conditions.
In addition, we found significant changes in the values ofA, R 0, and τc as the muscle made the transition from nonfatigue to fatigue conditions. Our approach to modeling fatigue was to develop a model that could modify the values of the force model parameters during fatigue. Studies have shown that, during sustained isometric contractions, the rate and amount of fatigue are related to the force-time integral produced by the muscle in response to the train used to fatigue the muscle (9). Therefore, we used instantaneous force as the driving function in the fatigue model.
We assume that there is one time constant (τfat) characterizing the rate of recovery. The fatigue model equations are Equation 3 Equation 4 Equation 5This fatigue model is governed by four parameters: αA (ms−2), αR0(ms−1 · N−1), ατc (N−1), and τfat (ms). A rest,R 0,rest, and τc,rest are the values when muscles are not fatigued. F is the isomeric force. Experimental forces were used in Eqs. 3-5 during parameter identification, but once the fatigue model was parameterized, forces predicted by Eqs. 1 and 2 were used to test the ability of the model to predict fatigue (seemethods and Fig. 3).
Experiments were conducted on human quadriceps femoris muscles of healthy subjects. The experimental setup was similar to that described previously (6). Briefly, subjects were seated on a force dynamometer with their hips flexed to ∼75° and their knees flexed to 90°. Two stimulating electrode pads were placed over the quadriceps muscle. The force transducer was positioned against the anterior aspect of the leg, proximal to the lateral malleolus. The subject's maximum voluntary isometric contraction (MVC) was first determined using a burst superimposition technique (6). The stimulus intensity was then set to elicit a force equal to 20% of the subject's MVC obtained during the first testing session. Isometric force data from the dynamometer's force transducer were digitized at a rate of 200 Hz.
Testing began by first potentiating the muscle with 35, 12-pulse CFTs with IPIs of 70 ms (CFT70; Fig. 2). One train was delivered every 5 s. Pilot work showed that this protocol maximally potentiated the muscle without causing a measurable decline in forces. Two seconds after the potentiation trains, a six-pulse VFT50 followed by 2 s of rest and then a six-pulse CFT100 (i.e., the V50-C100 pair that we previously showed was best for parameterizing the force model) were delivered to the muscle. These two trains were used to establish the initial, nonfatigue condition for the fatigue model. Two seconds after the V50-C100 pair, the fatigue protocol started. Different testing sessions used different fatigue protocols. Each train in the fatigue protocols had six pulses, and there was a 500-ms rest time between trains; so the only differences between stimulation trains were the pulse frequency and the train duration. One fatigue protocol had 75 V50-C100 pairs to fatigue the muscle (150 trains in total). Each of the other five fatigue protocols contained 10 cycles of trains; each cycle contained 13 fatigue-producing trains followed by the V50-C100 pair, giving a total of 150 trains. For each of these five fatigue protocols, the fatigue-producing trains were constant-frequency trains with IPIs of 25 ms (CFT25), 66 ms (CFT66), or 100 ms (CFT100), VFTs with an initial 5-ms doublet followed by constant IPIs of 80 ms (VFT80), or DFTs with a 5-ms doublet as first, third, and fifth IPIs and second and fourth IPIs of 155 ms (DFT155). The fatigue-producing trains had frequencies of 10 Hz (CFT100) to 40 Hz (CFT25) to cover the wide range of discharge frequencies reported during physiological activations of human quadriceps femoris muscle (29). The VFT80 and DFT155 were chosen because they had the same mean frequency as the CFT66 (15 Hz) but produced different force profiles (Fig. 1). The V50-C100 pairs included during the fatigue protocol were used to monitor changes in the force model parameters with time.
Eight healthy subjects were recruited to test the ability of the fatigue model to fit and predict forces during repetitive activation by stimulation trains with a wide range of frequencies and patterns. Each subject signed an informed consent form approved by University of Delaware Human Subject Review Board.
Three of the eight subjects were tested to study the variance of the experimental data. Each of these three subjects was tested with a CFT66 fatigue protocol (see Stimulation Protocols) on three different days. Peak forces were calculated for the force responses to each stimulation train throughout the testing. Fatigue profiles for peak forces were obtained by plotting the peak forces of each fatigue-producing train against their train number.R 2 values were used to compare the fatigue profiles between days (i.e., day 1 vs. day 2, day 2 vs. day 3, and day 1 vs. day 3). The range of R 2 values was 0.82–0.99 (mean 0.93). Next, the peak forces of the last 10 fatigue-producing trains were averaged to determine the fatigued forces for each day. For each subject, differences in the fatigued forces between days were averaged and then normalized to the averaged fatigued forces to obtain the percent difference. The experimental data showed that the percent differences in fatigued force between days ranged from 5.5 to 13.4% (mean 8.9%).
Six of the eight subjects were used to test the fatigue model by measuring the fatigue produced by six different stimulation patterns. Each subject was tested with all six patterns, only one pattern was tested during a session, and sessions were separated by ≥48 h. The order of sessions was randomized for each subject. The data were used to test the fatigue model.
Testing the Fatigue Model
Calculating the force model parameter values using experimental forces.
The force responses to the V50-C100 pair at the end of the potentiation protocol and during the fatigue protocols were used to calculate the sets of force model parameter values throughout the testing (76 sets of force model parameter values for VFT50-CFT100 fatigue protocol; 11 sets for the other 5 fatigue protocols). The V50-C100 pair delivered at the end of the potentiation protocol was used to calculate the initial, nonfatigued force model parameter values. Parameter τcwas set at 20 ms and τ1 was identified directly from the force record as described previously (15). The remaining three force model parameters (A, R 0, and τ2) were identified using the following objective function Equation 6 where G 1 was minimized using CFSQP (23), which employs feasibility set techniques coupled to sequential quadratic programming to identify the optimum values for the free variables numerically in an objective function (Fig.3 A).
To identify the force model parameter values during fatigue, an objective function similar to Eq. 6 was used, except τ1 and τ2 were fixed at the nonfatigued values and A, R 0, and τc were allowed to change (Fig. 3 B). The force model was fitted with the force responses to the V50-C100 pairs during fatigue by use of the following objective function Equation 7 G 2 was also minimized by CFSQP (23). There were 75 sets of (A,R 0,τc) for the V50-C100 fatigue protocol and 10 for each of the other 5 fatigue protocols.
Parameterizing the fatigue model.
Unlike previous physiology-based fatigue models that used EMG or metabolic changes as fatigue indexes, our fatigue model monitored the changes of our force model parameter values during fatigue (Eqs.3-5 ). To parameterize the fatigue model, the differences between the predicted and calculated force model parameter values were first normalized to the calculated force model parameter values. This normalization was performed to allow unbiased, concurrent evaluation of all three parameters. The objective function was then minimized Equation 8where q designates the number of the sets of force model parameter values and ranges from 1 to 76 for the V50-C100 protocol and from 1 to 11 for the other protocols;tq is 0 for the first set of force model parameter values and 855 ms after the start of the qth V50-C100 pair, which is the midpoint between the start of the V50 and the end of the C100; P j pred is the jth parameter value (j = 1, 2, and 3 forA, R 0 , and τc, respectively) predicted by Eqs. 3-5 at time tq , and is a function of τfat, the corresponding αj, and F(tq ) (the experimental force at time tq ).P j calc is the jth parameter value at time tq calculated as described above. This objective function was minimized using CFSQP (23) (Fig. 3 C).
Predicting forces by use of the parameterized fatigue model.
For each subject, the fatigue model was parameterized separately using the forces from each of the six fatigue protocols, resulting in six sets of fatigue model parameter values. Next, for each set of fatigue model parameter values, given the stimulation sequence of the fatigue protocol and initial nonfatigued force model parameter values, the forces produced in response to each train of each fatiguing protocol (including the protocol used to determine that set of fatigue model parameter values) were predicted using the force and fatigue models to integrate forward (Fig. 3 D).
For each of the six sets of fatigue model parameter values obtained (i.e., 1 from each fatigue protocol), we tested the ability of the fatigue model to predict the force responses to all six fatigue protocols, including the protocol used to parameterize the model. The evaluation of the fatigue model was begun by comparing subjectively the predicted with the calculated force model parameter values and by comparing the predicted with the experimental force responses to each stimulation train of the fatigue protocol (Figs.4 and 5). Then, the predicted and experimental peak force responses to each train of the fatigue protocol were plotted (Figs.6 and 7).R 2 values between the predicted and experimental forces were then calculated to compare the peak forces produced by each protocol. Next, the percent differences between the predicted and experimental peak forces at the end of fatigue test were compared. This measurement was included because R 2 is not sensitive to an offset between the predicted and experimental fatigue profiles. Finally, we compared the predicted with the experimental percent decline in forces produced by each fatigue protocol. Pairedt-tests were used to compare the predicted with the experimental percent decline in forces. Comparisons were significant ifP≤0.05.
When the fatigue model was parameterized with the V50-C100 protocol, it accurately predicted the changes in force model parameterA and τc (Figs. 4 and 5) when the muscles were fatigued by use of a wide range of fatigue protocols. However, although the data showed that R 0 appeared to change in a sigmoidal fashion, the model predicted that R 0changed exponentially. The fatigue model, parameterized with the V50-C100 protocol, generally predicted the force responses quite well, except for the beginning of the CFT25 and DFT155 protocols (Fig. 6).
When the fatigue model was used to predict the force responses to the same protocol that was used to parameterize the fatigue model (e.g., predict the force responses to CFT25 protocol when the model was parameterized by the CFT25 protocol), all the R 2values between the predicted and experimental peak forces were >0.9, except for the CFT100 (R 2 = 0.84) and DFT155 protocols (R 2 = 0.78; Fig. 7). Similarly, all the predicted peak forces at the end of fatigue testing were within 10% of the experimental peak forces, except for the CFT25 protocol (Fig. 8). The smallest percent differences were with the CFT100 and V50-C100 protocols (∼5%).
When the fatigue model was used to predict force responses to the protocols other than that used to parameterize the fatigue model (e.g., predict force responses to CFT100, CFT66, V50-C100, VFT80, and DFT155 protocols when the model was parameterized by the CFT25 protocol), the model generally accounted for >80% variance of the experimental forces for most protocols (Fig. 7). The fatigue model, however, did not do so well in predicting the force responses to the DFT155 protocols (R 2 = 0.52–0.79). The differences between the predicted and experimental peak forces at the end of the fatigue test were more variable than the R 2data, with results that ranged from 6% to 32% (Fig. 8). Overall, when the model was parameterized with the V50-C100 protocol, the fatigue model produced the smallest percent differences (9%).
Because the V50-C100 protocol was best at parameterizing the fatigue model among the six protocols tested, it was used to test the model's ability to predict the percent decline in peak forces during the fatigue test. When the V50-C100 protocol was used to parameterize the fatigue model, the predicted percent declines of the peak forces were within 7% of the experimental peak forces for all the fatigue protocols except the DFT155 protocol (46%) and the C100 trains in the V50-C100 protocol (15%). Significant difference was observed only for the DFT155 protocol (Fig. 9).
The fatigue model presented in this study was tested using data from human quadriceps femoris muscles with protocols consisting of trains with a wide range of frequencies (10–40 Hz) and activation patterns (CFT, VFT, and DFT). In general, subjective (Figs. 4-6) and objective evaluations (Figs. 7-9) suggested that the fatigue model predicted well the force responses to fatigue protocols with different stimulation patterns on different days.
Among the six protocols tested, the V50-C100 protocol was the best protocol to parameterize the fatigue model. We believe this was because this protocol gave the model the most information about changes in force. The V50-C100 protocol contained 75 pairs of the V50-C100 combination, giving the fatigue model 76 sets of force model parameter values to model. In contrast, there are only 11 pairs of the V50-C100 combination in the other protocols.
For the three force model parameters that were allowed to change during fatigue, the model fitted and predicted the changes in A and τc well (Figs. 4 and 5). However, the change inR 0 was not a simple exponential function of time, as we assumed in the fatigue model (Eq. 4 ).R 0 characterizes the magnitude of the nonlinear summation in Ca2+ produced by subsequent stimuli (for additional details see Ref. 15). The over- or underestimation ofR 0 by the model has a greater effect on the force responses to those stimulation trains containing closely spaced pulses, such as the CFT25, VFT80, and DFT155. As shown in Figs. 4 and5, the model overestimated R 0 in the beginning of the fatigue protocol, producing an increase in the predicted forces at the beginning of the fatigue tests for the VFT50 in the VFT50-CFT100 protocol and the DFT155 protocol (Fig. 6). This increase in force predicted by the model was most pronounced in the DFT155 protocol, because this pattern contained three doublets. In addition, we also observed that the fatigue model overestimated the force response to the first train in the CFT25 protocol and underestimated the forces to the first trains in the VFT80 and DFT155 protocols (Fig. 6). Because the fatigue model used the session's nonfatigued force model parameter values when the initial forces for that session were predicted, the over- or underestimation seen at the onset of the fatigue protocol was due to the inaccuracy of predictions in the force model. Modification of the force model and a more sophisticated modeling ofR 0 in the fatigue model, such as adding another differential equation to allow a second (shorter) time constant to capture the changes in R 0 at the beginning of fatigue, may be necessary to obtain better predictive ability in high-frequency trains or those containing doublets.
Previous studies on modeling muscle fatigue have not had stimulation frequency or pattern as one of the inputs and, therefore, were not designed to identify stimulation patterns to reduce fatigue during FES (11, 26, 28). The fatigue model presented in this study was governed by only four free parameters and was able to predict the forces of human quadriceps femoris muscles during repetitive activation with different fatigue protocols on different days, given only the stimulation protocol and the initial conditions of the muscle. When the model was parameterized by the best protocol (V50-C100), on average, it produced an R 2 of 0.90 and a percent difference of 9% between the predicted and experimental peak forces (Figs. 7 and8). Because the within-subject experimental variance in the fatigue profiles was 7% (R 2 = 0.93) and the percent difference in experimental forces at the end of the fatigue testing was 9% (see Subjects), this fatigue model was considered successful. In addition, the fatigue model also successfully predicted the effect of stimulation frequency and pattern on fatigue. Consistent with the experimental data, the model predicted that the CFT25 protocol would produce the most fatigue among the six protocols tested (Fig. 9). Also, among the three CFTs (CFT25, CFT66, and CFT100), the model predicted that the shorter the IPI, the greater the amount of fatigue; this was the same trend shown by the experimental data. Finally, for the three protocols that had the same mean frequency (CFT66, VFT80, and DFT155), the predicted and experimental data showed that the DFT was the least fatiguing train type.
In conclusion, this study is the first to develop a mathematical model that predicts isometric forces from human muscles during repetitive activation by different activation protocols. Unlike previous approaches to modeling muscle fatigue, we modeled the changes in our force model parameter values during fatigue, which led to the changes in the forces. The same approach could be applied to modeling fatigue during nonisometric muscle contractions once a successful nonisometric force model is developed. This model, although simple, appears very promising in predicting isometric forces across days, stimulation protocols, and levels of fatigue. Muscle fatigue depends on many factors, including stimulation parameters such as frequency, duty cycle, and activation patterns (2, 8, 21). We envision that mathematical models similar to that presented in this study could help elucidate the mechanisms of fatigue and facilitate the identification of stimulation patterns that minimize fatigue.
We thank Dr. William Rose for helpful comments concerning an earlier version of the manuscript.
This study was supported by National Institute of Child Health and Human Development Grant HD-36797 and the Office of Graduate Studies of the University of Delaware.
Address for reprint requests and other correspondence: S. A. Binder-Macleod, Dept. of Physical Therapy, 301 McKinly Laboratory, University of Delaware, Newark, DE 19716 (E-mail:).
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
- Copyright © 2000 the American Physiological Society