INTRODUCTION

TOP ABSTRACT INTRODUCTION DEVELOPMENT OF THE FATIGUE... METHODS RESULTS DISCUSSION REFERENCES

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.

View larger version (9K): [in this window] [in a new window]   Fig. 1.   Force responses to stimulation patterns. Data were collected from 1 typical subject's quadriceps femoris muscle under the nonfatigue condition. Vertical bars on the x-axis designate the stimulation pulses with pulse width of 600 µs. CFT25, CFT66, and CFT100, constant-frequency trains with interpulse intervals (IPIs) of 25, 66, and 100 ms; VFT50 and VFT80, variable-frequency trains with an initial doublet and subsequent IPIs of 50 and 80 ms; DFT155, stimulation pattern containing doublets throughout the train with IPI of 155 ms.

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

TOP ABSTRACT INTRODUCTION DEVELOPMENT OF THE FATIGUE... METHODS RESULTS DISCUSSION REFERENCES

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

(1)

with
and

(2)

The 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 Ca²⁺-troponin complex (C_N, unitless), in which R_i (unitless) accounts for the nonlinear summation of the Ca²⁺ transient in single muscle fibers when stimulated with two closely spaced pulses (16). This summation enhancement decays with IPI t_i  t_i  1 (ms). Thus, R_i  1 for a pulse that occurs a long time after the preceding pulse, and R_i  R₀ 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: C_N/(1 + C_N).

The force model is governed by only five free parameters: R₀, A, _c, ₁, and ₂ (Table 1). In our previous study, _c was fixed at 20 ms for nonfatigued and fatigued muscles and ₁ 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²) 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 Ca²⁺ transient with fatigued muscles has been observed by Westerblad and Allen (30), which might be due to depressed Ca²⁺ uptake (31), reduced rate of Ca²⁺ release (18), and elevated resting intracellular Ca²⁺ 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 ₁ and ₂. Conversely, the prolonged relaxation of the force (controlled by ₁ and ₂ in the force model) could also be modeled by the slower Ca²⁺ dynamics with the fatigued muscle. For simplicity, in this study for fatigued muscles, therefore, we fixed ₁ and ₂ at the values obtained in the nonfatigued condition and fit the values of the other three parameters (A, R₀, 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 of A, R₀, 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

(3)

(4)

(5)

This fatigue model is governed by four parameters: _A (ms²), _R0 (ms¹ · N¹), _c (N¹), 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 (see METHODS and Fig. 3).

	METHODS

TOP ABSTRACT INTRODUCTION DEVELOPMENT OF THE FATIGUE... METHODS RESULTS DISCUSSION REFERENCES

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.

Stimulation Protocols

View larger version (15K): [in this window] [in a new window]   Fig. 2.   Schematic representation of the experimental protocols. Each subject received 6 sessions of testing, and each session started with the same potentiation protocol and continued with 1 of the 6 fatigue protocols.

Subjects

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² 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² 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 _c was set at 20 ms and ₁ was identified directly from the force record as described previously (15). The remaining three force model parameters (A, R₀, and ₂) were identified using the following objective function

(6)

where G₁ 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. 3A).

View larger version (24K): [in this window] [in a new window]   Fig. 3.   Schematic representations of the procedure used to parameterize the force model for the nonfatigue condition (A) and during fatigue testing (B), to parameterize the fatigue model (C), and to predict the force responses to each fatigue protocol (D). See Table 1 for definition of abbreviations.

(7)

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

(8)

where 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; t_q 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 for A, R₀, and _c, respectively) predicted by Eqs. 3-5 at time t_q, and is a function of _fat, the corresponding _j, and F(t_q) (the experimental force at time t_q). P_j^calc is the jth parameter value at time t_q calculated as described above. This objective function was minimized using CFSQP (23) (Fig. 3C).

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. 3D).

Data Analysis

