## Abstract

Chronic spinal cord injury (SCI) induces detrimental musculoskeletal adaptations that adversely affect health status, ranging from muscle paralysis and skin ulcerations to osteoporosis. SCI rehabilitative efforts may increasingly focus on preserving the integrity of paralyzed extremities to maximize health quality using electrical stimulation for isometric training and/or functional activities. Subject-specific mathematical muscle models could prove valuable for predicting the forces necessary to achieve therapeutic loading conditions in individuals with paralyzed limbs. Although numerous muscle models are available, three modeling approaches were chosen that can accommodate a variety of stimulation input patterns. To our knowledge, no direct comparisons between models using paralyzed muscle have been reported. The three models include *1*) a simple second-order linear model with three parameters and *2*) two six-parameter nonlinear models (a second-order nonlinear model and a Hill-derived nonlinear model). Soleus muscle forces from four individuals with complete, chronic SCI were used to optimize each model's parameters (using an increasing and decreasing frequency ramp) and to assess the models' predictive accuracies for constant and variable (doublet) stimulation trains at 5, 10, and 20 Hz in each individual. Despite the large differences in modeling approaches, the mean predicted force errors differed only moderately (8–15 % error; *P* = 0.0042), suggesting physiological force can be adequately represented by multiple mathematical constructs. The two nonlinear models predicted specific force characteristics better than the linear model in nearly all stimulation conditions, with minimal differences between the two nonlinear models. Either nonlinear mathematical model can provide reasonable force estimates; individual application needs may dictate the preferred modeling strategy.

- spinal cord injury
- electrical stimulation
- muscle mechanics

it is well recognized that chronic spinal cord injury (SCI) induces a variety of muscle adaptations including atrophy (8, 10, 39), faster contraction time (28, 45), faster relaxation time (28, 45), increased fatigability (9, 28, 41, 45), and myosin heavy chain transformation to faster isoforms (7, 10, 30, 39, 40, 47) relative to normal muscle. Furthermore, the health status of individuals with SCI may be compromised by the loss of volitional muscle activity, including skin ulcerations, cardiovascular disease, kidney and urinary tract infections, and fractures secondary to osteoporosis. Traditional rehabilitation techniques for individuals with SCI typically focus on optimizing function after injury. However, SCI rehabilitation may increasingly include methods to preserve the integrity of the paralyzed limbs as a means to improve overall health quality. One approach is to induce isometric muscle contractions using neuromuscular electrical stimulation with the goal of preventing adverse musculoskeletal adaptations resulting from disuse (26). Moreover, if a cure to SCI is discovered, it may only be accessible to those whose bodies have not undergone substantial musculoskeletal deterioration. Even today, secondary complications associated with musculoskeletal deterioration after SCI can prove to be life threatening (44). The use of subject-specific mathematical predictions of electrically stimulated muscle force may aid in the development of individualized rehabilitation methodology to optimally stress the musculoskeletal system after SCI.

In the past decade, work has increasingly focused on the use of mathematical muscle models to predict muscle forces (4, 13, 23) for isometric exercise training regimens and/or mobility using functional electrical stimulation. Numerous muscle models have been built on the original phenomenological 1938 Hill model (25, 33) as well as the structural, cross-bridge approach introduced by Huxley in 1957 (25, 34), ranging from straightforward linear models (2, 3) to a variety of complex nonlinear models (11, 13, 23, 31, 46, 48, 49). Models of individual force production have been validated in animal muscle (5, 6, 29, 36), human able-bodied muscle (13, 16–18, 21), and paralyzed muscle (14). These models use one or two stimulation trains to fit the underlying subject-specific models' parameters, thereby predicting force for a range of potential stimulation trains. The models allow the determination of an optimal or ideal stimulation pattern for an individual without exposure to unnecessary stimulation bouts. However, to our knowledge, no between-model comparisons using paralyzed muscle have been reported. These are important investigations that are needed to better understand which muscle models may be best suited for an individual application. Moreover, the complexity required to adequately model human paralyzed muscle (e.g., nonlinear vs. linear) is not well understood.

