During FES, stimulation is delivered in the form of groups of pulses called trains. At any particular intensity of stimulation, both the stimulation frequency and pattern can be varied to control muscle force. Stimulation frequency can be varied by changing the duration of the inter-pulse intervals within a stimulation train. Stimulation trains that maintain a constant inter-pulse interval throughout a train are termed constant-frequency trains (CFTs). In contrast, trains with varying inter-pulse intervals within a train are called variable-frequency trains (VFTs) [16-18]. The most common type of VFTs that have been studied consist of two closely spaced pulses with 5 to 10-ms inter-pulse interval (doublet) at the onset of a CFT [16] (Figure ​(Figure1).1). Recently, trains consisting of regularly spaced doublets throughout the train, termed doublet-frequency trains (DFTs) have also been tested [16] (Figure ​(Figure1).1). VFTs and DFTs have been shown to augment muscle performance compared to CFTs of comparable frequencies, especially in fatigued muscles [16,19]. However, most commercial FES stimulators only deliver CFTs.

The generation of a sufficient isometric force level for a task is a prerequisite for effective performance of an FES-elicited task. For example, to manage foot drop using FES, the electrical stimulation parameters should elicit sufficient dorsiflexor muscle force to achieve ground clearance for numerous steps. However, the frequency or pattern of the stimulation train that generates the targeted performance may vary with the task, across individuals [18], between able-bodied and paralyzed muscles [20], and with the physiological condition of the muscle, such as fatigue or muscle length [21]. Thus, numerous measurements would be needed to identify the frequency and pattern that can generate the targeted forces during FES. Mathematical models that can predict the non-linear and time-varying relationships for each subject between stimulation parameters and electrically-elicited muscle forces can help reduce the number of testing sessions. When used in conjunction with a closed-loop controller, predictive mathematical models can enable FES stimulators to deliver customized, task-specific, and subject-specific stimulation patterns while continuously adapting these patterns to the changing needs of the patient [14,22].

Our laboratory has successfully developed a Hill-based [23] phenomenological mathematical model system that predicts muscle forces in response to stimulation trains of different patterns and a range of frequencies in able-bodied subjects [24,25] and individuals with spinal cord injury [26]. A recent comparative study [27] of muscle models that can be used in FES showed that our model [28] predicted electrically-elicited forces of the soleus muscles of individuals with chronic spinal cord injury as accurately as a 2^nd order nonlinear model [29] and with greater accuracy than a simple linear model. Another recent study [30] comparing 7 different muscle models showed our model [28], along with the Bobet-Stein model [29] provided the best fits for ankle dorsiflexor muscle forces over a range of joint angles in able-bodied individuals. However, the model has only been tested on able-bodied subjects and individuals with spinal cord injury. In addition, for our model to be successfully incorporated in a versatile FES-controller, it must predict force responses of a variety of lower extremity muscles in different patient populations. Therefore, our purpose was to test our model on the quadriceps femoris and ankle dorsiflexor and plantar-flexor muscles of individuals with hemiparesis following stroke. The three muscles tested in our study play an important role during functional activities such as ambulation [31,32] and are commonly impaired in individuals with post-stroke hemiparesis [33-37].

