Identification of a Continuous-Time Nonlinear Delay Dynamic System and Its Application to Neuromuscular Reflexes

CopYTight © IFAC System Identification, Copenhagen. Denmarlc. ) 994

IDENTIFICATION OF A CONTINUOUS-TIME NONLINEAR DELAY DYNAMIC SYSTEM AND ITS APPLICATION TO NEUROMUSCULAR REFLEXES LI-QUN ZHANG and W. ZEV RYMER Sensory Motor Performance Program. Rehabilitation Institute of Chicago Chicago, Illinois 60611. USA Abstract. The mechanical and reflex properties of the soleus muscle in decerebrate cats during white-noise random position perturbations were studied. The muscle's length and force were used to characterize the system dynamics of the muscle. A nonlinear dynamic system with internal delay was used to characterize the mechanical and reflex components independently. Continuous-time identification was used to estimate system parameters with genuine physical and physiological interpretations directly from the sampled input-output measurements. Results showed that the approach provided a useful way to characterize both the delayed (i.e. reflexive) and instantaneous mechanical contributions to net muscle force. Key words. Continuous-time identification; delays; neuromuscular reflexes; nonlinear systems

I. INTRODUCTION

instantaneous mechanical and delayed reflex actions. Continuous-time system identification was used to identify the model parameters from sampled data directly while retaining the parameters with their genuine physical or physiological interpretations. Results from both the soleus muscle of decerebrate cats and normal human elbow joint showed that the approach provided a helpful way to characterize the delayed reflex activation as well as the instantaneous mechanical properties of the nonlinear neuromuscular dynamics. Only the analysis on the single soleus muscle of decerebrate cats is presented here.

The role of muscle sensors in regulating muscle contraction has been studied intensively for several decades, yet the functional roles and efficacy of these regulatory systems remain unresolved due to difficulties encountered in studying them (Houk and Rymer, 1981). The difficulties have arisen for several different reasons. First, current models of neuromuscular dynamics generally do not differentiate muscle intrinsic and reflexive contributions adequately. Instead, they are usually lumped together in the models. Second, many quantitative descriptions of the components, such as muscle, muscle receptors, and spinal segmental neurons have relied heavily on linear systems approaches, which are largely inappropriate for the component characteristics. This is because muscle is a highly non-linear and state-dependent actuator with time-varying properties. Furthermore, many muscle receptors, especially those arising from muscle spindles are also highly nonlinear, making standard linear analyses inappropriate. Third, most analyses completely ignore the role of force sensors (particularly the Golgi tendon organs) and their central connections because these pathways have not been well quantified, making analysis and simulation difficult.

2. METHOD A nonlinear dynamic system with internal delays was used to study the nonlinear neuromuscular dynamics of the soleus muscle of decerebrate cats. First, nonlinear delay differential equations were used to characterize both the instantaneous mechanical and the delayed reflex actions of the muscles. The highly asymmetric property of the muscles was modelled non linearly . Second, continuous-time system identification was used to estimate the model parameters having physical or physiological interpretations directly from sampled data.

As a first step addressing these problems, a nonlinear continuous-time dynamic system with internal delay was used to characterize the neuromuscular dynamics (Zhang and Rymer, 1993). Nonlinear delay differential equations was used to differentiate the

2.1. Experimental Method To study both mechanical and neuro-sensory properties of a single muscle of decerebrate cats, 241

~~(s)

small amplitude random length perturbations were applied to the soleus muscle of decerebrate cats by a linear actuator under various background muscle activation levels modulated by the crossed-extension reflex. Muscle length, perturbing force, and EMG activity during the external perturbation were recorded. Muscle length and force were used to identify the neuromuscular dynamics directly at the mechanical level, and the EMG signal was used to check for reflex activation.