Three modeling approaches that are able to accommodate variable frequency stimulation patterns stand out as promising for the development of therapeutic models to stress the musculoskeletal system after SCI. Accordingly, the purpose of this novel investigation was to systematically evaluate the accuracy of three muscle models for predicting human chronically paralyzed muscle forces. We hypothesized that the more complex nonlinear models (the Hill-based and second-order nonlinear models) would provide the most accurate predictions of overall force as well as specific force properties [peak force (PF), work, time to peak, half-relaxation, fusion level, and “catch-like property” of muscle]. Hill-based models are frequently believed to have physiologically related model parameters, suggesting that a Hill-based model could provide the best estimates of specific force properties (PF, time to peak, etc.).

## METHODS

Three currently available muscle models that are capable of accommodating discrete, variable frequency inputs were chosen. The first and simplest model is based on a traditional engineering approach to modeling dynamic systems, a second-order linear model (for an overview of linear systems theory, please see Ref. 12). The second and third models are contemporary nonlinear mathematical models, each with six constant parameters.

#### Linear model.

A single differential equation (1) is used to represent force as f(*t*) with a train of input stimulation pulses [*r*(*t*)] represented by a series of pulses, or dirac delta functions [δ(*t*)], denoted by r(*t*) = Σ δ(*t*_{i} − *t*_{0}) for i = 1 to *n*. The δ(*t*) function equals 1/d*t* at the time (*t*_{i}) of each pulse, for a total area equal to 1.0, and zero otherwise. (1) Essentially, the forces resulting from each consecutive pulse are summed linearly to produce the overall force train. The three constant parameters involved, β, ζ, and ω_{n}, correspond to the gain of the system (Ns), the damping ratio (dimensionless), and the undamped natural frequency (rad/s), respectively (12). The gain provides a means to increase the force time area for a given input, influencing the maximum force production (27). The damping ratio, ζ, controls the oscillatory nature and the “responsiveness” (3) of the modeled force output and has little physiological meaning in isolation. It is mathematically useful for controlling the stability of a modeled system, i.e., a stable system will have oscillations that will remain constant or dampen out over time. A combination of ω_{n} and ζ determines the overall system time constant (22). For modeling human muscle, values for ζ are typically ≤1.0 (critically damped) (1, 3, 5). The natural frequency, ω_{n}, of the linear system may be more directly related to the physiological speed properties of a muscle (3, 27). Faster muscle speed properties in isolated animal preparations (e.g., contraction and relaxation times) have corresponded to higher values for ω_{n} (6).

#### Second-order nonlinear model.

The second-order nonlinear model proposed by Bobet and Stein (5) is represented by *Eqs. 2–5*. In addition to two first-order differential equations [solutions g(*t*) and F(*t*) shown as *Eqs. 2* and *4*], this model includes a saturation nonlinearity that limits force at higher levels [*Eq. 3*, where the intermediate output, *x*(*t*), is modulated by parameters *k* and *n*]. A variable rate constant parameter, b, defined by *Eq. 5*, varies with force and is modulated by parameters b_{0}, b_{1}, and B (5). When b_{1} is positive, this variable rate constant becomes slower with increasing force and faster when b_{1} is negative. Parameter B acts as a simple gain factor in this model (27). (2) (3) (4) (5) The final output from *Eq. 4*, F(*t*), is the time history of predicted force. These four equations use a total of six constant parameters, a, B, b_{0}, b_{1}, *n*, and *k*, to modulate the force output.

#### Hill-Huxley-type nonlinear model.

The third and most complex muscle model included in this study has been frequently described in the literature (13, 15–21, 37). This model was originally based on Hill's phenomenological model (48), but its evolution over the past several years has increased its resemblance to the distribution-moment model (50), a Huxley-based modeling approach. The original Hill elastic element terms are difficult to recognize in the model's current state. Consequently, it will be referred to as the Hill-Huxley-type nonlinear model for this discussion.