Typical measured and predicted force-time responses for the ankle dorsiflexor, and plantar-flexor muscles of a representative subject have been shown in Figures ​Figures33 and ​and4.4. Overall, the averaged FTI-frequency and PK-frequency relationships for CFTs, VFTs, and DFTs for the measured and the predicted force-time data matched well and there was consistency between the measured and predicted frequencies that generated the maximal FTI and PK for each of the muscles (See Figures ​Figures55 and ​and6).6). Interestingly, the model parameter values showed a high degree of variability across subjects and across the 3 muscles tested, with coefficients of variation ranging from 38% to 151% (Table ​(Table22).

The r² values comparing the shapes of the predicted versus measured force-time responses showed high levels of correlation between the predicted and measured forces (Figure ​(Figure7).7). For the quadriceps muscles, the r²values comparing the shapes of predicted and measured force-time responses were above 0.80 for all CFTs, VFTs, and DFTs (Figure ​(Figure7).7). For the dorsiflexor muscles, r² values were above 0.80 for all frequencies and patterns except the 10-Hz CFTs and 12.5-Hz DFTs (Figure ​(Figure7).7). For the plantar-flexor muscles, r² values were above 0.80 for all frequencies and patterns except the 12.5-Hz DFTs (Figure ​(Figure77).

ICCs comparing the measured versus predicted FTI and PK across all frequencies showed ICC values above 0.82 for the quadriceps, above 0.92 for the dorsiflexor muscles, and above 0.96 for the plantar-flexor muscles. In addition, scatter plots of predicted versus measured FTIs and PKs were plotted and the slopes of the trendlines with intercept set at zero were calculated. A perfect model would have ICC values and trendline slopes equal to one. In the current study, the trendline slopes for the 3 muscle groups tested ranged from 0.86 to 1.07 (Figure ​(Figure88).

The model error was within or below the +95% confidence interval of the physiological error for 91% of the comparisons between measured and predicted forces for the quadriceps, 94% of the comparisons for the dorsiflexor muscles, and 88% of the comparisons for plantar-flexor muscles (See Figure ​Figure9).9). The patterns for which the model errors was above the +95% confidence interval of the physiological error were: 25-Hz CFT PK, 20-Hz VFT PK, and 12.5-Hz DFT PK for the quadriceps; 10-Hz CFT PK and 10-Hz CFT FTI for the dorsiflexor muscles; 20-Hz CFT PK, 20-Hz VFT PK, and 12.5-Hz DFT PK and 12.5-Hz DFT FTI for plantar-flexor muscles (See Figure ​Figure99 for PK data).

Our laboratory previously developed a mathematical model system that successfully predicted responses to electrical stimulation during both isometric [26,39] and nonisometric contractions [44,45]. The present work was the first to test the mathematical force model on the muscles of individuals with hemiparesis following stroke. Unlike the 50-Hz CFT and 12.5-Hz VFT that were used to determine the parameter values in previous studies [26,39], in the present study we fit the measured forces in response to a 50-Hz CFT and 20-Hz DFT to determine the parameter values for each subject. Previously, we have found that fixing parameter τ_C at 20 ms for able-bodied individuals [39] and at 0.22 times each subject's half-relaxation time for individuals with spinal cord injury [26] enabled accurate predictions of isometric forces. Also, previously, parameter R₀ was determined using the relationship R₀ = K_m + 1.04 in healthy subjects [39], or set free in subjects with spinal cord injuries [26]. In this study, however, extensive preliminary testing showed that fixing τ_C at 11 ms and R₀ at 5 produced satisfactory predictions for all 3 muscles. As a result of fixing both τ_C and R₀ in the present study, the time course of C_N was identical across muscles and subjects, and only varied with the stimulation frequency and pattern (See equation (1)). Thus, parameter K_m was the only factor that determined the effect of C_N on force generation for different muscles and different subjects, and equation (2) played the primary role in prediction of muscle forces. The current version of our model, therefore, has the advantage of having only 4 free parameters, only equation (2) primarily governing force predictions, and the ability to predict muscle forces for multiple muscles of individuals with hemiparesis following stroke.

Muscles of individuals with stroke show changes in histochemical, morphometric, and structural properties compared to able-bodied individuals, with most studies reporting an increased percentage of type I fibers [46-48]. In addition, in the present study, reflex responses were often observed, especially for forces measured from the quadriceps and plantar-flexor muscle groups (See Figures ​Figures33 and ​and55 for examples). In light of marked differences in muscle physiology following stroke [46-48] and the presence of reflex activation, a modified interpretation of the model parameter values was needed in the present study. In our previous works, all measured muscle force responses were directly in response to electrical stimulation. However, in the present study, the measured forces were often a result of the combination of synchronous activation of motor units by the electrical stimulation and the asynchronous reflex activation. Reflex activation was manifested by the lack of return of force to baseline after the stimulation was turned off, and by the lack of one-to-one correspondence between the individual pulses in the stimulation trains and the shape of the force-time responses. Furthermore, because we identified the model parameter values by fitting the measured forces, the model parameter values reflected both the effects of forces in response to electrical stimulation and any reflex responses produced. In the present study, therefore, we slightly modified the physiological interpretation of 4 of our model parameters - C_N,, K_m, τ_1, and τ₂. Previously, C_N was defined as a unitless representation of the Ca²⁺-troponin complex; K_m was defined as the sensitivity of strongly bound cross-bridges to C_N ; τ₁ and τ₂were defined as time constants of force decline in the absence and presence of strongly bound cross-bridges, respectively [39]. In the present study, we reinterpreted C_N as the rate-limiting step before the myofilaments mechanically slide across each other and generate force; K_m as the sensitivity of force development to C_N ; τ₁ and τ₂ as time constants modeling the force decay, both electrically induced and reflexive, due to the visco-elastic components of muscle in the absence and presence of stimulation, respectively. Thus, the parameter definitions in the present study were more generic and applicable to force responses of whole muscles to electrical stimulation, both in the presence and absence of reflex activation.

Unlike our previous modeling studies [26,39], in the present study we found a high degree of variability in the model parameter values within each muscle tested across the individuals with stroke, with coefficients of variation ranging between 38% and 151% (See Table ​Table2).2). This variability in parameter values may be in part due to differences across subjects in the level of reflex activation of their muscles. Our model works under the assumption of synchronous activation of the whole muscle in direct response to electrical stimulation. However, in muscles of individuals with stroke, this assumption was violated due to the presence of reflex activation. Nevertheless, in spite of marked differences in properties of muscles of hemiparetic individuals and the presence of reflex activation, our model was able to predict isometric forces for 80% of the muscles studied.

A recent study comparing 3 different muscle models found that the 2^nd order nonlinear model developed by Bobet and Stein (1998) and our model [38], both consisting of 6 parameters, showed better predictions of muscle forces than a simple linear model, especially for higher frequency or variable frequency trains [27]. Furthermore, our model produced the least percentage error between measured and predicted among the three types of models compared [27]. In addition to the accuracy of our model demonstrated by two recent comparative studies of muscle models [27,30], the current version of the model has the added advantage of only 4 free parameters. Interestingly, in a recent sensitivity analyses of 3 muscle models, Frey Law and Shields [49] suggested that due to the influences of parameter τ_C on muscle force properties predicted by our model, we should perhaps keep the value of parameter τ_C free. However, both in our previous work [39] and in the present study, we found that our model accurately predicted muscle forces despite fixing τ_C. Furthermore, during the preliminary testing of the model in the present study, we found that for all 3 muscles, the model predicted FTIs and PKs with greater accuracy when τ_C was fixedat 11 and R₀ was fixed at 5, compared to either when both τ_C and R₀ were free or when τ_C was fixed at 20 and R₀ was fixed at 2 (values employed for able-bodied subjects) [25,39]. Law and Shields suggested that interactions between parameters are likely to exist [49], hence we hypothesize that the parameters τ₁, τ₂, and Km were able to compensate for the fixed τ_C . Because the fewest parameters are desirable for an ideal feedforward model in FES [22], and because preliminary testing showed that fixing both τ_C and R₀ did not compromise the predictive ability of the model, we decided to reduce the number of free parameters in our model by keeping the activation dynamics (See equation 1) constant for a given frequency and pattern across subjects and muscle groups.

A practical FES system needs both a feedforward model, which designs subject-specific and task-specific stimulation patterns, and a feedback controller, which corrects errors by informing the feedforward model when changes in stimulation patterns are needed [50]. When used in conjunction with a closed-loop controller, mathematical models can allow FES stimulators to deliver patient-specific and task-specific stimulation patterns that can adapt to the actual needs of the patient in real-time [14,26]. The use of customized stimulation patterns can reduce the energy expenditure and improve the speed at which functional tasks are performed during FES [14,22,26]. A model with the fewest parameters that can accurately predict the PK and FTI in response to a wide range of frequencies and patterns is desirable for a feedforward model in FES-systems [22,50]. The present model can accurately predict forces of 3 different muscles, which are important muscles for ambulation [31,32] and are commonly impaired in individuals with post-stroke hemiparesis [33-37]. Thus, because the present model can predict accurately and consists of only 4 free parameters, it has potential for use as a feedforward model in FES-systems for individuals with hemiparesis following stroke.