Abstract
Because muscles must be repetitively activated during functional electrical stimulation, it is desirable to identify the stimulation pattern that produces the most force. Previous experimental work has shown that the optimal pattern contains an initial highfrequency burst of pulses (i.e., an initial doublet or triplet) followed by a low, constantfrequency portion. Pattern optimization is particularly challenging, because a muscle's contractile characteristics and, therefore, the optimal pattern change under different physiological conditions and are different for each person. This work describes the continued development and testing of a mathematical model that predicts isometric forces from fresh and fatigued muscles in response to brief trains of electrical pulses. By use of this model and an optimization algorithm, stimulation patterns that produced maximum forces from each subject were identified.
 functional electrical stimulation
 fatigue
 human quadriceps femoris muscle
 doublets
 variablefrequency trains
electrical stimulation can be used to evoke skeletal muscle contractions that enable patients with central nervous system damage to perform functional movements such as ambulation (2, 19, 21). The practical use of functional electrical stimulation (FES) is limited by muscle fatigue, which develops far more rapidly during this exogenous activation of skeletal muscle than during activation by the central nervous system (25). The frequency of the stimulation train affects the forces produced (1) and the rate of muscle fatigue (3). Traditionally, during FES, skeletal muscles are activated with constantfrequency trains (CFTs), where the pulses within each train are separated by regular interpulse intervals (IPIs; Fig.1). More recently, however, it has also been shown that changing the pattern of the pulses within a brief train can markedly affect force production and fatigue (4, 7). Therefore, it is important to identify the stimulation pattern that maximizes force.
One stimulation pattern that has been shown to produce more force from a fatigued muscle than any CFT is a variablefrequency train (VFT) that begins with a doublet (i.e., 2 pulses separated by a brief IPI; Fig. 1) (6). However, the specific pattern of activation that maximizes the force varies from person to person (20) and, even for the same individual, depends on the physiological conditions of the muscle, such as the level of fatigue (6) or length of the muscle (24). Numerous tests, therefore, would be required to find the optimal patterns for each patient experimentally. One way to assist the search for the optimal pattern is to use mathematical models that can predict forces under a range of physiological conditions. Such models could markedly reduce the number of tests needed to identify the optimal pattern for each individual.
Although many models have been developed to explore different aspects of muscle contraction (911, 17, 28, 29, 33), most are complex (8,10, 28) and few have the ability to predict force responses to a wide range of stimulation patterns (13, 32). Our laboratory has recently developed a twostep model that predicts force responses to a wide range of stimulation patterns under various physiological conditions (13). This model is governed by two differential equations with four free parameters. It successfully predicted the forces in response to CFTs and VFTs with various IPIs for human quadriceps femoris muscles. However, a major limitation of this model is that separate sets of parameter values were needed to predict the forces produced in response to CFTs and VFTs when the muscles were fatigued. This resulted in a nonunique representation of each muscle by the model. The model, therefore, cannot be used to identify the pattern that produces the greatest forces.
The present study first outlines the modification of our twostep model to characterize each muscle with one set of parameter values regardless of stimulation pattern. Next, this new model is used to identify the activation pattern that produces the greatest forces. This new model suggested an activation pattern that we have not previously tested as the optimal for fatigued muscles. Finally, experimental results are presented that verify the model's predictions.
METHODS
Development of R Model
Modification of the twostep model.
Previously, we developed a twostep mathematical model for predicting isometric muscle forces that contained two differential equations and four free parameters (13). Briefly, the two differential equations are
As one of the simplest mathematical representations of muscle force, our previous twostep model overlooked many physiological details, such as the nonlinear summation of Ca^{2+} transients in single muscle fibers stimulated with two closely spaced pulses (15). Because VFTs contain doublets (closely spaced pulses), this could explain why the twostep model needed separate sets of parameter values to predict CFTs and VFTs. Similarly, we noted that our previous model did not accurately predict responses to highfrequency CFTs, which contain closely spaced pulses.
Duchateau and Hainaut (16) investigated the force summation from human adductor pollicis muscles triggered by paired stimuli at different IPIs ranging from 5 to 200 ms. They showed that the forces produced by the twopulse trains were greater than the sum of two individual twitches and that this force enhancement from the second pulse was highest when the IPI was 5 ms and declined exponentially with increases of the IPI. The enhanced force of the paired stimuli was suggested to be due to the intensified release of Ca^{2+} by the second pulse (15). In our new model, a factor R
_{i} is introduced to account for this nonlinear summation. We, therefore, modified the model by adding a factor modifying Eq. 1
Parameter identification.
Data from a previous study of the human quadriceps femoris muscles (13) were used to develop and test our new R model. Pilot studies showed that a fixed value of 20 ms for τ_{c} was sufficient for human quadriceps muscles under different physiological conditions. The remaining four free parameters must be identified for each individual under each physiological condition. We believe that trains or train combinations that exhibited tetanic forces, as well as a clear rise and fall in force for each pulse, would give the model the most information. Ten trains [CFTs with IPIs of 20 (CFT20), 30 (CFT30), 40 (CFT40), 70 (CFT70), and 100 ms (CFT100) and VFTs with IPIs of 10 (VFT10), 70 (VFT70), 80 (VFT80), 90 (VFT90), and 100 ms (VFT100)] were tested to determine which train or train combinations would be best to parameterize the model. We found that two VFTs, one with an initial 5ms IPI followed by constant IPIs of 10 ms (VFT10) and the other with an initial 5ms IPI followed by constant IPIs of 100 ms (VFT100), were the best train combinations to parameterize the model for nonfatigued muscles among the trains or train combinations tested. The combination of CFT100 and VFT100 was found to be the best for fatigued muscles. In our previous work (13,32), parameter identification was performed using SAAM II, a compartmentmodeling software package developed by the University of Washington. In this work we developed our own technique. First, τ_{1} was identified as described previously (13). Briefly,Eq. 2
became dF/dt = −F/τ_{1}by setting C_{N} = 0. By taking the force decay at the end of the forcetime response and performing a linear regression fitting of lnF vs. t, τ_{1} was obtained from the slope of the fit. The remaining three parameters were identified by employing the following objective function
We next tested the model's ability to predict the sixpulse stimulation pattern that produced the greatest forcetime integral. CFSQP was also used to maximize the forcetime integral by varying the IPIs separating each pulse within the sixpulse train. The forcetime integral was approximated by
Testing of the R Model
Twelve subjects were recruited for this part of study. Each subject signed an informed consent form approved by the University of Delaware Human Subject Review Board and participated in one testing session.
The experimental setup was similar to that described previously (5). Briefly, subjects were seated on a force dynamometer with their hips flexed to ∼75°, and their knees were 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 stimulus intensity was set to elicit a force equal to 20% of the subject's maximum voluntary isometric contraction when stimulated with a sixpulse, 100 pulse/s train. Force data from the force dynamometer's force transducer were entered into a personal computer via an analogtodigital converter board and digitized at a rate of 200 Hz.
Stimulation protocols.
Isometric muscle performance was tested using CFTs, VFTs, and DFTs when the muscles were fresh and fatigued (Fig. 1). All trains had six pulses (5 IPIs), because studies have shown that most functional movement employed a brief activation pattern (18). To cover a wide range of frequencies, the nine CFTs had all IPIs equal to 10, 30, 50, 70, 90, 100, 110, 130, or 150 ms. These trains are referred to as CFT10, CFT30, CFT50, CFT70, CFT90, CFT100, CFT110, CFT130, and CFT150, respectively. Comparable VFTs had an initial IPI of 5 ms, and the remaining pulses were equally spaced by 10, 30, 50, 70, 90, 100, 110, 130, or 150 ms. These trains are referred to as VFT10, VFT30, VFT50, VFT70, VFT90, VFT100, VFT110, VFT130, and VFT150, respectively. Six DFTs had the three 5ms doublets separated by 30, 50, 70, 90, 110, and 130 ms (i.e., the 2nd and 4th IPIs were long). These trains are referred to as DFT30, DFT50, DFT70, DFT90, DFT110, and DFT130, respectively.
The nonfatigue protocol consisted of the 24 testing trains presented in random order, which were then repeated in reverse order, for a total of 48 trains. Trains were delivered once every 10 s to avoid muscle fatigue. The forces produced in response to the second occurrence of each train were compared with the forces produced in response to the first occurrence to determine whether the muscles became fatigued during the nonfatigue protocol. None of the subjects showed a substantial (>5%) decline in the force of the second occurrence. This nonfatigue testing lasted ∼10 min. The subjects were then allowed 5 min of rest before the fatigue test.
To produce fatigue, a 13pulse CFT with IPIs of 25 ms (CFT25) was delivered to the muscle once every 1 s for a total of 100 contractions. Testing commenced the next second after the fatigueproducing stimulation. Trains continued to be delivered in 1s intervals. The 24 testing trains were delivered with each train preceded by two fatigueproducing stimulations (i.e., two 13pulse CFT25). This procedure was used to control for prior activation history of the muscle and to maintain a consistent level of fatigue throughout the fatigue testing (Fig. 2). Thus the fatigue protocol consisted of the 100 fatigueproducing trains, 24 testing trains, and 48 intervening fatigueproducing trains for a total of 172 trains and lasted ∼3 min.
Data management.
For the nonfatigued data, the two occurrences of each of the 24 stimulation patterns were averaged to minimize the effects of previous activation history on the muscle's responses to each train type.
TESTING THE MODEL'S ABILITY TO PREDICT FORCES.
The force responses to all CFTs, VFTs, and DFTs were used to test the model. For each subject, the model was first parameterized, as described above using G _{1} (Eq. 6 ) and CFSQP employing the VFT10VFT100 combination for the nonfatigued muscles and CFT100VFT100 combination for the fatigued muscles. The parameterized model was then used to predict the force responses to all other stimulation patterns. Two tests were used to examine the model's ability in predicting the force responses to a wide range of stimulation patterns. First, to compare the shape of the forcetime response to each stimulation pattern, a correlation coefficient was calculated by comparing the measured and predicted forces at 5ms intervals. Any phase shift or disagreement in shape between the predicted and measured data lowers the correlation coefficient. Second, for each subject, the forcetime integral was calculated for each train tested. The measured forcetime integrals were then plotted against the predicted forcetime integrals for each condition (i.e., nonfatigued and fatigued). Linear regression trend lines, with intercepts of zero, were used to compare force response across all the stimulation patterns and determine how well the model predicts the forcetime integrals. A perfectly accurate model would have a slope and R ^{2}of 1 (the identity line).
TESTING THE MODEL'S ABILITY TO PREDICT OPTIMAL PATTERNS.
Then, with the model parameterized for each subject under each condition, the optimizer (Eq. 7 ) was used to predict the activation pattern that produces the greatest forcetime integral. The upper limit of each IPI was set at 200 ms, and the lower limit was set at 5 ms. Because the program varies each IPI in a continuous manner and can test every possible combination of IPIs while searching for the optimal pattern, we first categorized the predicted optimal pattern to allow comparisons of the predicted results to our experimental data. Predicted optimal trains that had all five IPIs of similar value (standard deviation <10% of the mean) were called CFTs. Trains that had a short IPI (≤10 ms) in the beginning of the train and the other four IPIs of similar value were labeled VFTs. Trains that had three short IPIs separated by two long IPIs (>10 ms) were labeled DFTs. All short IPIs of the predicted optimal trains were equal to 5 ms except one (7 ms), and all predicted optimal trains fell into one of these three categories. Next, the long IPIs within each train (i.e., all 5 IPIs for CFT, the last 4 IPIs for VFT, and the 2 long IPIs for DFT) were averaged, and this average was used to designate the optimal IPI. The optimal IPI, along with the train type, was assigned as the optimal train pattern. The model's ability in predicting the activation pattern that produced the greatest forcetime integral was examined by comparing the train type and the IPI of the predicted optimal pattern with those of the experimental optimals.
RESULTS
The model was tested with data collected from nonfatigued (n = 12) and fatigued (n = 12) human quadriceps femoris muscles. One subject's fatigued data were not included into group analysis, because her fatigued forces were too low, resulting in an unacceptably high signaltonoise ratio. Subjective evaluations of the forcetime responses showed that the model was able to predict the force responses of the muscle to a wide range of stimulation patterns in the nonfatigued and fatigued states (Figs. 3and 4). The correlation coefficients between the measured and predicted forcetime responses, averaged across subjects, generally showed that the model more accurately predicted the force profiles when the muscles were nonfatigued (median = 0.94, SD = 0.06) than when they were fatigued (median = 0.88, SD = 0.05; Table 1). In general, for all three train types when the muscles were nonfatigued, the higher the average frequency (i.e., the shorter the IPI), the greater the correlation coefficient. When fatigued, CFTs and VFTs with 70 to 100ms IPIs showed the highest correlation coefficients.
For nonfatigued and fatigued conditions, the bestfit trend lines for modeled forcetime integrals vs. experimental data had slopes close to 1. The model accounted for ∼87 and 75% of the variance of the forcetime integrals of the experimental data for nonfatigued and fatigued muscles, respectively (Fig. 5).
The activation pattern that produced the greatest forcetime integral was predicted for each subject under each condition (i.e., nonfatigued and fatigued) and compared with the experimental data (Table2). For the experimental data, because of the physiological variability in the force responses of a muscle tested repeatedly with the same stimulation pattern, trains were considered to be producing the maximum forcetime integral from the muscle if they produced forces that were within 5% of the greatest force produced by any of the trains. For 20 of 23 comparisons (i.e., data from 12 nonfatigued and 11 fatigued subjects), the model and the experimental measurements identified the same train type as producing the maximum forcetime integral (Table 2). In addition, for 18 of the 20 matches described above, the model predicted the optimal IPI within 10 ms of the experimentally observed optimum (Fig.6).
Similar to our twostep model, fatigue produced a significant decrease in the parameter value A (4.8 ± 1.5 and 1.1 ± 0.5 for nonfatigued and fatigued muscles, respectively) but significant increases in τ_{1} (47.3 ± 7.6 and 89.2 ± 18.6 ms for nonfatigued and fatigued muscles, respectively) and τ_{2}(703.4 ± 383.6 and 1,564.5 ± 646.5 ms for nonfatigued and fatigued muscles, respectively). The new parameter,R _{0}, showed a significant increase when muscles were fatigued (0.9 ± 0.4 and 7.2 ± 3.7 for nonfatigued and fatigued muscles, respectively).
DISCUSSION
This new R model predicted accurately the force response to a wide range of stimulation patterns when muscles were fresh or fatigued. Subjective and objective evaluations of the force response showed that the predicted force responses closely matched the forces that were measured (Figs. 35, Table 1). In addition, the R model successfully identified the optimal pattern for each individual under different physiological conditions (Table 2, Fig. 6).
The present model was refined from our previous twostep model (13). Its structural simplicity and predictive ability remained intact, but the effect of shortterm activation history on the force production (R _{i} in Eq. 3 ) was taken into consideration. Previous study in our laboratory has shown that for nonfatigued and highly potentiated human quadriceps femoris muscles the force produced by the doublet in the VFT was slightly less than that produced by the first two pulses in the CFT, whereas when muscles were fatigued the doublet produced more forces than the first two pulses in the CFT (6). Thus values of R _{0} close to 1 for nonfatigued muscles and much greater than 1 for fatigued muscles were expected. The differences in the values of R _{0}between the nonfatigued and fatigued conditions are also consistent with the fact that, for most subjects, CFTs were optimal for the nonfatigued condition and DFTs or VFTs were optimal for the fatigued condition (Table 2). The major improvement over the twostep model is that the R model could characterize each individual muscle with only one set of parameter values regardless of stimulation pattern, allowing the model to predict the optimal patterns for each person under each condition. In addition, the present R model predicted more accurately the responses to highfrequency CFTs than our previous model.
An interesting finding was the identification of the DFT as the optimal pattern for the activation of the majority of subjects when muscles were fatigued. As an example, for one subject's fatigued muscle, the model and experimental data suggested that DFT130 produced the greatest forcetime integral among all trains tested (Fig. 6 B). The firing pattern of DFTs (i.e., repetitive doublets) has been observed in motoneuron discharge patterns of human trapezius muscles during voluntary contractions (22). In our study the minimum value allowed for the doublets was 5 ms, because pilot testing showed that trains with doublet values <5 ms produced equivalent or less force than trains with doublet values of 5 ms. In addition, Karu and colleagues (20) showed that a 5ms interval produced the greatest forcetime integral for twopulse trains. Maxwell and colleagues (26) examined the isometric force outputs from human quadriceps muscles stimulated with doublets varying from 2.5 to 20 ms. All DFTs tested produced greater forcetime integrals and peak forces than a CFT50, and the best DFTs had a doublet of 5 ms.
In this study we also explored the optimal interdoublet interval for DFTs that produced the greatest forcetime integral. Our experimental data showed that the DFT110, which had a train duration similar to the CFT50 and VFT50 (235, 250, and 205 ms, respectively), was the most common optimal activation pattern for the fatigued muscles (Table 2). Interestingly, Kudina and Alexeeva (22) showed that repetitive doublets fired most frequently at an interdoublet interval of 110 ms for human trapezius muscles during voluntary contractions. In contrast, Karu and colleagues (20) suggested 50–80 ms as the interdoublet interval to produce the greatest forcetime integral per pulse for human quadriceps muscles. The DFTs in their study had a fixed train duration (3 s) but a different number of pulses. In contrast, our present study used only sixpulse DFTs. Different optimal patterns may have emerged if we used longer trains or controlled the train duration rather than the number of pulses.
Two mechanisms have been proposed for the enhancement in force seen with doublets: increased Ca^{2+} transient (15) and increased stiffness (30) by the second pulse in the doublet. Impaired excitationcontraction coupling (31) and decreased muscle stiffness (12) with fatigue could make the positive effects of a doublet on force production more effective. In addition, unlike VFTs that have a doublet only in the beginning of the train, DFTs have doublets throughout the train. The positive effect of doublets would be used to advantage with each doublet.
Conclusion
The present model predicts isometric forces for human quadriceps femoris muscles in response to a wide range of stimulation patterns by use of brief trains of electrical pulses and under a variety of physiological conditions. Interestingly, this model predicted that DFTs produced greater forcetime integrals than traditionally used patterns under fatigued condition. The predictions of the model were verified with experimental data, demonstrating the predictive power and usefulness of this mathematical muscle model.
Acknowledgments
This study was supported by National Institute of Child Health and Human Development Grant HD36797.
Footnotes

Address for reprint requests and other correspondence: S. A. BinderMacleod, 315 McKinly Laboratory, University of Delaware, Newark, DE 19716 (Email: sbinder{at}udel.edu).

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. §1734 solely to indicate this fact.
 Copyright © 2000 the American Physiological Society