The most recent model version for nonfatigued muscle is comprised of two differential equations (*Eqs. 6* and *7*) and a third equation defining a variable parameter (*Eq. 8*) (17, 18). A total of six constant parameters, A, τ_{1}, τ_{c}, τ_{2,} R_{0}, and *k*_{m}, modulate this Hill-Huxley modeling approach. The first differential equation (*Eq. 6*) has been reported to represent muscle contraction calcium kinetics (both the release/reuptake of Ca^{2+} as well as troponin binding) (16, 18, 21). This equation is modulated by the variable parameter, R_{i} (*Eq. 8*), which decays as a function of each interpulse interval (*t*_{i} to *t*_{i-1}) from a maximum value of R_{0} for very rapid pulses to a minimum of 1 for widely spaced pulses. An analytical solution for the first differential equation (C_{N}) is shown in *Eq. 9* (16, 18, 21). However, no analytical solution to the differential force equation (*Eq. 7*) exists; thus iterative numerical estimation techniques must be used. This iterative process to estimate an approximate solution adds to the challenge of solving this model. Furthermore, many optimization techniques involve gradient-based (derivative) procedures, which is problematic with numerical estimation techniques. Although this model also has only six parameters, it takes substantially more time and iterative steps to solve than does the second-order nonlinear model, making it the most computationally complex of the three. (6) (7) (8) (9)

#### Experimental setup.

Four individuals with chronic complete SCI were recruited for this study (see Table 1 for subject characteristics). All subjects were otherwise healthy with no known medical conditions other than SCI. All experimental procedures were conducted according to protocols approved by the University of Iowa Institutional Review Board. Written, informed consent, as approved by the institutional review board, was obtained before inclusion in this study.

Subjects remained seated in their own wheelchairs with their left foot placed on a custom-made ankle force plate. This measurement system has been previously described in detail (41). Briefly, the foot was firmly secured to the force measurement device, with the knee and ankle joints positioned at ∼90°. Force was measured perpendicularly to the metatarsal heads (15-cm moment arm) using a Genisco AWU-250 force transducer (Genisco), amplified with a gain of 500 (200 for *subject 1*). Error range was determined to be 0.01–1.49% full scale (defined as 0–840 N), with maximum capacity of 1,112 N and overall mean absolute error of 4.7 N. Muscle electromyographic measurements of the soleus were obtained using preamplified, 0.8-cm diameter, bipolar silver-silver chloride surface electrodes with a fixed 2-cm interelectrode distance (Therapeutics Unlimited). Preamplification gain was 35 times, and further differential amplification was performed with a gain of 5,000.

A custom dual-pronged nerve electrode was used to stimulate the tibial nerve in the popliteal fossa, using computer-driven square waves generated by a custom-made constant-current stimulator using 250-μs pulsewidths (41). The placement of the electrode was optimized to produce the largest soleus M-wave peak-to-peak amplitude. Stimulation intensities were then increased to ∼1.5 times greater than required for a maximum M wave to ensure supramaximal stimulation throughout the protocol.

The force, electromyographic, and stimulation trigger signals were recorded simultaneously on magnetic (VCR) tape via a pulse code modulator (model 300, Vetter Digital), digitized at 1,000 Hz, and later analyzed using Matlab 6.0 software (Mathworks).

#### Stimulation protocol.