View larger version (13K): [in this window] [in a new window]   Fig. 4.   Example of predicted forces to VFT50-CFT100 combination (V50-C100) fatigue protocol when the fatigue model was parameterized by the V50-C100 fatigue protocol for a typical subject. Force model parameter values were calculated (individual data points in A-C), and then fatigue model parameter values were identified (see values for A, R0, and c in A-C). The fat for this subject was 46.9 s. After the fatigue model was parameterized, force model parameter values (solid lines in A-C) and forces (D) were predicted given only initial force model parameter values (Arest = 11.9, R0,rest = 0.12, c,rest = 20) and the stimulation protocol (V50-C100).

View larger version (13K): [in this window] [in a new window]   Fig. 5.   Example of predicted forces to DFT155 fatigue protocol when the fatigue model was parameterized by the V50-C100 fatigue protocol for a typical subject (same subject as in Fig. 4). Force model parameter values were calculated (individual data points in A-C). After the fatigue model was parameterized as described in Fig. 4, force model parameter values (solid lines in A-C) and forces (D) were predicted given only initial force model parameter values (Arest = 11.0, R0,rest = 0.45, c,rest = 20) and the stimulation protocol (DFT155).

View larger version (21K): [in this window] [in a new window]   Fig. 6.   Experimental and predicted peak forces for a typical subject (same subject as in Fig. 4). The fatigue model was parameterized using data collected when the muscle was fatigued by the V50-C100 fatigue protocol (Fig. 4). The peak forces produced in response to the V50-C100 (A), CFT25 (B), CFT66 (C), CFT100 (D), VFT80 (E), and DFT155 (F) fatigue protocols were plotted. Although V50-C100 trains were included within each protocol, responses to these trains are not plotted.

View larger version (28K): [in this window] [in a new window]   Fig. 7.   Group (n = 6) R2 values (means ± SE) between experimental and predicted peak forces in response to each train for each fatigue protocol when the fatigue modeled was parameterized with force responses to the CFT25 (A), CFT66 (B), CFT100 (C), V50-C100 (D), VFT80 (E), and DFT155 (F) fatigue protocols. Arrows, model's prediction for the fatigue protocol that was used to parameterize the model.

	RESULTS

TOP ABSTRACT INTRODUCTION DEVELOPMENT OF THE FATIGUE... METHODS RESULTS DISCUSSION REFERENCES

When the fatigue model was parameterized with the V50-C100 protocol, it accurately predicted the changes in force model parameter A 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₀ appeared to change in a sigmoidal fashion, the model predicted that R₀ changed 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² values between the predicted and experimental peak forces were >0.9, except for the CFT100 (R² = 0.84) and DFT155 protocols (R² = 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%).

View larger version (25K): [in this window] [in a new window]   Fig. 8.   Group (n = 6) percent differences (means ± SE) between experimental and predicted peak forces at the end of fatigue protocol when the fatigue modeled was parameterized with force responses to the CFT25 (A), CFT66 (B), CFT100 (C), V50-C100 (D), VFT80 (E), and DFT155 (F) fatigue protocols. Arrows, model's prediction for the fatigue protocol that was used to parameterize the model.

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² = 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² data, 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).

View larger version (16K): [in this window] [in a new window]   Fig. 9.   Group (n = 6) percent decline (mean ± SE) in peak forces of fatigue-producing trains between experimental and predicted data. Predicted data were obtained by the fatigue model when the model was parameterized using the V50-C100 protocol. Paired tests were used to compare predicted and experimental percent declines. *Comparison was significant, P0.05.

	DISCUSSION

TOP ABSTRACT INTRODUCTION DEVELOPMENT OF THE FATIGUE... METHODS RESULTS DISCUSSION REFERENCES

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 in R₀ was not a simple exponential function of time, as we assumed in the fatigue model (Eq. 4). R₀ characterizes the magnitude of the nonlinear summation in Ca²⁺ produced by subsequent stimuli (for additional details see Ref. 15). The over- or underestimation of R₀ 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 and 5, the model overestimated R₀ 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 of R₀ in the fatigue model, such as adding another differential equation to allow a second (shorter) time constant to capture the changes in R₀ 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² of 0.90 and a percent difference of 9% between the predicted and experimental peak forces (Figs. 7 and 8). Because the within-subject experimental variance in the fatigue profiles was 7% (R² = 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.