= {I + g(A )e- s,,,} / {sB(A) + s(rp(A)~ p(s)

N(s) +rn(A )~n (s»e- S1",

+ K(A) + K d(A )e-n", } (1)

.

.

B(A) ~(t) + K(A)~~(t) + rp(A )Rp(~(t - t dl ), Dp )

.

+rn(A)Rn(~(t - t dl ), Dn) + Kd(A )~~(t - t dl ) = N(t) + g(A )~f(t - t d2 ) + e(t)

The perturbation sequences were band-limited Gaussian, white noise, which impose random length perturbations around a designated muscle length. The bandwidth, defined as the range of frequencies within which 95% of the power resides, of the perturbation signal was from about 0-15Hz to 0-90Hz, depend on individual trials. If the reflex action needed to be emphasized for the purpose of the study, perturbations with a relatively narrower bandwidth was used. On the other hand, if the reflexive action needed to be minimized or even eliminated, perturbations with a very wide bandwidth were used.

(2)

~t)

Fig. I. A block diagram of neuromuscular dynamics of a single muscle.

Twelve normally distributed stochastic length perturbations with combinations of three peak-lopeak amplitudes (±O.3mm, ±1.5mm, and ±3.0mm) and three frequency bandwidths (0-1 5Hz, 0-35Hz, and 0-90Hz) were applied to the muscle via the length servo. Muscle length, force, and EMG signals were sampled throughout the trial at 333Hz after anti-alias filtering (100Hz cutoff). Background muscle force modulated by the crossed-extension reflex was maintained constant during each trial. If this background muscle force was not steady during the trial, only the steady segment of the data was analyzed. Different background muscle forces were used for different trials.

where A. is the operating state which includes the mean muscle length, background muscle force, average fusimotor activity, and perturbation bandwidth. s is the Laplace transfonn operator. 1;(t) is the muscle length and f(t) is the external perturbing force. ~1;(t) and M(t) are their deviations from their mean values during the small-amplitude random perturbations respectively. B(A.) is the instantaneous muscle mechanical viscous coefficient. K(A.) is the instantaneous muscle mechanical elastic stiffness. Rp(e) is a positive half-wave rectification function with a threshold Op and time delay id I . ~(s) is the corresponding Laplace transform operation. Rp(e) and the differential operator s block together describe the reflex activation evoked by muscle dynamic lengthening. idl corresponds to the delay from muscle length change (mainly through the muscle spindle activation) to the reflex contribution to muscle force. Time delays on the afferent and efferent pathway, and the neuromechanical delay were merged and represented by the single time delay id I . Op represents the threshold under which the muscle dynamic lengthening causes virtually no reflexive muscle force (not muscle spindle discharge which may still be invoked under the threshold Op). rp(A.) characterizes the corresponding strength of such reflexive activation. Because muscle spindles have highly asymmetric properties (lengthening vs. shortening) and muscle shortening may produce little or no response (Houk and Rymer, 1981), a similar tenn, rn(A.)Rn(e), is used to describe the dynamic reflexive activation caused by shonening the muscle. Kd(A.) describes the reflex activation caused by static muscle length deviation from its mean value. Therefore, the well-recognized velocity (dynamic) and position (static) sensitivities of stretch reflex are characterized by rp(A.) and rn(A.), and Kd(A.)

2.2 Nonlinear Delay Differential EQuations Muscle force is determined by both the mechanical and reflex properties of the muscle. The forces generated by the different components were lumped together and measured by the force transducer used in the experiment. A major difference in their activation mechanisms was used to differentiate the two kinds of components from the total measurement-the mechanical responses were instantaneous, while the reflex responses could not happen until both the afferent and efferent delays had elapsed. These delayed interactions fonned a component of our evaluation of the role of length and force sensors in regulating muscle mechanical behavior. Specifically, the block diagram in Fig. I was used to describe the neuromuscular dynamics involved. The corresponding system transfer function and the nonlinear delay differential equations are, respectively,

242

respectively. The block G(S)=g(A)e- st d2 describes the force feedback originated from the Golgi tendon organs with td2 being the overall time delay involved. g(A) represents the gain of force feedback or other more complex function. Because force feedback involves interneurons as well as the motoneurons and hence longer synaptic delays, ~2 is generally several ms longer than ~ I (Houk et al. 1981). Ms(t) and Mg(t) correspond to the reflexive forces modulated through muscle spindles and Golgi tendon organs respectively. e(t) is the modeling error.

'---.:,:.~--EJ-i

r-----'fOl I I I I 0 I I

I I

+

I I1

~l)

SP _

Motoneuron pool

Fig. 2. Block diagram of the single muscle neuromuscular dynamics showing the motoneuron pool. f(t) is the external penurbing force. E;(t) is the muscle length. M(t) and A~(t) are their deviations from their mean values during the smallamplitude penurbations respectively. Block T represents tendon organ feedback. P the spindle feedback. and M the muscle active contraction. Block C represents the external penurbation controller and actuator. SP represents the neurological feedback signal from muscle spindles. TO the feedback signal from Golgi tendon organs. and D the input to the motoneurons from other sources (the background excitatory level). AI is the operating state.

The block G(S)=g(A)e- st d2 in Fig. I (corresponding to the T*M blocks in Fig. 2) describes the force feedback but it does not form a direct feedback loop. This is because M(t) which is sensed by the Golgi tendon organs is the external perturbing force while the tendon organ feedback affects only the muscle contraction directly. Nonetheless the G block in Fig. I provides a quantitative measure of the Golgi tendon organ force feedback. Since length perturbation is applied to the muscle with muscle length being the controlled variable, ~ f(t) is determined by the perturbation controller (the C block in Fig. 2) based on ~~(t) which is affected by both the active and passive muscle forces. Therefore M(t) is indirectly affected by the Golgi tendon organ force feedback. Furthermore, with the system parameters identified, M(t) can be decomposed into the reflexive and mechanical forces which can provide further details of the force and length feedback mechanisms.

The motoneuron pool in Fig. 2 sums the neurological inputs SP, TO and D, and sends out the efferent signal E to the muscle. The active contraction of the M block produces force fa(t). The external perturbing force f(t) overcomes the active muscle contracting force fa(t) and passive muscle resistance (expressed as the muscle viscosity Band elastic stiffness K) to perturb the muscle length to ~(t). Since both the active and passive muscle forces are transmitted through the muscle tendons, the total force f(t) is sensed by the Golgi tendon organs and fed back to the motoneuron pool. Mathematically, the relationship can be described as: fp=f-fD=f-(D+P~-Tf)M

(3)

If small-amplitude perturbation around an operating state is used and assume that the background excitatory level D is constant during the smallamplitude perturbation, Eq. (3) is modified to

Compared to muscle length feedback, the muscle force feedback mechanism is much less well understood. It is not addressed in the first stage of the study, and the g(A) in Eqs. (I) and (2) was simply taken as zero for the study presented here. Further study on muscle force feedback is currently under investigation and will be reported elsewhere.

N p = !if -(P~~ -

T~f)M

(4)

Assuming the muscle contraction M block in Fig. 2 is linear around the operating state under the smallamplitude perturbations, superposition, the fundamental property of linear systems, can be applied and Eq. (4) becomes

2.3. Physiological Intewretation The block diagram in Fig. I was derived for the purpose of building a mathematical model characterizing both the muscle mechanical and reflex properties which can be identified from easily measurable data. It was not intended and was also not necessary to model the exact underlying physiological structures. For example, the upper summation operation of forces in Fig. I does not have a direct physiological correspondent. However, indirectly it characterizes the neurological signal summation at the motoneuron pool and the subsequent muscle contraction. To help comprehend the block diagram in Fig. I, the corresponding block diagram with the summation occurring at the motoneuron pool is shown in Fig. 2 (cf. Houk and Rymer, 1981, Houk et al. 1970).

!ifp

= !if -

PM~~

+ TM!if = N -

S~~

+ GN

(5)

Comparing Eq. (5) with Fig. I, it is readily seen that the above equation is the operation of the upper summation in Fig. I. The G and S in Eq. (5) correspond to the Golgi tendon organ force-feedback and spindle length-feedback respectively. Specifically, G=TM=g(A)e- st d2 and S= PM=(sR(s)+~)e-stdJ.The effects of the background excitatory level D are described by the operating state A in Fig. I. The reason for using the block diagram in Fig. I instead of that in Fig. 2 is that the former presents a model structure which is readily identifiable. In fact, assuming op and on, the force summations in Fig. I means the nonlinear model is linear in parameters-a parameter estimation problem which has been studied thoroughly (Golub and Van Loan, 1989; Ljung, 1987). Studying these 243

different force components as well as the parameters with genuine physical or physiological interpretations will help us get insight into the underlying system structure and mechanisms.

rectification was done according to the polarity checked. The threshold Op and on can be processed in a similar way by combining the threshold checking with the polarity checking operation. Since only the polarity of the velocity instead of the numerical values of the velocity itself was used, the problem associated with numerical approximation of differentiation was alleviated substantially. The problem was more serious for the threshold checking. However, the thresholds did not seem to be very significant factors.

2.4. Continuous-Time System Identification System identification has been done overwhelmingly in discrete time and continuous-time modeling was almost completely overshadowed until the past few years (Rao and Sinha, 1991). However, many physical models, like Eq. (2), are established in continuous-time and the model parameters have genuine physical or physiological interpretations. Generally, parameters of a linear continuous-time model can be estimated by the so-called "indirect" method-converting the differential equations into the corresponding discrete-time counterpartdifference equations, followed by estimating the discrete-time model and converting the discrete-time model parameters back to continuous time. However, difference equation is only an approximation of the original continuous-time differential equation because differentiation cannot be truly realized in digital environment (Rao and Sinha, 1991; Schoukens and Pintelon, 1991).

~jVY~~\J o

200

400

600

800

1000

1200

1400

1600

ms LonglIl (+veIoci1y)

~_~V~L~ 200

400

600

800 1000 ms Longlh (..,oIocily)

1200

1400

1800

~jr-=ru~\ ~ o 200

400

600

800 ms

1000

1200

1400

1600

velociry

j-~f\&\~\/ j o 200

400

600

800 ms

1000

1200

1400

1600

Fig. 3. Muscle length. its portions corresponding to positive and negative velocities respectively. and velocity.

The problem is especially important in our case where the major interest is the analysis of the continuous-time parameters with physical or physiological interpretations instead of using the identified system in a control application. More seriously, Eq. (2) is a nonlinear equation with items being a nonlinear function of past derivatives of muscle length. The nonlinearity makes it impractical to convert the continuous-time model into its discrete-time correspondent. The same problem exists for other "indirect" continuous-time system identification methods such as the frequency domain identification method (Schoukens and Pintelon, 1991) and the operator transformation method (Johansson, 1993).

The measured muscle length and force data over a period of about 6.2s were segmented into many pieces. Two consecutive segments overlapped each other by 50% of the segment length. Model parameters were estimated over each segment and trimmed ensemble averages were calculated for each of the parameters studied. The parameter estimations presented in the next section were all obtained in this way.

3. RESULTS Since muscle mechanical properties have been discussed by many investigators (Houk and Rymer, 198 1; Kearney and Hunter, 1990), only the results related to the reflex contributions are discussed here. The values of the reflex delay Idl were initially taken from literature (Houk et aI., 1981) and then tested in the neighborhood. The one which gave the minimum simulation error was selected. As a result, 36ms was used for the following results. Op and on were taken as zero.

The so-called "direct" methods was used to identify the above nonlinear continuous-time system from sampled data while retaining the system parameters with their physical or physiological interpretations (Rao and Sinha, 1991). Specifically, both sides of the system dynamic equations such as Eq. (3) were integrated over various sampling intervals using linear integral filtering and finite impulse response filtering (FIR) (Sagara and Zhao, 1990; Zhao et aI., 1991).

Both the model structure with the Kct(A.) term in Eq. (2) and the one without it were studied. Fig. 4 shows the results when both K(A.) and Kct(A.) were modelled. Each point in the figure was the ensemble average of II estimations. Since the background muscle force modulated through the crossed-extension reflex may not be very stable during a trial, there may be some "outliers" in the parameter estimates. To reduce the effects of such "outliers" and obtain robust estimates, trimmed means were calculated by

Nonlinearity of the rectification operation can be dealt with by combining the integration with the rectification. First, the muscle length was divided into portions with positive and negative velocities respectively, as shown in Fig. 3. Second, the length signal was low-pass filtered and numerically differentiated. Third, the polarity of the past velocity d~(t-td I )/dt was checked, then the integration corresponding to either the positive or negative

244

excluding the highest and lowest 5% of the estimates. Each curve in Fig. 4 was obtained by fitting the individual K(A.) or 1
1.5

As shown in Fig. 4, when the 1
., 1imo1MC)