The stimulation protocol began with a series of 1-Hz twitches to determine the optimal electrode placement and ensure supramaximal levels of activation. Five warm-up contractions and seven-pulse trains (20-Hz trains, 330 ms on, 670 ms off) applied at 1 train/s were then given to ensure the soleus muscle was in a more potentiated state without inducing fatigue. Following the warm-up, a preprogrammed set of 10 stimulation trains was performed with 2-s rest intervals between each train for a total duration of ∼30 s. The set consisted of a doublet ramp of 15 pulses with interpulse intervals of 200, 6, 127, 100, 67, 67, 50, 50, 67, 67, 6, 94, 133, and 200 ms (1,234-ms duration; Fig. 1); three constant-frequency trains (CT) at 5, 10, and 20 Hz (8 pulses, 1,400-ms duration; 10 pulses, 900-ms duration; and 12 pulses, 550-ms duration, respectively); three single doublet trains (DT) at the base frequencies of 5, 10, and 20 Hz, with an additional doublet pulse 6 ms immediately after the first pulse; and three dual doublet trains (DDT) at the base frequencies of 5, 10, and 20 Hz, with two added pulses (doublets) to the CTs, one 6 ms after the first pulse and a second 6 ms after the middle pulse of each train (the 5th, 6th, and 7th pulse for 5-, 10-, and 20-Hz DDTs, respectively). See Fig. 1 for a schematic representation of the different force train patterns used.

These diverse input patterns represent a relatively wide range of possible stimulation approaches across the force-frequency curve of paralyzed soleus muscle (41–43, 45) by using a small number of trains. To minimize the effects of fatigue, a short overall duration of stimulation is necessary (∼30 s) in untrained paralyzed muscle. The mean on-off time ratio for the 10 stimulation trains was ∼1:2. A maximum stimulation frequency of 20 Hz, in concert with 167-Hz doublets, minimized the risk of injury to each subject, which has occurred with higher frequency stimulation (32).

#### Model parameterization.

All three models were parameterized using a single predetermined stimulation pattern, the doublet ramp (Fig. 1), modified from Bobet and Stein (5). Previously, the models have been parameterized using different types of muscle trains (5, 13, 16, 48) and a priori (5) vs. trial-and-error (13, 16, 17) methodologies. To appropriately compare results between models, one approach was necessary. The doublet ramp incorporates a range of frequencies (5, 7.5, 10, 15, and 20 Hz, with 6-ms doublets = 167 Hz) over a relatively short duration (1,234 ms) and provides approximately equal time at the low and high frequencies to minimize bias from any one frequency.

The model parameters were mathematically optimized using the Matlab optimization toolbox function, lsqnonlin, a least-squares approach using the Levenberg-Marquardt method, i.e., the lsqnonlin function minimizes the sum of the squares of the differences between model and experimental forces at each point in time. The Hill-Huxley-type nonlinear model required a fixed-step (vs. the more common variable step) Runge-Kutta algorithm for the optimization routine to be able to make estimates of the gradient and converge to an optimal solution. The variable-step numerical techniques frequently converged to local minima that varied with parameter starting values (e.g., not true optimal solutions). Although a global minimum cannot be mathematically guaranteed using most optimization techniques, further investigation of the possible solution space for each model supported that the minima found in this manner were good estimates of the global minima.

#### Model error estimates.

Using the optimal parameter values determined for each subject and model, force-train predictions for nine stimulation pattern/frequency combinations (CT, DT, and DDT trains at 5, 10, and 20 Hz) were produced. An overall error estimate, mean percent error (% error), was calculated as the mean absolute difference between the predicted and experimental force time profiles for each stimulation pattern/frequency (sampled at 1,000 Hz) normalized by the experimental PF.

