Journal of Applied Physiology Journal of Neurophysiology
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

J Appl Physiol 100: 1027-1036, 2006. First published November 23, 2005; doi:10.1152/japplphysiol.00935.2005
8750-7587/06 $8.00
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF) Free
Right arrow All Versions of this Article:
100/3/1027    most recent
00935.2005v1
Right arrow Submit a response
Right arrow Alert me when this article is cited
Right arrow Alert me when eLetters are posted
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Frey Law, L. A.
Right arrow Articles by Shields, R. K.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Frey Law, L. A.
Right arrow Articles by Shields, R. K.

Predicting human chronically paralyzed muscle force: a comparison of three mathematical models

Laura A. Frey Law and Richard K. Shields

Graduate Program in Physical Therapy and Rehabilitation Science, The University of Iowa, Iowa City, Iowa

Submitted 1 August 2005 ; accepted in final form 22 November 2005


   ABSTRACT
TOP
 ABSTRACT
METHODS
 RESULTS
 DISCUSSION
 GRANTS
 REFERENCES
 
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
TOP
 ABSTRACT
 METHODS
RESULTS
 DISCUSSION
 GRANTS
 REFERENCES
 
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 [{delta}(t)], denoted by r(t) = {Sigma} {delta}(ti – t0) for i = 1 to n. The {delta}(t) function equals 1/dt at the time (ti) of each pulse, for a total area equal to 1.0, and zero otherwise.

Formula 1(1)
Essentially, the forces resulting from each consecutive pulse are summed linearly to produce the overall force train. The three constant parameters involved, beta, {zeta}, and {omega}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, {zeta}, 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 {omega}n and {zeta} determines the overall system time constant (22). For modeling human muscle, values for {zeta} are typically ≤1.0 (critically damped) (1, 3, 5). The natural frequency, {omega}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 {omega}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 b0, b1, and B (5). When b1 is positive, this variable rate constant becomes slower with increasing force and faster when b1 is negative. Parameter B acts as a simple gain factor in this model (27).

Formula 2(2)

Formula 3(3)

Formula 4(4)

Formula 5(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, b0, b1, 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, {tau}1, {tau}c, {tau}2, R0, and km, 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 Ca2+ as well as troponin binding) (16, 18, 21). This equation is modulated by the variable parameter, Ri (Eq. 8), which decays as a function of each interpulse interval (ti to ti-1) from a maximum value of R0 for very rapid pulses to a minimum of 1 for widely spaced pulses. An analytical solution for the first differential equation (CN) 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.

Formula 6(6)

Formula 7(7)

Formula 8(8)

Formula 9(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.


View this table:
[in this window]
[in a new window]
  		Table 1. Subject characteristics

 
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.


Figure 1
View larger version (11K):
[in this window]
[in a new window]
  		Fig. 1. Schematic representation of the stimulation patterns, showing only the 10-Hz series and the doublet ramp used for parameter determination for each model. CT, constant-frequency train; DT, single doublet train; DDT, dual doublet train.

 
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 (RT1/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, RT1/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 RT1/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 x 2 x 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 RT1/2, RFI, and %DDiff were analyzed using a two-way (3 x 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 {alpha} = 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
TOP
 ABSTRACT
 METHODS
 RESULTS
DISCUSSION
 GRANTS
 REFERENCES
 
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.


Figure 2
View larger version (31K):
[in this window]
[in a new window]
  		Fig. 2. Representative examples of the raw data from a single subject (subject 4) resulting from the doublet ramp (A), the 5-Hz DT (B), the 10-Hz DDT (C), and the 20-Hz CT (D). EMG, electromyogram. *Doublet pulses.

 

Figure 3
View larger version (26K):
[in this window]
[in a new window]
  		Fig. 3. Representative example of the optimized best fits for the 3 models to the experimental doublet ramp (solid line) train for subject 1: linear (dotted line), 2nd-order nonlinear (NL; dash-dot line), and the Hill-Huxley-type NL (dashed line).

 

Figure 4
View larger version (35K):
[in this window]
[in a new window]
  		Fig. 4. Experimental force (first row) and modeled forces (second to fourth rows) for chronically paralyzed soleus muscle (subject 4) at 5, 10, and 20 Hz. The linear model is shown in row 2, 2nd-order NL model in row 3, and Hill-Huxley-type model in row 4. CT and DDT stimulation patterns are displayed at left and right, respectively.

 

View this table:
[in this window]
[in a new window]
  		Table 2. Optimized parameter values

 
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. 5A) 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. 5B] as well as significant interactions for model by frequency [F(4,12) = 13.18, adjusted P value = 0.0127; Fig. 5C], model by stimulation pattern [F(2,6) = 15.47, adjusted P value = 0.0161; Fig. 5D], 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.


Figure 5
View larger version (25K):
[in this window]
[in a new window]
  		Fig. 5. Overall %error between modeled and experimental forces. A: mean (SD) %error by subject. B: mean (SE) overall %error. C: mean (SE) %error by frequency. D: mean (SE) %error by stimulation pattern.

 

View this table:
[in this window]
[in a new window]
  		Table 3. Error comparisons

 
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 x frequency, model x stimulation pattern, and model x frequency x stimulation pattern) resulted for all directional and absolute %PF and %FTI errors, except the model x 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.