Fig. 5. Simulation results from the estimated system parameters. The solid line represents the measured muscle length. The dotted line corresponds to the simulation result based on measured muscle force.

Fig. 6 shows the estimated parameters related to the dynamic stretch reflex. As shown by the three upper curves in Fig. 6, for the soleus muscle of a decerebrate cal under Gaussian random perturbation of 15Hz bandwidth, mechanical viscosity B (the solid lines) was generally larger than the reflex parameters rp (the dashed lines) and rn (the dotted lines), which is actually a stability condition of the system when rp and rn are equal (Kuang 1993). rp was generally higher than rn. which means stretching the muscle generated stronger reflexive activation than did the shortening muscle. This is true for both the models with the 1
1500

E 1000

Z

500

o

•

+

Kd

.---~

5 6 Force (N)

9

10

Fig. 4. The mechanical elastic stiffness K(A) (the solid curve) and reflex static length feedback Kd(A) (the dashed curve,) contributions of the soleus muscles, shown as a function of the background muscle force (the abscissa). The data were obtained under the perturbation of tl.6mm and 15Hz bandwidth.

The insignificance of the 1
8O'r---~--~--~--~----,

B 70 60 50

.

_. lP

B

30

Fig. 5 shows the comparison between the measured muscle length (the solid line) and the simulated length (the dotted line) obtained with the estimated model parameters (K=1951 N/m, B=16.0 N-s/m, rp=3.93 N-slm and rn=2.14 N-slm). The simulation was done with known input f(t) (muscle force) and zero initial condition. Measured muscle length information 1;(t) was not used in the simulation.

