Functional electrical stimulation can assist paralyzed individuals to perform functional movements, but muscle fatigue is a major limitation to its practical use. An accurate and predictive mathematical model can facilitate the design of stimulation patterns that optimize aspects of the force transient while minimizing fatigue. Solution nonuniqueness, a major shortcoming in previous work, was overcome with a simpler model. The model was tested on data collected during isometric contractions of rat gastrocnemius muscles and human quadriceps femoris muscles under various physiological conditions. For each condition tested, parameter values were identified using the force response to one or two stimulation trains. The parameterized model was then used to predict forces in response to other stimulation patterns. The predicted forces closely matched the measured forces. The model was not sensitive to initial parameter estimates, demonstrating solution uniqueness. By predicting the force that develops in response to an arbitrary pattern of stimulation, we envision the present model helping identify optimal stimulation patterns for activation of skeletal muscle during functional electrical stimulation.
- functional electrical stimulation
- mathematical model
- catchlike property
- muscle length
functional electrical stimulation (FES) can assist paralyzed individuals to perform functional movement (18, 20). Muscle fatigue, however, is a major limitation to its practical implementation (18, 20). It has long been known that the stimulation frequency affects the forces produced and the rate of muscle fatigue (8). More recently, however, it has been shown that, in addition to the mean frequency, the pattern of the pulses can markedly affect force production and fatigue (4, 5). Studies showed that the activation pattern that produces the greatest force varies from person to person (19) and, even for the same individual, changes for different physiological conditions [e.g., muscle length (21)]. Therefore, it is important and difficult to identify activation patterns that maximize force and minimize fatigue for individual patients.
The development of a mathematical model that predicts the force responses of muscles to electrical stimulation would enable clinicians to identify the optimal activation patterns for individual subjects. Several mathematical models of muscle force have been developed, including pure mathematical models (9, 12), Huxley-type models (17, 26,27, 32), and Hill-type models (11, 13, 16, 22, 31). Different aspects of the muscle have been modeled. Our long-term goal is to develop a model that can predict forces of muscle during FES. A reliable isometric force model is a crucial step toward this goal. The only model able to predict isometric force response to a wide range of stimulation patterns incorporates the Hill paradigm (31). Hill-type models are phenomenological, originating from experimental observations of the muscle’s mechanical behavior, and describe the whole muscle as a contractile element in series with a spring, together with various modifications such as dampers placed in parallel or series with other components.
In our previous work we decomposed the contractile response of skeletal muscle into three steps: 1) Ca2+ release and uptake by the sarcoplasmic reticulum, 2) Ca2+ reaction with troponin, and3) force mechanics modeled by a linear spring, damper, and motor in series (31). Steps 1 and 2 make the model less phenomenological than Hill’s original model (16). This three-step model consisted of three ordinary differential equations and nine parameters that were identified by fitting the model to forces measured in response to a single constant-frequency train (CFT) stimulation. Then, using the parameter values identified from this fit, we tested the fidelity of the model by comparing its predictions with responses measured from other stimulation patterns. For rat gastrocnemius and soleus muscles, the model accurately fit and predicted the force responses during stimulation with a wide range of stimulation frequencies and patterns. The primary shortcoming of the model was its sensitivity to different initial parameter estimates. Numerous fitting trials were necessary to identify the initial parameter values, leading to identified values that reliably predicted force responses. In addition, subsequent analysis of these previously published results showed that the model was only able to predict the stimulation pattern that produced the greatest force-time integral for 9 of the 15 gastrocnemius muscles studied. All the muscles for which the model was successful had twitch contraction times <25 ms. The purpose of this study, therefore, was to develop a model with greater reliability and range of applicability.
The simpler model presented here consists of only two ordinary differential equations and four parameters. It has been tested on isometric contractions from rat and human muscles in response to a wide range of brief, physiologically relevant stimulation patterns. Forces used to test the model were collected from fresh and fatigued muscles and at long (near-optimal) and short (suboptimal) muscle length.
FORMULATION OF THE MODEL
In the new muscle model the contractile response is decomposed into two distinct physiological steps: cross-bridge activation and force generation.
Pilot studies with our previous model (31) showed that the ability of the model to predict forces is relatively insensitive to the specific shape and magnitude of the Ca2+and Ca2+-troponin complex transient. This implies that steps 1and 2 in our previous model could be combined into one step. Therefore, we use one unitless factor,CN , the dynamics of which are governed by a time constant, τc, to qualitatively describe the rate-limiting step before the myofilaments mechanically slide across each other and generate force (1; see for details). The differential equation describing this dynamic is Equation 1 and the chosen analytic solution is Equation 2 where the transient behavior ofCN characterizes the rate-limiting step in force generation,t (ms) is the time since the beginning of the stimulation train,ti (ms) is the time of the ith stimulation pulse since the beginning of the stimulation train,n is the number of stimulation pulses before time t in the train (its maximum value is 6 in this study), and τc (ms) is the time constant controlling the transient shape ofCN .
Similar to our previous model (31), the force generation is modeled mechanically by a linear spring, damper, and motor in series, which represent the elastic, viscous, and contractile elements of the muscle, respectively (Fig. 1). The damper represents the viscous resistance of the contractile and connective tissue (22). The force (F) exerted by the damper is given by Equation 3 whereb is the damping coefficient,x is the length of the spring, andV is the contractile velocity of the motor.
The motor represents the contractile element, or the sliding of actin and myosin filaments, and is driven by strongly bound cross bridges (1). Because there is a sigmoidal force-pCa relationship (2), we tested for a similar relationship betweenCN and the strongly bound cross bridges and found that a simple Michaelis-Menten function, CN /(1 + CN ), sufficed. The contractile capability of the muscle,V, therefore, is Equation 4 in which B is a proportionality constant.
A linear spring represents the series elastic components of the muscle and tendon. The force exerted by the spring is given by Equation 5 whereK is the spring constant. Differentiating Eq. 5 with respect to time and using Eq. 3 to eliminate dx/dtand Eq. 4 to eliminateV gives Equation 6The term b/K represents the time constant over which the force decays. The friction between the actin and myosin filaments increases with an increase in bound cross bridges, so we set b/K = τ1 + τ2 CN /(1 + CN ), where τ1 is the value of the time constant in the absence ofCN and τ2 is the additional frictional component due to the cross-bridge chemical bonds. Substituting forb/K and replacingKB with a new constant,A, in Eq.6 gives Equation 7 The isometric force dynamics are governed by Eqs.1 and 7 , which describe the transient behavior of the two state variables,CN and F, subject to the four parameters τc,A, τ1, and τ2 (Table1).
Rat and human data were taken from previously published or unpublished studies in our laboratory, which were approved by the University of Delaware’s Animal and Human Subjects Review Board. Data were collected from 16 nonfatigued and 8 fatigued rat gastrocnemius muscles. The experimental methods are similar to those previously described (31). Briefly, the rats were deeply anesthetized, and a hook stimulating electrode was positioned around the nerve leading to the gastrocnemius muscle. Force was recorded directly from the Achilles’ tendon. The gastrocnemius muscle was stimulated with a series of nine six-pulse trains. During traditional FES, skeletal muscle is activated with short CFTs, where the pulses within each train are separated by regular intervals. Five of the trains in our study, therefore, were CFTs with interpulse intervals (IPIs) of 10, 20, 30, 40, or 50 ms (Fig.2 A): CFT10, CFT20, CFT30, CFT40, and CFT50, respectively. The other four trains were variable-frequency trains (VFTs), with an initial IPI of 10 ms, and the remaining four IPIs within the train were 20, 30, 40, or 50 ms (Fig. 2 B): VFT20, VFT30, VFT40, and VFT50, respectively. The VFTs were chosen, because recently it has been shown that VFTs that take advantage of the catchlike property of skeletal muscle can markedly increase the force output from fatigued muscle compared with CFT stimulation (4, 5, 8). The catchlike property is the tension enhancement seen when an initial brief high-frequency burst of pulses (2–4 pulses) at the onset of a subtetanic CFT is used to activate the muscle (6, 10) and is a fundamental muscle property that is not attributable to the motor axon or neuromuscular junction (4, 10). To collect the responses to nonfatiguing stimulation, ∼10 s separated each train. To fatigue the muscle, 3 random sequences of these 9 trains were used to form a cycle of 27 trains, and 8 such cycles (216 trains) were used. Each train was separated by ∼1 s. The force responses to the last cycle of trains were used to represent the force responses from the fatigued muscles.
The experimental setup for collection of human data was similar to that previously described (8). Briefly, subjects were seated on a force dynamometer with their hips flexed to ∼75° and their knees flexed to 90° for the long-muscle-length study or 15° for the short-muscle-length study. Two stimulating electrodes were placed over the quadriceps femoris muscle. The force transducer was positioned against the anterior aspect of the leg, proximal to the lateral malleolus. Nonfatiguing and fatiguing protocols, similar to those used for rat muscles, were used for human subjects, except the initial IPI of the VFTs was 5 ms, and the IPIs of the CFTs were 10–120 ms in steps of 10 ms, for a total of 24 stimulation patterns. Nonfatiguing force responses were collected for long-muscle-length (n = 12) and short-muscle-length (n = 11) studies. Fatigued force data were collected only at long muscle length (n = 11).
Reliability of the model was first investigated for rat and human muscles by using the model to fit the force data for each animal or subject by incorporation of a wide range of initial parameter values. Convergence to the same parameter values for each rat or human subject under each condition would suggest a unique solution of the model within the range of initial parameter values tested.
SAAM II, a modeling system developed by the University of Washington and the National Institutes of Health, was used to identify the optimum parameter values and predict forces for other stimulation trains (24). For each muscle, τ1 was first calculated using the numerical module of SAAM II. The procedure to calculate τ1 was described in detail previously (31). Briefly, it was found by settingCN = 0, so thatEq. 7 became dF/dt = −F/τ1. By taking the force decay at the end of the tetany and performing a linear regression of lnF vs. t, we obtained the value of τ1 from the slope of the fit. The other three parameter values were then identified using the compartment module of SAAM II, which accepted as input1) Eqs.1 and 7 expressed compartmentally, 2) force measured for a contraction, 3) the stimulation pattern, and 4) initial parameter estimates. SAAM II performs a multidimensional fit between the measured and predicted force data to identify the best parameter values.
Preliminary studies were conducted to find the best train(s) for parameter identification. For rat and human muscles, trains that produce nonfused tetanic force responses were tested first, because they exhibited a clear rise and fall of force for each pulse and a buildup of force over multiple pulses. CFTs with IPIs of 30 ms (CFT30, Fig. 3 C) and 70 ms (CFT70, see Fig. 6 C) were initially used to fit the model for rat gastrocnemius and human quadriceps femoris muscles, respectively.
For rat nonfatigued or fatigued muscles, the 30-ms CFT did not always lead to a set of parameter values that successfully predicted the forces from other stimulation patterns. In fact, no other train was found that could satisfactorily predict forces from all other stimulation patterns. Further testing showed that only a combination of two stimulation trains, one that produced a nearly fused tetanic response (i.e., CFT20) and another that produced an unfused tetanic response (i.e., CFT50), identified a reliable set of parameter values (see Fig. 3 for examples of force responses to CFT20 and CFT50). The CFT20 and CFT50, therefore, were used to parameterize the model for the rat muscles.
Unlike the rat muscles, parameter values for human quadriceps femoris muscles obtained by fitting the combination of a high- and a low-frequency CFT (e.g., CFT30 and CFT100) did not improve the accuracy of the predicted responses. When muscles were at a long length and not fatigued, the CFT with IPI of 70 ms (see Fig.6 C) yielded parameter values that best predicted the forces. However, when muscles were at long length and fatigued or at short length and not fatigued, we found that the CFT70 could only be used to identify parameter values to predict forces to other CFTs and that the VFT70 (see Fig.7 H) best identified parameter values for the VFTs.
For each muscle, with use of the parameterized model, the state variable CN and the force responses to other stimulation patterns were predicted by simply changing the simulation pattern input to the parameterized compartment module. We then tested the fidelity of the model by comparing the predicted forces with the measured forces.
To validate the model, several comparisons were made between the predicted and measured data. To compare the shape of the forces, a correlation coefficient was calculated by comparing the measured and predicted forces at 2-ms-interval time steps. Any phase shift or disagreement in shape between the predicted and measured data lowers the correlation coefficient. Next, for each animal or subject, the force-time integral was compared for the measured and the predicted responses to each stimulation pattern. For each condition, the measured force-time integrals were plotted against the predicted ones. Linear regression trend lines, with intercepts of zero, were used to compare force response across all the stimulation patterns and determine how well the model predicts the force-time integrals. A perfectly accurate model would have a slope and correlation coefficient of 1. Also within each condition, percent differences were calculated between the average measured and predicted force-time integral for each stimulation pattern. If the percent differences were greater than a chosen criterion, 5%, a paired t-test was used to check whether the difference was significant across the subjects. We chose 5%, because we believed that this was the minimal difference that would be clinically significant.
The best-fit parameter values were pooled for all subjects within each condition, and means ± SD were calculated. For each parameter, a paired t-test was used to compare the parameter values for the same animal or subject across conditions, with nonfatigued and long-muscle-length-nonfatigued data as the control condition for rat and human studies, respectively.
Pilot studies on reliability of the model showed convergence of the parameter values within the range of values tested for each of the four rat and six human muscles tested, suggesting that the model was able to fit the data with a unique solution. Subjective evaluation of the force responses showed that the model accurately fit and predicted the responses to a wide range of stimulation patterns for the nonfatigued and fatigued rat muscles (Figs. 3 and 4, respectively). The median correlation coefficients were 0.98 for the nonfatigued and fatigued muscles (Table 2), which is consistent with the subjective impression that the predicted response accurately matched the shape of the measured data. Typical trend lines for rat gastrocnemius force-time integrals had slope and correlation coefficient close to 1 (Fig. 5,A andD), which also suggested an accurate fit and prediction with the model. For nonfatigued rat muscles, the model significantly overestimated the CFT10 force-time integral by 16% (Fig. 5 B). For fatigued rat muscles, although the model significantly overestimated the VFT30, VFT40, and VFT50, the differences between the measured and predicted forces were small (8.1, 6.7, and 8.5%, respectively; Fig.5 F). The model successfully predicted the activation patterns that produced the greatest force-time integrals for seven of the nine rat muscles with twitch contraction times <25 ms and for seven of the eight muscles with twitch contraction times >25 ms.
For human studies, subjective evaluation of the force responses also showed that, in general, the model accurately fit and predicted the responses to a wide range of stimulation patterns and conditions (Figs.6-8). The correlation coefficients between the measured and predicted forces, averaged across subjects, had a median of 0.96, 0.96, and 0.97 (Table2) for nonfatigued long-length, fatigued long-length, and nonfatigued short-length conditions, respectively, suggesting good agreement between the shape of the measured and predicted forces.
When muscles were not fatigued and were held at long length, the trend lines for force-time integral had slopes close to 1, and the model accounted for 65–68% of the variability in the measured data (Fig. 9). The model significantly overestimated the force produced in response to the CFTs and VFTs with short IPIs (10–60 ms) and underestimated those with long IPIs (80–110 ms).
Surprisingly, when muscles were fatigued the model did a better job of predicting the force output than when muscles were not fatigued. The slopes of force-time integral trend lines were only slightly <1, and the model accounted for 90–96% of the variability in the measured data (Fig. 10). The model tended to underestimate the force-time integrals of CFTs with short IPIs (20–40 ms), the percent error being worst (25%) for CFT10. The model slightly overestimated the force-time integrals of the VFTs with short IPIs (20–40 ms) and underestimated the VFTs with long IPIs (90–110 ms).
When muscles were not fatigued and were at short length, trend lines of the force-time integral had slope values close to 1 and correlation coefficients of 0.89–0.94 (Fig. 11). The model significantly overestimated the force-time integral of VFTs with short IPIs (20–40 ms; Fig. 11). However, these differences were <8%.
Overall, the model did an excellent job of predicting the pattern of activation that produced the maximum force-time integral. The only exceptions were the human nonfatigued data collected at long muscle length (Fig. 9). In this condition the error between the maximum predicted and measured force-time integrals was, however, only ∼5%. The model predicted that the maximum force-time integral is obtained with IPIs that are short compared with the measurements, leading to an actual force-time integral that is 10% lower than the measured maximum.
Fatigue produced changes in all four parameter values that were consistent for rat and human muscles (Table 3). Pairedt-test showed a significant increase in τ1 and a significant decrease in A when muscles were fatigued or held at short length for rats and humans. Only for humans, however, were significant differences seen for τ2 and τc across conditions.
The current model predicted isometric forces for rat gastrocnemius and human quadriceps femoris muscles with reasonable accuracy during brief subtetanic and tetanic isomeric contractions by using CFTs and VFTs. It successfully predicted nonfatigued force responses when muscles were held at long or short length and fatigued force responses at long length. Subjective and objective evaluations of the force response showed that the shape of the predicted force responses closely matched the responses that were measured. Comparisons between the predicted and measured force-time integrals also suggested excellent agreement between the predicted and measured forces. Although significant differences between the predicted and measured force-time integrals were observed with some stimulation patterns, most differences were small.
Comparison with previous work.
The present model differs from most previously developed models, in that it is simpler in structure and has fewer parameters (11, 17, 26,31). With only four free parameters, the model allows precise and unique parameter estimation from the force data. This conclusion was based on the convergence to the same parameter values starting from different initial estimates. In addition, this model is the first to be tested on muscles under different physiological conditions.
The process of muscular contraction at the molecular level is very complicated. Many steps, however, may not contribute significantly to force generation and would not be detectable from analysis of force data alone (27). Previous investigators have tried other approaches to a simpler model. Schultz and colleagues (25) developed a simple Hill-type model, also with four independent parameters. Their model was designed to predict the time course of the muscle force during isometric tetanic and isovelocity shortening contractions. Unlike the model presented in this study, however, the model developed by Schultz and colleagues cannot predict the force-time response to unfused tetanic, subtetanic, or variable-frequency stimulations. Recently, Shames and colleagues (26) developed a simple version of a Huxley-type model, also with four parameters, to explore the Ca2+-force relationship during isometric tetanic contractions. With use of the Ca2+ transient as input, their model successfully fit the forces produced during isometric tetanic contractions of various durations and isometric twitch contractions. However, this model has limited predictability, because the Ca2+ transient is not typically available.
An advantage of the present model over our previously reported model (31) is that it can predict force responses from fatigued muscles. The previous model failed to predict the forces from fatigued muscles by fitting the force in response to a single CFT. Subsequent studies with our previous model also showed that using the responses to more than one frequency of stimulation (e.g., the force produced in response to a high- and a low-frequency train, seemethods,Stimulations) to identify the parameter values of the model did not improve the accuracy of the model to predict forces. The present model was also successful with human muscles at long and short lengths, although the short-muscle-length condition was not tested on the previous model. The quality of the present model’s predictions was independent of the twitch contraction times of the muscles. In contrast, the previous model failed to predict forces of all six nonfatigued rat gastrocnemius muscles with twitch-contraction times >25 ms.
The advantage of Hill- and Huxley-type models over pure mathematical models is that the parameter values may be physiologically meaningful and the interpretation of the changes of those values with muscle condition can provide useful insight into the mechanisms of muscle contraction and fatigue (12). The parameters of this model, derived from Hill’s original model and based on the mechanisms of muscle contraction, may also be helpful in elucidating characteristics of muscle physiology.
Of the four model parameters, A is directly responsible for the magnitude of the force. When the muscle is fatigued or at short length, the force declines. Thus a decrease in the values of A for fatigued muscle or muscle placed at short length is expected. Our results are consistent with this expectation (Table 3).
The relaxation of the force is primarily determined by τ1, because the contribution of τ2 diminishes rapidly with the decline of CN , which always precedes the force relaxation (23). Fitch and McComas (14) observed a slower muscle relaxation in fatigued muscles, whereas Gandevia and McKenzie (15) observed faster relaxation of nonfatigued shortened human muscle. The increases and decreases in the fitted values of τ1 (Table 3) for fatigued and short-length muscles, respectively, are consistent with these previous findings.
The parameter τc is the characteristic time forCN to reach its peak and for it to decay to zero. The interpretation of the changes in τc across conditions is difficult. Many underlying physiological steps, such as Ca2+ release from and uptake by the sarcoplasmic reticulum, Ca2+binding to and dissociation from troponin, and cycling between the strongly bound and weakly bound cross bridges, may affect the time course of CN . Also, it is unclear which of the steps is rate limiting (27). Until we have a more complete understanding of cross-bridge activation, we cannot properly interpret τc. However, some previous experimental findings could, in part, explain our results. The prolonged Ca2+transient with fatigued muscles observed by Westerblad and Allen (30) could be one of the causes for an increased τc for fatigued muscles (Table3). Reduced Ca2+ sensitivity to troponin (3, 28) when muscles are at short length could cause less troponin-Ca2+ complex to be formed and, consequently, lead to activation of fewer cross bridges. Thus the decrease in τc with muscles at short length could also be explained.
As the characteristic time in the force equation (Eq.7 ), the sum of τ1 and some fraction of τ2 determines the overall rate of rise of force. When muscles were activated with CFTs, Binder-Macleod et al. (8) observed a slowing in the rate of rise of force as the muscles fatigued. In contrast, there was no difference in the rise time between the fatigued and the nonfatigued muscles when they were stimulated with VFTs (8). The increases of τ1 and τ2 for CFTs and increase of τ1 and decrease of τ2 for VFTs as the muscles become fatigued (Table 3) are, therefore, consistent with previous findings.
The major limitation of this model is that, for human muscles at short length and fatigued conditions, separate sets of parameter values were needed for CFTs and VFTs. In contrast, when human muscles were held at long lengths and not fatigued, only one set of parameter values was needed. However, the model did not always identify the stimulation pattern that produced the greatest force-time integral for this condition (Fig. 9). For this condition, CFT60 and VFT60 were identified as producing the greatest force-time integral. In contrast, the measured data showed that CFT90 and VFT100 produced the greatest force-time integrals among CFTs and VFTs, respectively. Similarly, the model consistently underestimated the forces at low frequencies (long IPIs) and overestimated the forces at high frequencies (short IPIs; Figs. 6 and 9). All these shortcomings could be due to some physiological steps not being taken into consideration, such as unequal amounts of Ca2+ released in response to each pulse of the train.
The mathematical model presented here is unique in its ability to predict isometric forces from rat and human skeletal muscles in nonfatigued and fatigued conditions at various muscle lengths. Modification of the present model are needed to allow a single set of parameter values to predict forces produced by CFTs and VFTs and to improve the model’s accuracy under various physiological conditions such as fatigued or nonfatigued muscle and short or long muscle length. In addition, the model will need to be modified to predict forces during repetitive stimulation of muscles that are allowed to shorten or lengthen if the model is to be used to predict forces during FES. By predicting the force that develops in response to an arbitrary pattern of stimulation, we envision mathematical models similar to the present one helping identify the optimal stimulation patterns for activation of skeletal muscle during FES.
The authors thank Samuel C. K. Lee, Lisa Landis, April Fritz, Lorin Kucharki, and Michelle Gerdom for assistance in data collection.
Address for reprint requests: S. A. Binder-Macleod, 315 McKinly Laboratory, University of Delaware, Newark, DE 19716.
This study was supported by National Institutes of Health Grant AR/HD-41254.
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. §1734 solely to indicate this fact.
- Copyright © 1998 the American Physiological Society
Muscle force has been assumed to be proportional to the amount of Ca2+-troponin complex in some models (11, 27, 31) or the amount of strongly bound cross bridges in others (26, 27, 29), suggesting that these factors play a similar role in muscle activation. Therefore, we used one unitless factor,C, to represent the single rate-limiting step before the myofilaments mechanically slide across each other. No experimental data have elucidated the dynamics of the Ca2+-troponin complex or the strongly bound cross-bridge transients during muscle contraction. However, studies of Chou and Hannaford (11) and Wexler et al. (31) suggest a Ca2+-transientlike curve for the dynamics of the Ca2+-troponin complex. Therefore, we assumed a similar transient shape forC, which can be expressed as Equation A1wheren is the number of pulses beforetime t in the train and has a maximum value of 6 in this study,ti represents the times that the electrical stimulation was applied, and and are the time constants controlling the shape and height ofC.
Pilot studies with this model showed that the force production was not sensitive to the magnitude of C, because the effect of increased height ofC could be compensated for by adjustment of parameter values in the force equation. Therefore, we normalized Eq. EA1 by dividing the right and left side by /( + ), which gives Equation A2 whereCN represents the normalized C and is equal toC( + )/ . Pilot studies showed a relationship between and and that the model worked best when was much larger than . Therefore,Eq. EA2 can be further simplified to Equation A3 There are a number of linearly dependent analytic solutions toEq. EA3 . The one chosen to represent the transient behavior ofCN is Equation A4 because, among the choices tested, this solution provided the best fits to the force data. Examples of the changes inCN are shown in Fig. 12.