Figure 6
View larger version (19K):
[in this window]
[in a new window]
  		Fig. 6. Overall %peak force (PF) error between modeled and experimental forces. A: mean (SD) PF %error by subject. B: mean (SE) overall PF %error. C: mean (SE) PF %error by frequency. D: mean (SE) PF %error by stimulation pattern. Note that mean errors may be approaching zero but have clearly visible variances.

 

Figure 7
View larger version (19K):
[in this window]
[in a new window]
  		Fig. 7. Overall %force time integral (FTI) error between modeled and experimental forces. A: mean (SD) FTI %error by subject. B: mean (SE) overall FTI %error. C: mean (SE) FTI %error by frequency. D: mean (SE) FTI %error by stimulation pattern. Note that mean errors may be approaching zero but have clearly visible variances.

 
Significant between-model differences were observed for the absolute and directional errors for the force time properties, TPT, RT1/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. 8A, 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 RT1/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. 8B, Table 3). The only significant interaction for RT1/2 (model x frequency; see Fig. 8B) occurred for the directional RT1/2 error [F(4,12) = 8.85, adjusted P value = 0.0477]. Force fusion errors (RFI; Fig. 8C) 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. 8C). 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 x 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.


Figure 8
View larger version (21K):
[in this window]
[in a new window]
  		Fig. 8. Specific force time errors between modeled and experimental forces by frequency. A: mean (SE) time to peak tension (TPT) for 5 Hz only. B: mean (SE) half relaxation time (1/2RT). C: mean (SE) relative fusion index (RFI). D: mean (SE) doublet difference %error (DDiff), an indicator of the catch-like property of muscle.

 
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. 8D). The linear model produced mean %DDiff errors that far exceeded the nonlinear models' errors (Table 3, Fig. 8D). 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 x 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. 8D, 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
TOP
 ABSTRACT
 METHODS
 RESULTS
 DISCUSSION
GRANTS
 REFERENCES
 
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, RT1/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 RT1/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 RT1/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, Ri, 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 R2 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 R2 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, RT1/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
TOP
 ABSTRACT
 METHODS
 RESULTS
 DISCUSSION
 GRANTS
REFERENCES
 
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
 