~--

.0

m

20

10

.. ''0'

.,

lP

_Cj_o~ + -:.~..-;:;~ m °0~L§JE:.II;~s:..:.;:'::='~-==~6:'::""::--:'::""-~--.J'0 Force (N)

Fig. 6. The reflex contributions of a cat soleus muscle shown as a function of the background muscle force (the abscissa). The upper three curves are the parameter B (* and the solid line), rp (+ and the dashed line) and r n (0 and the dotted line) under the

245

penurbation of ± 1.6mm peak-to-peak amplitude and bandwidth of 15Hz. The lower three curves are the B. rp and r n parameters for the penurbation of ±1.6mm and 90Hz bandwidth. Each curve was obtained by fining the individual B. r p or r n values obtained under different background muscle forces with a second-order polynomial.

Kuang, Y. (1993). Delay Differential Equations, Academic Press. Ljung, L. (1987). System Identification: Theory for the User, Prentice-Hall, New Jersey. Rao, G. P. and N. K. Sinha (1991). "ContinuousTime Models and Approaches", in N. K. Sinha and G. P. Rao (ed.), Identification of ContinuousTime Systems - Methodology and Computer Implementation, pp. 1-15, Kluwer Academic, Boston. Sagara, S. and Z. Zhao (1990). "Numerical Integration Approach to On-line Identification of Continuous-time Systems," Automatica, Vol. 26, pp. 63-74.

