Development of a mathematical model that predicts optimal muscle activation patterns by using brief trains

Jun Ding, Anthony S. Wexler, Stuart A. Binder-Macleod


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 high-frequency burst of pulses (i.e., an initial doublet or triplet) followed by a low, constant-frequency 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
  • variable-frequency 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 constant-frequency 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.

Fig. 1.

Schematic representation of 3 stimulation train types. Bottom train (CFT30), constant-frequency train with all interpulse intervals (IPIs) equal to 30 ms; middle train (VFT30), variable-frequency train with an initial doublet of 5 ms and remaining pulses equally spaced by 30 ms; top train (DFT30), DFT with 5-ms doublets separated by intervals of 30 ms. Each train's name is based on duration of its longer IPI (e.g., VFT50 contains an initial doublet of 5 ms followed by pulses equally spaced by 50 ms).

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; 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 (9-11, 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 two-step 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 two-step 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.


Development of R Model

Modification of the two-step model.

Previously, we developed a two-step mathematical model for predicting isometric muscle forces that contained two differential equations and four free parameters (13). Briefly, the two differential equations aredCNdt=1τci=1nexpttiτcCNτc Equation 1anddFdt=ACN1+CNFτ1+τ2CN1+CN Equation 2representing the chemical and physical events that result in force development.Equation 1 characterizes the rate-limiting step to the formation of the Ca2+-troponin complex (C N, unitless); i identifies each stimulation pulse within the train (e.g., 1st, 2nd, or 3rd pulse); τc (in ms) is the time constant controlling the rise and decay of C N; t (in ms) is the time since the first pulse of the train; ti (in ms) is the time of the ith stimulation pulse; nis the total number of stimulation pulses before t (its maximum value is 6 in this study). F (in N) in Eq. 2 represents the whole muscle force, the development of which is driven byC N/(1 + C N), the force producing cross bridges, and is scaled by parameter A(in N/ms). The second term in Eq. 2 accounts for the force decay over two time constants, τ1 (in ms) and τ2 (in ms); τ1 represents the time constant in the absence of actin-myosin cross bridges, and τ2 represents the extra friction between actin and myosin resulting from the presence of actin-myosin cross bridges. In Eqs.1 and 2 , C N and F are the two dependent variables and τc, A, τ1, and τ2 are free parameters.

As one of the simplest mathematical representations of muscle force, our previous two-step model overlooked many physiological details, such as the nonlinear summation of Ca2+ transients in single muscle fibers stimulated with two closely spaced pulses (15). Because VFTs contain doublets (closely spaced pulses), this could explain why the two-step 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 high-frequency 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 two-pulse 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 Ca2+ 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 dCNdt=1τci=1nRiexpttiτcCNτc Equation 3whereR i is a scaling term that accounts for the differences in the degree of activation by each pulse relative to the first pulse of the train. The analytic solution of Eq. 3 isCN=i=1nRittiτcexpttiτc Equation 4whereRi=1+(R01)exptiti1τc Equation 5fori > 1 and R 1 = 1. The magnitude of the enhancement is characterized by R 0, and its duration is characterized by τc. This new model, termed the R model, has two differential equations (Eqs. 2 and 3 ) and is governed by five free, constant parametersR 0, τc, A, τ1, and τ2.

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 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), 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 compartment-modeling 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/τ1by setting CN = 0. By taking the force decay at the end of the force-time 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 functionG1(A,R0,τ2)=i[Fpred(ti;A,R0,τ2)Fobs(ti)]2 Equation 6where F pred, the forces predicted by Eqs. 2 and 3 at time t i, are a function of the four parameters A, R 0, τ1 ,and τ2. Fobs are the observed forces attime t i. G 1 was minimized using CFSQP (23), which numerically identifies the optimum values for the three parameters (A, R 0, τ2) by employing feasibility set techniques coupled to sequential quadratic programming.

