Mathematical Model That Predicts Lower Leg Motion In Response To Electrical Stimulation

 

Perumal R 1, Wexler AS 2, Binder-Macleod SA 3

 

1 Department of Mechanical Engineering, University of Delaware, Newark, DE-19716, USA.

2 Mechanical and Aeronautical Engineering, One Shields Ave., UC-Davis, CA-95616, USA.

3 301 McKinly Lab, Department of Physical Therapy, University of Delaware, Newark, DE-19716, USA.

 

Email: perumal@me.udel.edu


Abstract

Direct electrical activation of skeletal muscles of patients with upper motor neuron lesions can restore functional movements, such as standing or walking. Because responses to electrical stimulation are highly nonlinear and time varying, accurate control of muscles to produce functional movements is very difficult. Accurate and predictive mathematical models can facilitate the design of stimulation patterns and control strategies that will produce the desired force and motion. The purpose of this study was to validate our previously developed model during nonisometric contractions when the leg was allowed to move freely in response to electrical stimulation. Our results showed that the model could accurately predict the angular position and velocity of the lower leg when the muscle was stimulated with wide range of clinically relevant stimulation patterns and with different loads placed around the ankle joint.

1        Introduction

Functional Electrical Stimulation (FES) is the coordinated electrical excitation of paralysed or weak muscles in patients with upper motor neuron lesions to produce functional movements such as sit-to-stand or walking. Despite many technical advances, FES has had a relatively minor impact on rehabilitation because the physiological and biomechanical processes involved in the generation of FES-elicited movements are highly non-linear and time varying. Hence, it is difficult to determine the appropriate stimulation patterns necessary to produce the desired muscle force and limb motion. However, the use of mathematical muscle models can improve and accelerate the development of FES systems for practical use.

Although a number of mathematical models capture the behaviour of muscles [1][2], to date none of the models have succeeded in FES systems, because they lacked the ability to predict muscle forces, moments, or motions in response to a wide range of stimulation frequencies, patterns, and physiological conditions. Previously, we developed a Hill-type model that predicted the force responses of the quadriceps femoris muscle to a wide range of clinically relevant stimulation patterns at different muscle lengths and velocities [3][4]. The purpose of this study was to validate our model during nonisometric contractions when the lower leg was allowed to move freely in response to electrical stimulation. Validation was done by comparing the model predictions to the experimental data obtained during stimulation of the quadriceps femoris muscle with a range of stimulation patterns and loads placed around the ankle joint.

2           Methods

The model has three differential equations:             ,   (1)

where  for i = 1, and   for i >1;                                                                                           

             (2)

where, and ;     

and for general nonisometric movements, hence

,       (3a)

and for isovelocity movements ,                                                                                        (3b)     

 .                               (3b)

Equation (1) represents the dynamics of the rate-limiting step leading to the formation of the Ca2+-troponin complex, CN. t (ms) is the time since the beginning of stimulation, ti (ms) is the time when the ith pulse is delivered, τc is the time constant controlling the rise and decay of CN, and R0 characterizes the magnitude of enhancement in CN when there are two closely spaced pulses.  Equation (2) represents the development of the force due to stimulation, F (N) . A40 (N/ms) is the scaling factor for force at 40° of knee flexion, a (N/ms-deg2) and b (N/ms-deg) are scaling factors to account for force at each knee flexion angle, q (deg) is the knee flexion angle,  (deg/s) is the angular velocity of the shank, V1 (N/deg2) is a constant for each subject to be determined at -200°/s, Km is the sensitivity of strongly bound cross-bridges to CN, t1 (ms) is the time constant of force decline in the absence of strongly bound cross-bridges, and t2 (ms) is the time constant of force decline when there are strongly bound cross-bridges. Equation (3) is the equation of motion for the shank in response to electrical stimulation [4]. Two separate equations (3a and 3b) were used to represent the motion of the shank during general nonisometric and isovelocity movements.(deg/s2) is the acceleration of the shank, L is the distance from the knee joint centre of rotation to the location of F above the ankle joint, I (kg-m2) is the mass moment of inertia of the lower leg and the external load wrapped above the ankle joint, M (N) is the resistance to knee extension that includes the visco-elastic resistance of the knee joint and the weight of the lower leg,  FKC (N) is the experimental force measured by the KinCom dynamometer during isometric and isovelocity experiments, FLOAD (N) is the external load wrapped above the ankle joint, and (deg) is the resting angle.

Six healthy subjects were recruited for this study. Isometric, constant shortening velocity (isovelocity) force, and varying velocity (nonisometric) angle data were collected of the human quadriceps femoris muscle in response to electrical stimulation of the human quadriceps femoris muscle. For the isometric and isovelocity contractions, force and angle data were collected using the KinCom. In contrast, for the nonisometric contractions the dynamometer arm was replaced by a low frictional resistance custom built arm that measured knee joint angle.

 