4. COLCLUSION The nonlinear delay differential equations presented in this study provide us a useful tool to characterize the muscle mechanical and reflexive activations simultaneously and independently. The nonlinear delay differential equations also characterize the highly asymmetric muscle properties in a nonlinear way. Several variations of the model structures were discussed. These model structures are not overly complex and hence provide us a tractable method to characterize muscle properties reliably.

Schoukens, J. and R. Pintelon (1991). Identification of Linear Systems, Pergamon Press, Oxford. Zhang, L.-Q. and W. Z. Rymer (1993). "Modeling of Muscle Mechanical and Reflexive Properties Using Nonlinear Delay Differential Equations", Proc. 15th Ann. Int. Con! IEEE Eng. Med. BioI., San Diego, pp. 1165-1166.

Further studies can be done to model force feedback as well as length (velocity) feedback simultaneously. Furthermore, combinations of these tractable models over different operating ranges will provide us more comprehensi ve method to address the very complex neuromuscular dynamics.

Zhao, Z.-Y., S. Sagara, and K. Wada (1991). "BiasCompensating Least Squares Methods for Identification of Continuous-Time Systems From Sampled Data", Int. 1. Control, Vol. 53, pp. 445-461.

ACKNOWLEDGMENT The authors would like to thank R. Kirsch and D. Boskov for collecting the cat data, and acknowledge the support of NIH and the Whitaker Foundation. REFERENCES Golub, G. H. and C. F. Van Loan (1989). Matrix Computations, 2nd ed., Johns Hopkins Univ. Press. MD. Houk, J. C. and W. Z Rymer (1981). "Neural Control of Muscle Length and Tension", in V. B. Brooks (ed.), Handbook of Physiology, Chap. 8, pp. 257-323. Houk, J. c., P. E. Crago and W. Z Rymer (1981). "Function of Spindle Dynamic Response in Stiffness Regulation - a Prediction Mechanism Provided by Nonlinear Feedback", in A. Taylor and A. Prochazka (ed.), Muscle Receptors and Movement, Macmillan, London, pp. 299-309. Houk, J. c., J. J. Singer, and M. R. Goldman (1970), "An Evaluation of Length and Force Feedback to Soleus Muscles of Decerebrate Cats", J. Neurophysiol., Vol. 33, pp. 784-811. Johansson, R. (1993). System Modeling and Identification, Prentice Hall, NJ. Kearney, R. E. and I. W. Hunter (1990). "System Identification of Human Joint Dynamics", CRC Critical Review in Biomed. Eng., Vol. 18, pp. 55-87. Kearney, R. E. and I. W. Hunter (1988). "Nonlinear Identification of Stretch Reflex Dynamics", Ann. Biomed. Eng., Vol. 16, pp. 79-94. 246