Address for reprint requests and other correspondence: R. K. Shields, 1-252 Medical Education Bldg., The Univ. of Iowa, Iowa City, IA 52242 (e-mail: richard-shields{at}uiowa.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. Section 1734 solely to indicate this fact.


   REFERENCES
TOP
 ABSTRACT
 METHODS
 RESULTS
 DISCUSSION
 GRANTS
 REFERENCES
 

  1. Aaron SL and Stein RB. Comparison of an EMG-controlled prosthesis and the normal human biceps brachii muscle. Am J Phys Med 55: 1–14, 1976.[Web of Science][Medline]
  2. Baratta RV, Zhou BH, and Solomonow M. Frequency response model of skeletal muscle: effect of perturbation level, and control strategy. Med Biol Eng Comput 27: 337–345, 1989.[CrossRef][Web of Science][Medline]
  3. Bawa P and Stein RB. Frequency response of human soleus muscle. J Neurophysiol 39: 788–793, 1976.[Abstract/Free Full Text]
  4. Bobet J. Can muscle models improve FES-assisted walking after spinal cord injury? J Electromyogr Kinesiol 8: 125–132, 1998.[CrossRef][Web of Science][Medline]
  5. Bobet J and Stein RB. A simple model of force generation by skeletal muscle during dynamic isometric contractions. IEEE Trans Biomed Eng 45: 1010–1016, 1998.[CrossRef][Web of Science][Medline]
  6. Bobet J, Stein RB, and Oguztoreli MN. A linear time-varying model of force generation in skeletal muscle. IEEE Trans Biomed Eng 40: 1000–1006, 1993.[CrossRef][Web of Science][Medline]
  7. Burnham R, Martin T, Stein R, Bell G, MacLean I, and Steadward R. Skeletal muscle fibre type transformation following spinal cord injury. Spinal Cord 35: 86–91, 1997.[CrossRef][Web of Science][Medline]
  8. Castro MJ, Apple DF Jr, Hillegass EA, and Dudley GA. Influence of complete spinal cord injury on skeletal muscle cross-sectional area within the first 6 months of injury. Eur J Appl Physiol Occup Physiol 80: 373–378, 1999.[CrossRef][Medline]
  9. Castro MJ, Apple DF Jr, Rogers S, and Dudley GA. Influence of complete spinal cord injury on skeletal muscle mechanics within the first 6 months of injury. Eur J Appl Physiol 81: 128–131, 2000.[CrossRef][Web of Science][Medline]
  10. Castro MJ, Apple DF Jr, Staron RS, Campos GE, and Dudley GA. Influence of complete spinal cord injury on skeletal muscle within 6 mo of injury. J Appl Physiol 86: 350–358, 1999.[Abstract/Free Full Text]
  11. Chou CP and Hannaford B. Dual stable point model of muscle activation and deactivation. Biol Cybern 66: 511–523, 1992.[CrossRef][Web of Science][Medline]
  12. Close CM and Frederick DK. Modeling and Analysis of Dynamic Systems. New York: John Wiley & Sons, 1995.
  13. Ding J, Binder-Macleod SA, and Wexler AS. Two-step, predictive, isometric force model tested on data from human and rat muscles. J Appl Physiol 85: 2176–2189, 1998.[Abstract/Free Full Text]
  14. Ding J, Lee SC, Johnston TE, Wexler AS, Scott WB, and Binder-Macleod SA. Mathematical model that predicts isometric muscle forces for individuals with spinal cord injuries. Muscle Nerve 31: 702–712, 2005.[CrossRef][Web of Science][Medline]
  15. Ding J, Storaska JA, and Binder-Macleod SA. Effect of potentiation on the catchlike property of human skeletal muscles. Muscle Nerve 27: 312–319, 2003.[CrossRef][Web of Science][Medline]
  16. Ding J, Wexler AS, and Binder-Macleod SA. Development of a mathematical model that predicts optimal muscle activation patterns by using brief trains. J Appl Physiol 88: 917–925, 2000.[Abstract/Free Full Text]
  17. Ding J, Wexler AS, and Binder-Macleod SA. A mathematical model that predicts the force-frequency relationship of human skeletal muscle. Muscle Nerve 26: 477–485, 2002.[CrossRef][Web of Science][Medline]
  18. Ding J, Wexler AS, and Binder-Macleod SA. Mathematical models for fatigue minimization during functional electrical stimulation. J Electromyogr Kinesiol 13: 575–588, 2003.[CrossRef][Web of Science][Medline]
  19. Ding J, Wexler AS, and Binder-Macleod SA. A predictive fatigue model–I: predicting the effect of stimulation frequency and pattern on fatigue [erratum appears in IEEE Trans Neural Syst Rehabil Eng 11: 86, 2003]. IEEE Trans Neural Syst Rehabil Eng 10: 48–58, 2002.[CrossRef]
  20. Ding J, Wexler AS, and Binder-Macleod SA. A predictive fatigue model—II: predicting the effect of resting times on fatigue. IEEE Trans Neural Sys Rehabil Eng 10: 59–67, 2002.[CrossRef]
  21. Ding J, Wexler AS, and Binder-Macleod SA. A predictive model of fatigue in human skeletal muscles. J Appl Physiol 89: 1322–1332, 2000.[Abstract/Free Full Text]
  22. DiStefano JJ III, Stubberud AR, and Williams IJ. Feedback and Control Systems. New York: McGraw-Hill, 1990.
  23. Dorgan SJ and O'Malley MJ. A nonlinear mathematical model of electrically stimulated skeletal muscle. IEEE Trans Rehabil Eng 5: 179–194, 1997.[CrossRef][Medline]
  24. Durfee WK and Palmer KI. Estimation of force-activation, force-length, and force-velocity properties in isolated, electrically stimulated muscle. IEEE Trans Biomed Eng 41: 205–216, 1994.[CrossRef][Web of Science][Medline]
  25. Epstein M and Herzog W. Theoretical Models of Skeletal Muscle. Chicester, UK: John Wiley & Sons, 1998.
  26. Frey Law LA and Shields RK. Femoral loads during passive, active, and active-resistive stance after spinal cord injury: a mathematical model. Clin Biomech 19: 313–321, 2004.[CrossRef][Medline]
  27. Frey Law LA and Shields RK. Mathematical models use varying parameter strategies to represent paralyzed muscle force properties: a sensitivity analysis. J Neuroeng Rehabil 2: 12, 2005.[CrossRef]
  28. Gerrits HL, De Haan A, Hopman MT, van Der Woude LH, Jones DA, and Sargeant AJ. Contractile properties of the quadriceps muscle in individuals with spinal cord injury. Muscle Nerve 22: 1249–1256, 1999.[CrossRef][Web of Science][Medline]
  29. Gollee H, Murray-Smith DJ, and Jarvis JC. A nonlinear approach to modeling of electrically stimulated skeletal muscle. IEEE Trans Biomed Eng 48: 406–415, 2001.[CrossRef][Web of Science][Medline]
  30. Grimby G, Broberg C, Krotkiewska I, and Krotkiewski M. Muscle fiber composition in patients with traumatic cord lesion. Scand J Rehabil Med 8: 37–42, 1976.[Web of Science][Medline]
  31. Hannaford B. A nonlinear model of the phasic dynamics of muscle activation. IEEE Trans Biomed Eng 37: 1067–1075, 1990.[CrossRef][Web of Science][Medline]
  32. Hartkopp A, Murphy RJ, Mohr T, Kjaer M, and Biering-Sorensen F. Bone fracture during electrical stimulation of the quadriceps in a spinal cord injured subject. Arch Phys Med Rehabil 79: 1133–1136, 1998.[CrossRef][Web of Science][Medline]
  33. Hill VA. The heat of shortening and the dynamic constants of muscle. Proc R Soc Lond B Biol Sci 126: 136–195, 1938.[Free Full Text]
  34. Huxley AF. Muscle structure and theories of contraction. Prog Biophys Biophys Chem 7: 255–318, 1957.[Medline]
  35. Huynh H and Feldt LS. Estimation of the box correction for degrees of freedom from sample data in the randomized block and split plot designs. J Edu Stat 1: 69–82, 1976.
  36. Mannard A and Stein RB. Determination of the frequency response of isometric soleus muscle in the cat using random nerve stimulation. J Physiol 229: 275–296, 1973.[Abstract/Free Full Text]
  37. Perumal R, Wexler AS, Ding J, and Binder-Macleod SA. Modeling the length dependence of isometric force in human quadriceps muscles. J Biomech 35: 919–930, 2002.[CrossRef][Web of Science][Medline]
  38. Phillips CA, Repperger DW, Neidhard-Doll AT, and Reynolds DB. Biomimetic model of skeletal muscle isometric contraction: I. An energetic-viscoelastic model for the skeletal muscle isometric force twitch. Comput Biol Med 34: 307–322, 2004.[CrossRef][Web of Science][Medline]
  39. Round JM, Barr FM, Moffat B, and Jones DA. Fibre areas and histochemical fibre types in the quadriceps muscle of paraplegic subjects. J Neurol Sci 116: 207–211, 1993.[CrossRef][Web of Science][Medline]
  40. Scelsi R, Marchetti C, Poggi P, Lotta S, and Lommi G. Muscle fiber type morphology and distribution in paraplegic patients with traumatic cord lesion. Histochemical and ultrastructural aspects of rectus femoris muscle. Acta Neuropathol (Berl) 57: 243–248, 1982.[CrossRef][Medline]
  41. Shields RK. Fatigability, relaxation properties, and electromyographic responses of the human paralyzed soleus muscle. J Neurophysiol 73: 2195–2206, 1995.[Abstract/Free Full Text]
  42. Shields RK. Muscular, skeletal, and neural adaptations following spinal cord injury. J Orthop Sports Phys Ther 32: 65–74, 2002.[Web of Science][Medline]
  43. Shields RK and Chang YJ. The effects of fatigue on the torque frequency curve of the human paralyzed soleus muscle. J Electromyogr Kinesiol 7: 3–13, 1997.
  44. Shields RK and Dudley-Javoroski S. Musculoskeletal deterioration and hemicorporectomy after spinal cord injury. Phys Ther 83: 263–275, 2003.[Abstract/Free Full Text]
  45. Shields RK, Law LF, Reiling B, Sass K, and Wilwert J. Effects of electrically induced fatigue on the twitch and tetanus of paralyzed soleus muscle in humans. J Appl Physiol 82: 1499–1507, 1997.[Abstract/Free Full Text]
  46. Stein RB and Oguztoreli MN. A model of whole muscles incorporating functionally important nonlinearities. In: Nonlinear Phenomena In Mathematical Sciences, edited by Lakshmikantham V. New York: Academic Press, 1982, p. 749–766.
  47. Talmadge RJ, Castro MJ, Apple DF Jr, and Dudley GA. Phenotypic adaptations in human muscle fibers 6 and 24 wk after spinal cord injury. J Appl Physiol 92: 147–154, 2002.[Abstract/Free Full Text]
  48. Wexler AS, Ding J, and Binder-Macleod SA. A mathematical model that predicts skeletal muscle force. IEEE Trans Biomed Eng 44: 337–348, 1997.[CrossRef][Web of Science][Medline]
  49. Zahalak GI. The two-state cross-bridge model of muscle is an asymptotic limit of multi-state models. J Theor Biol 204: 67–82, 2000.[CrossRef][Web of Science][Medline]
  50. Zahalak GI and Ma SP. Muscle activation and contraction: constitutive relations based directly on cross-bridge kinetics. J Biomech Eng 112: 52–62, 1990.[Web of Science][Medline]



This article has been cited by other articles:


Home page
 J. Appl. Physiol.Home page
R. K. Shields, S. Dudley-Javoroski, and K. R. Cole
Feedback-controlled stimulation enhances human paralyzed muscle performance
J Appl Physiol, November 1, 2006; 101(5): 1312 - 1319.
[Abstract] [Full Text] [PDF]

Home page
 J. Appl. Physiol.Home page
R. K. Shields, S. Dudley-Javoroski, and A. E. Littmann
Postfatigue potentiation of the paralyzed soleus muscle: evidence for adaptation with long-term electrical stimulation training
J Appl Physiol, August 1, 2006; 101(2): 556 - 565.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF) Free
Right arrow All Versions of this Article:
100/3/1027    most recent
00935.2005v1
Right arrow Submit a response
Right arrow Alert me when this article is cited
Right arrow Alert me when eLetters are posted
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Frey Law, L. A.
Right arrow Articles by Shields, R. K.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Frey Law, L. A.
Right arrow Articles by Shields, R. K.

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Visit Other APS Journals Online
Copyright © 2006 by the American Physiological Society.