The specific force characteristic errors determined for each model included percent PF (%PF); percent force-time integral (FTI); time to peak tension (TPT; defined as the time to reach 90% of PF from force initiation); half relaxation time (RT_{1/2}; defined as the time for force to decay from 90 to 50% of the final force peak); relative fusion index [RFI; a measure of tetany (fusion) defined as the mean difference between the last four force peaks (maxima) and their preceding force minima, normalized by the force maxima (100% indicates full fusion, 0% constitutes a series of twitches)]; and the “catch-like” property of muscle [% doublet difference (%DDiff); defined as the additional normalized force due to an extra doublet pulse]. The mean absolute force difference between the DDT minus DT trains and the DT minus CT trains, normalized by the maximum experimental doublet difference, was used to calculate the %DDiff. The %PF and %FTI were calculated for the CT and DDT trains; the TPT, RT_{1/2}, and RFI were determined using the CT trains (5 Hz only for TPT); and the %DDiff used the CT, DT, and DDT trains. Both absolute and directional errors were calculated for each of the specific force property error measures, with the exception of only absolute error for %DDiff. Absolute error provides a measure of the error magnitude only, whereas directional errors can indicate mean over- vs. under-predictions.

#### Repeatability.

Between-day testing resulted in good repeatability of the experimental force time profiles (*n* = 3). Mean errors in the doublet ramp (the parameterization train and most complex pattern) included 4.4% overall error, 7.6% PF, 2.0% force time integral, 3.3-ms time-to-peak tension, and 1.0-ms RT_{1/2}. Mathematical models should not be expected to predict forces more accurately than can be obtained experimentally.

#### Statistical methods.