We next tested the model's ability to predict the six-pulse stimulation pattern that produced the greatest force-time integral. CFSQP was also used to maximize the force-time integral by varying the IPIs separating each pulse within the six-pulse train. The force-time integral was approximated byG2(IPI¯)=i[Fpred(ti;IPI¯)+Fpred(ti+1;IPI¯)]2(ti+1ti) Equation 7 I̅P̅I̅is the set of five IPIs in the model stimulation and (t i + 1t i) is the time between the ith and the (i + 1)th pulse. The minimum duration allowed for each IPI was 5 ms, because pilot testing showed that reducing the IPI to <5 ms had no effect on force or produced a decline in force. The present model, however, is not capable of demonstrating this decline in force at very short IPIs. In addition, pilot testing of the model showed that for fatigued muscles a stimulation pattern that contained doublets throughout (i.e., not only in the beginning of the train) should produce the greatest forces for some muscles. This pattern, denoted as a DFT (Fig. 1), has rarely been tested previously (20, 26). The following experiments, therefore, were conducted to verify the predictions of the model by use of a wide range of six-pulse stimulation patterns.

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 six-pulse, 100 pulse/s train. Force data from the force dynamometer's force transducer were entered into a personal computer via an analog-to-digital 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 5-ms 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 13-pulse 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 fatigue-producing stimulation. Trains continued to be delivered in 1-s intervals. The 24 testing trains were delivered with each train preceded by two fatigue-producing stimulations (i.e., two 13-pulse 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 fatigue-producing trains, 24 testing trains, and 48 intervening fatigue-producing trains for a total of 172 trains and lasted ∼3 min.

Fig. 2.

Schematic representation of stimulation protocols. All testing trains were randomized.

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.


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 VFT10-VFT100 combination for the nonfatigued muscles and CFT100-VFT100 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 force-time response to each stimulation pattern, a correlation coefficient was calculated by comparing the measured and predicted forces at 5-ms intervals. Any phase shift or disagreement in shape between the predicted and measured data lowers the correlation coefficient. Second, for each subject, the force-time integral was calculated for each train tested. The measured force-time integrals were then plotted against the predicted force-time 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 force-time integrals. A perfectly accurate model would have a slope and R 2of 1 (the identity line).


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 force-time 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 force-time integral was examined by comparing the train type and the IPI of the predicted optimal pattern with those of the experimental optimals.


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 signal-to-noise ratio. Subjective evaluations of the force-time 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 force-time 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 100-ms IPIs showed the highest correlation coefficients.

Fig. 3.

Measured and predicted forces from a typical subject's nonfatigued quadriceps muscle. A: model was first parameterized by fitting to force responses to VFT10 and VFT100. exp, Experimental; mod, modeled. B–J: experimental and predicted force responses to other CFTs, VFTs, and DFTs.

Fig. 4.

Measured and predicted forces from a subject's fatigued quadriceps muscle (same subject presented in Fig. 3). A: model was first parameterized by fitting to force responses to CFT100 and VFT100. B–J: experimental and predicted force responses to other CFTs, VFTs, and DFTs.

View this table:
Table 1.

Averaged correlation coefficients

For nonfatigued and fatigued conditions, the best-fit trend lines for modeled force-time integrals vs. experimental data had slopes close to 1. The model accounted for ∼87 and 75% of the variance of the force-time integrals of the experimental data for nonfatigued and fatigued muscles, respectively (Fig. 5).

Fig. 5.

Comparisons of force-time integrals between measured and predicted forces of 12 nonfatigued (A) and 11 fatigued (B) human quadriceps muscles.

The activation pattern that produced the greatest force-time 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 force-time 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 force-time 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).

View this table:
Table 2.

Measured and predicted optimal activation pattern

Fig. 6.

Model's ability to predict optimal pattern for 1 typical subject under nonfatigued (A) and fatigued (B) conditions. ○, □, and ⋄, Experimental force-time integrals produced by CFTs, VFTs, and DFTs, respectively, plotted against long-duration IPI within each train. Arrows, optimal patterns predicted by model (i.e., CFT66 and DFT130, respectively).

Similar to our two-step 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).


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. 3-5, 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 two-step model (13). Its structural simplicity and predictive ability remained intact, but the effect of short-term 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 0between 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 two-step 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 high-frequency 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 force-time 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 5-ms interval produced the greatest force-time integral for two-pulse 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 force-time 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 force-time 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 force-time 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 six-pulse 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 Ca2+ transient (15) and increased stiffness (30) by the second pulse in the doublet. Impaired excitation-contraction 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.


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 force-time 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.


This study was supported by National Institute of Child Health and Human Development Grant HD-36797.


  • Address for reprint requests and other correspondence: S. A. Binder-Macleod, 315 McKinly Laboratory, University of Delaware, Newark, DE 19716 (E-mail: sbinder{at}

  • 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.


View Abstract