Figure 1: Schematic representation of the three stimulation patterns tested. Bottom train (CFT50) is a constant-frequency train with all interpulse intervals equal to 50 ms; middle train (VFT50) is a variable-frequency train with an initial doublet of 5 ms and remaining pulses equally spaced by 50 ms; and top train (DFT50) is a doublet-frequency train with 5-ms doublets separated by interdoublet interval of 50 ms. The train’s name is based on the duration of the longest interpulse interval within that train. Each train has a maximum duration of 1 sec.

For each subject all data were collected in a single testing session. Each testing session consisted of three parts. First, subjects were tested isometrically at knee flexion angles of 15°, 40°, 65°, and 90°. Then, subjects were tested at an isovelocity speed of -200°/s (all shortening velocities are assigned negative values in this study). Finally, nonisometric tests were performed with weights of 0, 4.54 and 9.08 kg wrapped just above the ankle joint. A rest period of 5 min was provided between testing each angle, velocity, and load to avoid fatigue. The stimulation intensity was set to activate  ~40% of the muscle. For isometric and isovelocity testing the muscle was stimulated with one-second long VFT20 and VFT80 (Figure. 1), while under nonisometric conditions the muscle was stimulated with one-second long CFTs, VFTs, and DFTs of interpulse intervals (IPIs) ranging from 10ms to 100ms. The experimental procedures for isometric and isovelocity studies are detailed in [4]. For nonisometric testing, stimulations were initiated when the subject’s leg was in a resting position and were truncated when either the leg had gone through a 45° excursion or all pulses within the one-second long train were delivered. The above trains were delivered every 10 sec to avoid fatigue.

First, parameters a, b, A40, τ1, τ2, and KM were identified under isometric conditions, then parameters V1 and M were identified at a constant shortening velocity of -200°/s, and finally parameters , M, and  were identified under nonisometric conditions. Under isometric and isovelocity conditions parameters were identified by fitting the force response to the VFT20-VFT80 train combination. For nonisometric contractions parameters were identified by fitting the velocity response to VFT20 train at each load condition.


The accuracy of the predictions of the response to the nonisometric condition were tested by comparing measured and modeled knee flexion angles and velocities in response to the CFTs, VFTs, and DFTs for all the six subjects tested. We used the RMS angle error measure to see how well the model predicted leg position in response the stimulation at different loads.  The RMS angle error was defined as


                               (4)

where qmod,i is the modeled knee flexion angle, qmeas,i is the measured knee flexion, i is the increment time interval of 5 ms, and N is the total number of data points from the beginning of the stimulation to the measured knee flexion angle at the end of extension.

 

3        Results

Subjective evaluations showed that the model was able to predict the knee flexion angle and the angular velocity in response to different stimulation IPIs and patterns (Figure 2). The averaged angle errors were generally lower at the two higher loads and were less than 8° for most of the stimulation IPIs tested (Figure 3). In addition, the averaged RMS angle errors were generally higher at longer IPIs for each of the two loads tested.

 

4        Discussion and Conclusions

The model developed and tested here accurately predicted the angular position and velocity of the shank when the muscle was stimulated with wide range of stimulation frequencies and patterns, and loads. Compared to other models [1][2], the current model requires very few parameters to describe the nonisometric motion of the limb in response to electrical stimulation. This makes the current model a suitable candidate for adaptive control algorithms that can track parameter variations online and provide a better control of motion during FES.

 

 

Figure 2: Predicted and measured knee flexion angles and shortening angular velocities during extension from a typical subject at load of 9.08 kg for four of the 16 stimulation trains.

Figure 3: Bar graphs showing the averaged angle errors (+ standard error) for six subjects in response to the 16 stimulation trains at 0 and 9.08 kg different loads.

References

[1]         Riener R, Quintern J, Schmidt G. Biomechanical model of the human knee evaluated by neuromuscular stimulation. J. Biomech. 29, 1157-1167, 1996.

[2]         Dorgan SJ, O’Malley MJ. A mathematical model for skeletal muscle activated by N-let pulse  trains. IEEE Trans. Rehab. Eng. 6, 286-299, 1998.

[3]         Perumal R, Wexler AS, Ding J, et al. Modeling the length dependence of isometric force in human quadriceps muscles. J. Biomech. 35, 919-930, 2002.

[4]         Perumal R, Wexler AS, Ding J, et al. Mathematical Model For Electrically Elicited Isovelocity Human Muscle Contractions —I: Development Of The Model. ASME J.Biomech. Eng. In Review.   

Acknowledgements

This study was supported by the National Institutes Health Grant HD 36797.