The overall %error, %PF, and %FTI error measures were statistically analyzed using a three-way (3 × 2 × 3) univariate repeated-measures ANOVA for within-subject analyses of model (Hill-Huxley-type nonlinear, second-order nonlinear and linear), stimulation pattern (CT vs. DDT), and stimulation frequency (5, 10, and 20 Hz) for simple and interaction effects (SAS 9.0, SAS Institute, Cary, NC). The error magnitudes for RT_{1/2}, RFI, and %DDiff were analyzed using a two-way (3 × 3) repeated-measures ANOVA for model and frequency. TPT was evaluated using a one-way repeated-measure ANOVA (for model). Repeated-measures *P* values were adjusted using the Huynh-Feldt epsilon when necessary to correct for nonsphericity or heterogeneous variance and intrasubject correlation (35). Significance was set at α = 0.05. The Bonferroni correction for multiple comparisons (*n* = 3) was used for all follow-up tests to preserve the error level for each measure. Pilot analyses predicted 80% power would be achieved to detect an effect size (Cohen's d) of 2.0 for %error using three subjects in a repeated design.

## RESULTS

Representative examples of the experimental data, including stimulation pulses, soleus electromyographic signals, and the force output for the doublet ramp, 5-Hz single DT, 10-Hz DDT, and 20-Hz CT are shown in Fig. 2. The two nonlinear models were able to more closely fit (via optimization) the experimental doublet ramp pattern used to determine the parameter values than the linear model (see Fig. 3). Using the optimized parameter values, a representative example (*n* = 1) of several predicted force trains (CT and DDT) are displayed in Fig. 4 relative to the experimental force trains. Note the ability for all three models to represent paralyzed muscle forces reasonably well, particularly the two nonlinear models. Each subject's optimal model parameter values are provided in Table 2.

#### Overall %error.

These three muscle models were able to predict human paralyzed muscle force trains with mean errors ranging from 8.1 to 15.2% (Fig. 5) with the two nonlinear models more closely predicting the actual muscle forces than the linear model. The individual subject results (Fig. 5*A*) demonstrate similar between-model error trends (Hill-Huxley-type ≤ second order nonlinear < linear); the individual error magnitudes of the three models varied between subjects. Overall, %error resulted in significant main effects for model [*F*(2,6) = 28.9, Huynh-Feldt adjusted *P* value = 0.0042; Fig. 5*B*] as well as significant interactions for model by frequency [*F*(4,12) = 13.18, adjusted *P* value = 0.0127; Fig. 5*C*], model by stimulation pattern [*F*(2,6) = 15.47, adjusted *P* value = 0.0161; Fig. 5*D*], and model by frequency by stimulation [*F*(4,12) = 18.65, adjusted *P* value < 0.0001]. The between-model follow-up tests (Table 3) indicated the Hill-Huxley-type model produced the least overall %error and the linear model produced the greatest %error. However, the difference between the nonlinear models, although statistically significant, was small (<2%), and the differences between the linear and the nonlinear models were only 5–7% (Table 3). Post hoc power analyses revealed that this study was able to detect pairwise effect sizes (Cohen's d) of 1.23, 3.79, and 3.02 for the two nonlinear models, the Hill-Huxley-type vs. linear models, and the second-order nonlinear vs. linear models, respectively, despite the small sample size.

#### Specific force property errors.

Significant between-model differences resulted for both directional and absolute %PF and %FTI errors [*F*(2,6) values from 8.82 to 45.92, with Huynh-Feldt adjusted *P* values from 0.0404 to 0.0002; see Figs. 6 and 7]. Both nonlinear models were significantly better able to predict PF than the linear model (Table 3), but for FTI predictions only the second-order nonlinear model was significantly better than the linear model (Table 3). No differences between the two nonlinear models were observed for either PF or FTI. Significant two- and three-way interactions (model × frequency, model × stimulation pattern, and model × frequency × stimulation pattern) resulted for all directional and absolute %PF and %FTI errors, except the model × stimulation pattern interaction for the absolute PF [*F*(2,6) = 0.28, *P* = 0.7575]. See Figs. 6, *C* and *D*, and 7, *C* and *D*, to more clearly see how PF and FTI errors varied across frequency and with stimulation pattern.

Significant between-model differences were observed for the absolute and directional errors for the force time properties, TPT, RT_{1/2}, and RFI (adjusted *P* values from <0.001 to 0.0187). The second-order nonlinear model produced force too rapidly, resulting in the largest TPT errors (under-predictions) of the three models (Fig. 8*A*, Table 3). No interaction effects were included for TPT, as only the 5-Hz CT force train was used to determine this property. The overall mean RT_{1/2} was most accurately predicted by the nonlinear models: first the Hill-Huxley-type model followed closely by the second-order nonlinear model (all differences significant; Fig. 8*B*, Table 3). The only significant interaction for RT_{1/2} (model × frequency; see Fig. 8*B*) occurred for the directional RT_{1/2} error [*F*(4,12) = 8.85, adjusted *P* value = 0.0477]. Force fusion errors (RFI; Fig. 8*C*) also differed significantly between all three models overall [*F*(2,6) = 8.31, adjusted *P* value = 0.0187 and 114.49, adjusted *P* value <0.0001, absolute and directional, respectively]. The linear model over-predicted the RFI, whereas the two nonlinear models significantly under-predicted the fusion level (Fig. 8*C*). However, despite a significant between-model effect using the repeated-measures ANOVA, individual pairwise comparisons of RFI error magnitude (absolute error) did not reach significance using the conservative Bonferroni correction (Table 3). The model × frequency interaction for RFI was significant (directional and absolute, adjusted *P* values from <0.0001 to 0.007), indicating between-model differences varied with frequency. For example, the linear model, although showing the largest fusion errors overall, produced the lowest mean RFI error at 10 Hz.

The catch-like property of muscle was evaluated using the overall %DDiff error (mean of the normalized initial and mid-doublet differences). All three models had difficulty predicting muscle forces resulting from doublet stimulations located either at the start or within a train of stimuli, particularly at frequencies resulting in partial or full force fusion (e.g., 10 and 20 Hz; Fig. 8*D*). The linear model produced mean %DDiff errors that far exceeded the nonlinear models' errors (Table 3, Fig. 8*D*). This difference between the nonlinear and linear models was driven predominantly by the large linear model error at 20 Hz, where the linear model was incapable of scaling down the additional doublet force, as occurred experimentally (see Fig. 4). This is supported by the significant model × frequency interaction [*F*(4,12) = 40.78, adjusted *P* value = 0.0044]. This was the only specific force property that the second-order nonlinear model provided significantly better predictions than the Hill-Huxley-type nonlinear model (Fig. 8*D*, Table 3), but the mean difference (∼2%) was small relative to the observed mean errors (23 and 25% for the second-order nonlinear and Hill-Huxley-type nonlinear models, respectively).

## DISCUSSION

The primary finding of this study was that chronically paralyzed muscle force-time profiles were most accurately predicted by the nonlinear models [Hill-Huxley-type nonlinear model (18) followed closely by the second-order nonlinear model (5)], supporting our initial hypothesis. However, significant two- and three-way interactions indicate differential errors with stimulation frequency and pattern. The nonlinear models predicted most of the specific force properties (PF, FTI, RT_{1/2}, and catch-like property of muscle) better than the linear model, although TPT and RFI were equally well represented by the linear model. Interestingly, the more complex Hill-Huxley-type nonlinear model predicted only two specific force properties, the TPT and RT_{1/2}, better than the less-complex second-order nonlinear model, contrary to the hypothesis that a Hill-based model might better represent physiological force properties.

#### Specific force property predictions.

Each model had its individual predictive strengths and weaknesses. The similar PF and FTI errors across all models at the lower frequencies suggest that a linear approximation of muscle force magnitude may be reasonable for unfused force trains. Beyond the fusion frequency, the linear model would likely have increasing difficulty adjusting for the nonlinear force saturation. However, for those with chronic paralysis, higher frequencies may be associated with an increased risk of fracture (32), making lower frequency stimulation a more desirable activation strategy.

The frequencies chosen in this study spanned a range of the force frequency curve for human chronically paralyzed soleus muscle (43), with mean experimental RFIs of 8.4, 82.1, and 98.0% (where 100% indicates no discernable drop in force between pulses). Although fusion is a nonlinear relationship, the linear and nonlinear models were equally able to predict RFI for frequencies below full fusion.

Surprisingly, TPT was accurately predicted by the simple linear and the Hill-Huxley-type nonlinear models, whereas the nonlinearities introduced by the second-order nonlinear model (relative to the linear model) actually worsened the TPT estimates (Fig. 8). The linear model had substantially more difficulty predicting RT_{1/2} suggesting muscle force decay involves a more complex, nonlinear process than does force initiation.

Experimentally, the doublet produced considerably greater additional force at 5 Hz than at 20 Hz (see Fig. 4). Both nonlinear models were able to attenuate doublet force output with increasing frequency, unlike the linear model. The second-order nonlinear model predicted this catch-like property of muscle better than the Hill-Huxley-type nonlinear model despite the inclusion of the force-modulating variable parameter, R_{i}, in the Hill-Huxley-type model (15, 16). The frequency-dependent variable parameter, b, in the second-order nonlinear model, although not previously described for this purpose, is able to represent this nonlinear force summation phenomenon. Accordingly, both modeling approaches (i.e., using a force-dependent or a frequency-dependent variable parameter) adequately model the nonlinear nature of the catch-like property of muscle.

#### Previous error findings.

Our observed errors are relatively consistent with previous studies. Ding and colleagues reported *R*^{2} values for the Hill-Huxley-type model ranging from ∼0.64 to 0.90 for PF and FTI comparisons (37) and 0.80 to 0.99 (17) for overall force profiles in human able-bodied quadriceps muscle activated at ∼20% of maximum. We observed similar *R*^{2} values for the Hill-Huxley-type model, ranging from 0.74 to 0.99. In paralyzed muscle, the Hill-Huxley-type model previously produced <10% errors in PF and FTI for 81 and 89% of the test trials (each the mean of two contractions) (14), but the mean or range of error values were not reported. We observed similar errors of 7.2 ± 1.2 and 7.6 ± 1.5% for PF and FTI, respectively.

The second-order nonlinear model produced from 2.4 to 5.5% error in isolated cat muscle (5), which is relatively consistent with our observations in humans (5.4–16.0%). Differences in experimental paradigms, isolated animal muscle vs. in vivo human testing and nonparalyzed vs. paralyzed muscle, may explain these error discrepancies.

Second-order linear models have been used to fit experimental animal (36) and human (1, 3) muscle forces but have been described as having poor predictive capabilities (5, 6, 24). No specific error values were provided for these conclusions. Conversely, we observed reasonable linear predictions of chronically paralyzed human muscle force for low CTs, with poor linear force predictions with higher frequency or variable stimulation.

#### Limitations.

Although spasticity can be problematic in chronically paralyzed muscle, consistent, reproducible force profiles were readily obtained at each frequency with supramaximal nerve stimulation (e.g., Fig. 4), eliminating the need to average trials to reduce noise. These results, however, are limited to supramaximally stimulated, isometric contractions in paralyzed soleus muscle at a single joint angle and may not extend to different joint angles (3), potentiation states, dynamic motions, and/or fatigue (21). The parameter values and errors obtained for the soleus muscle may not be representative of those for other paralyzed muscles, i.e., large, multi-joint muscles. Individualized parameterization techniques for each model may improve their individual abilities to predict paralyzed muscle force; however, using a single method provided a common baseline for comparison. The small effect size between the two nonlinear models for PF, FTI, and RFI may represent the predictive similarities between these two nonlinear models or be a result of inadequate power for these select comparisons. The purpose of this study was to determine the predictive modeling errors associated with a fresh muscle state, but the test protocol stimulations induced a minimal but measurable level of fatigue (≤10 % mean decrease in PF of the doublet ramp train).

#### Considerations for an ideal muscle model.

Specific force properties (PF, TPT, RT_{1/2}, etc.) have been associated with Hill model parameter values more frequently than with other types of modeling strategies. Therefore, one might also expect specific force properties to be most accurately predicted by a Hill-based model. Our findings do not support this assertion. The systems-based second-order nonlinear model, described as a second-order filter by Bobet and Stein (5), was nearly equivalent to the Hill-Huxley-type model for all specific force property measurements, except for TPT. Equating parameter values with physiological muscle properties does not appear to assist in making the most accurate force predictions for human paralyzed muscle.

Although model accuracy is a primary consideration for choosing an appropriate modeling strategy, model complexity should be considered as well. Based on the law of parsimony (38), a preferred model would use the simplest means to produce a reasonable force representation. Given the similar predictive capabilities of the two nonlinear models relative to their difference in complexities, coupled with the larger errors associated with the linear model, the simplest model to best predict human chronically paralyzed muscle force is the second-order nonlinear model. However, the linear model's inaccuracies could be tolerable if its simplicity is required, i.e., a complex control strategy is involved (24). Furthermore, the relatively small differences in model accuracy between the two nonlinear models may be important depending on the intended application (e.g., where accurate TPT is crucial).

In conclusion, the force predictions of three different muscle modeling approaches were relatively similar, suggesting that physiological muscle force can be adequately represented mathematically in multiple ways. Indeed, the interactions between model, frequency (5, 10, and 20 Hz), and stimulation pattern (CT vs. DDT) indicate the difficulty in assessing the differences between models, as the models' accuracies varied between stimulation trains. Any one of the three models could prove to be an ideal choice, depending on the intended therapeutic application, and thus the frequency and/or stimulation patterns used. Considering the three models' predictive errors, we conclude that *1*) the linear model is best able to predict paralyzed muscle forces using constant, low-frequency stimulation trains; *2*) the two nonlinear models provide nearly equivalent force predictions despite their differences in modeling strategies; and *3*) the nonlinear models are most accurate when using higher frequency or variable stimulation trains for predicting electrically activated muscle force in chronically paralyzed human muscle.

## GRANTS

This work was supported in part by an award (R01-HD-39445) from the National Center for Medical Rehabilitation Research and the Christopher Reeve Foundation.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2006 the American Physiological Society