JN Fuel your research with LabChart
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


J Neurophysiol 96: 1772-1788, 2006. First published May 3, 2006; doi:10.1152/jn.00868.2005
0022-3077/06 $8.00
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow All Versions of this Article:
96/4/1772    most recent
00868.2005v1
Right arrow Alert me when this article is cited
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 ISI Web of Science (6)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Mileusnic, M. P.
Right arrow Articles by Loeb, G. E.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Mileusnic, M. P.
Right arrow Articles by Loeb, G. E.

Mathematical Models of Proprioceptors. I. Control and Transduction in the Muscle Spindle

Milana P. Mileusnic1, Ian E. Brown3, Ning Lan2 and Gerald E. Loeb1

1Department of Biomedical Engineering, Alfred E. Mann Institute for Biomedical Engineering and 2Department of Biokinesiology and Physical Therapy, University of Southern California, Los Angeles, California; and 3Center for Neuroscience Studies, Queen's University, Kingston, Ontario, Canada

Submitted 18 August 2005; accepted in final form 1 March 2006


    ABSTRACT
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX 1
 ACKNOWLEDGMENTS
 REFERENCES
 
We constructed a physiologically realistic model of a lower-limb, mammalian muscle spindle composed of mathematical elements closely related to the anatomical components found in the biological spindle. The spindle model incorporates three nonlinear intrafusal fiber models (bag1, bag2, and chain) that contribute variously to action potential generation of primary and secondary afferents. A single set of model parameters was optimized on a number of data sets collected from feline soleus muscle, accounting accurately for afferent activity during a variety of ramp, triangular, and sinusoidal stretches. We also incorporated the different temporal properties of fusimotor activation as observed in the twitchlike chain fibers versus the toniclike bag fibers. The model captures the spindle's behavior both in the absence of fusimotor stimulation and during activation of static or dynamic fusimotor efferents. In the case of simultaneous static and dynamic fusimotor efferent stimulation, we demonstrated the importance of including the experimentally observed effect of partial occlusion. The model was validated against data that originated from the cat's medial gastrocnemius muscle and were different from the data used for the parameter determination purposes. The validation record included recently published experiments in which fusimotor efferent and spindle afferent activities were recorded simultaneously during decerebrate locomotion in the cat. This model will be useful in understanding the role of the muscle spindle and its fusimotor control during both natural and pathological motor behavior.


    INTRODUCTION
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX 1
 ACKNOWLEDGMENTS
 REFERENCES
 
The muscle spindle is a sense organ found in most vertebrate skeletal muscles. In a typical mammalian lower limb muscle, several tens (or even hundreds) of muscle spindles can be found lying in parallel with extrafusal fibers and experiencing length changes representative of muscle length changes (Barker 1962Go; Eldred et al. 1962Go). The spindle has been found to play a dominant role both in kinesthesia and in reflexive adjustments to perturbations (Hulliger 1984Go; Matthews 1981Go). Its sensory transducers (primary and secondary afferents) provide the CNS with information about the length and velocity of the muscle in which the spindle is embedded. The spindle provides the main source of proprioceptive feedback for spinal sensorimotor regulation and servocontrol. At the same time that the spindle supplies the CNS with afferent information, it also receives continuous control through specialized fusimotor efferents (static and dynamic fusimotor efferents) whose task is to shift the spindle's relative sensitivities over the wide range of lengths and velocities that occur in various natural tasks (Banks 1994Go; Matthews 1962Go).

The spindle consists of three types of intrafusal muscle fibers: long nuclear bag1 and bag2 fibers and shorter chain fibers (Fig. 1A) (Boyd et al. 1977Go). Typically, one bag1, one bag2, and about four to 11 chain fibers lie in parallel within a spindle (Boyd and Smith 1984Go). The bag1 fiber is the only fiber that has dynamic fusimotor efferent endings located on it and is primarily responsible for the velocity sensitivity of the spindle. The bag2 and chain fibers receive static fusimotor control and contribute mainly to length sensitivity. Located in the equatorial region of all three types of intrafusal fibers are primary afferent endings responsible for supplying the CNS with the length and velocity information from the muscle. The more slowly conducting secondary afferent has its endings located more eccentrically, on only the bag2 and chain fibers, and is responsible for supplying CNS with primarily length information.


Figure 1
View larger version (37K):
[in this window]
[in a new window]
 
FIG. 1. A biological muscle spindle and the structure of the spindle model. A: a biological muscle spindle consists of 3 types of intrafusal fibers that receive several fusimotor inputs (gamma static and dynamic) while giving rise to primary (Ia) and secondary (II) afferent (modified from Bakker 1980Go). B: spindle model consists of 3 intrafusal fiber models; it receives 3 inputs (fascicle length, in terms of optimal length L0, and static and dynamic fusimotor drives) to produce primary and secondary afferent firing.

 
Because of the difficulty of accurately recording afferent and especially fusimotor activity during motor behavior, theories of motor control usually rely on assumptions about spindle activity. Throughout the years, several attempts have been made to formalize these assumptions in mathematical models capable of accurately capturing spindle activity over the wide range of kinematic and fusimotor conditions in which spindles naturally operate. These modeling approaches involved either transfer functions (Chen and Poppele 1978Go; Matthews and Stein 1969Go), nonlinear functions based on curve-fitting (Houk et al. 1981Go; Maltenfort and Burke 2003Go), or reduction to constituent anatomical components (Hasan 1983Go; Lin and Crago 2002Go; Rudjord 1970Go; Schaafsma et al. 1991Go). As a consequence of the complex nature of spindle responses, the earliest models attempted to model primarily afferent activity in the absence of fusimotor activation, with the exception of Hasan (1983)Go. Only recently have several more complex spindle models, incorporating the fusimotor effects on afferent activity, been developed (Lin and Crago 2002Go; Maltenfort and Burke 2001; Schaafsma et al. 1991Go). Although each new model has enhanced our understanding of how the spindle operates, there are still certain aspects of spindle behavior that these models fail to capture (such as occlusion between transduction regions and temporal dynamics of changing fusimotor input; see DISCUSSION).

We constructed a physiologically realistic model of the spindle that is composed of mathematical elements closely related to the anatomical components found in the biological spindle and their physiological properties. Its parameters were optimized using cat soleus muscle afferent data recorded during a variety of kinematic and fusimotor conditions. The emergent behavior of the model is shown to account well for the complete range of complex phenomena described in various experiments. The model was validated against recently reported data on medial gastrocnemius afferent activity during decerebrate locomotion of the cat that include direct recordings of related fusimotor efferent activity. Preliminary reports of this study were previously published (Mileusnic et al. 2001Go, 2002Go).


    METHODS
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX 1
 ACKNOWLEDGMENTS
 REFERENCES
 
In this section we will initially concentrate on the equations and individual terms that describe the spindle model and how they are related to the current understanding of the biological system. Afterward, the implementation of the model will be discussed and will include the criteria for selecting the afferent records and techniques used to perform the parameter optimization.

The spindle model

The spindle model is composed of three intrafusal fiber models (nuclear bag1 and bag2 fibers and chain fiber) and two afferent firing summation models (primary and secondary afferent firing models) (see Fig. 1). The same basic function was used to model all three intrafusal fiber types, with different coefficient values to account for their different physiology and effects. Each intrafusal fiber model responds to two inputs: fascicle length (L; in units of L0, which represents the optimal muscle fascicle length) and the relevant fusimotor drive [in the case of bag1 fiber it is dynamic fusimotor drive ({gamma}dynamic), whereas in the case of bag2 it is chain static fusimotor drive ({gamma}static) (Fig. 1B)]. The spindle model generates two outputs: primary (Ia) and secondary (II) afferent activity. The primary afferent response results from the contributions of all three intrafusal fiber models on which the primary afferent receptor has transduction endings. The secondary afferent receives inputs from only the bag2 and chain intrafusal fiber models.

In the next two major sections we will develop the spindle model in stages. First we will describe the general intrafusal fiber model. This section of model development is further subdivided into the three components that deal with fusimotor activation, the mechanics of stretch within the intrafusal fiber, and finally sensory transduction from stretch to afferent endings. Second, we describe the afferent firing model that deals with nonlinear summation between the intrafusal fibers' transduction regions.

The intrafusal fiber model

All the intrafusal fiber models share the same general structure, a modified version of McMahon's spindle model (McMahon 1984Go), which has its origins in the earlier work by Crowe (1968Go, 1970Go). The relative importance of model parameters differs for three intrafusal fiber models to account for the different properties of three fiber types (see Table 1 and Fig. 2). The form of the intrafusal fiber model was also influenced by recent improvements in the modeling of extrafusal mechanics (Brown and Loeb 2000Go; Brown et al. 1996bGo, 1999Go).


View this table:
[in this window]
[in a new window]
 
TABLE 1. Spindle model parameters

 

Figure 2
View larger version (8K):
[in this window]
[in a new window]
 
FIG. 2. Intrafusal fiber model. All intrafusal fibers consist of polar and sensory zones with qualitatively identical mechanical components. For each fiber model the stretch in the sensory and polar regions is calculated to determine its contribution to the firing of each afferent. Model parameters for each intrafusal fiber type are provided in Table 1.

 
BACKGROUND. The linear spindle model suggested by McMahon (1984)Go makes no distinction between different intrafusal fibers, but rather the fibers are lumped together into one spindle model. Similar to the structure of intrafusal fibers found in the biological spindle, McMahon's spindle model consists of two regions, sensory and polar (a modified version is shown in Fig. 2). The sensory (transduction) region contains afferent endings wrapped around it. Stretch of this region results in distortion of afferent endings, depolarization of their membranes, and increase in the rate of action potential firing (Boyd and Smith 1984Go). McMahon modeled transduction as a spring whose stretch is proportional to afferent firing. The rest of the intrafusal fiber on either side of the sensory region is called the polar region, constituting essentially a striated muscle fiber innervated by fusimotor endings; for simplicity, the two polar regions of the biological spindle are combined into one polar region for the model. McMahon used a typical Hill model for the polar region, including a passive spring in parallel with a contractile element, which further consists of an active force generator and a damping element. The contractile element was designed to represent the properties of the spindle that change in response to fusimotor activation, although these coefficients were never explicitly modeled by McMahon.

FUSIMOTOR ACTIVATION. A biochemical Hill-type equation was used to convert the actual fusimotor frequency [{gamma}dynamic or {gamma}static in pulses per second (pps)] to an equivalent activation level (fdynamic or fstatic, defined within range 0 to 1). This allows each intrafusal fiber model to capture the saturation effects that take place at high fusimotor stimulation frequencies as observed in isolated, identified fibers (bag1: freqbag1 = 100 pps; bag2: freqbag2 = 100 pps; chain: freqchain = 150 pps; Boyd 1976Go). In addition, the model incorporates the different temporal properties of intrafusal fiber responses that were measured previously for individual intrafusal fibers in response to step changes in fusimotor activation (Boyd et al. 1977Go). These different temporal properties are thought to arise from differences in the spread of activation in twitch muscle fibers that propagate action potentials (including chain fibers) versus tonic muscle fibers where synaptic depolarization spreads electrotonically (including bag fibers; Boyd 1976Go), as well as being related to differences in calcium kinematics. To model these differences, low-pass filters between the fusimotor inputs and the equivalent activation levels were introduced for the two relatively slow bag intrafusal fibers (see APPENDIX 1, A). The following equations were used to convert the actual fusimotor frequency ({gamma}dynamic or {gamma}static) to an equivalent activation level (fdynamic or fstatic) for three types of intrafusal fibers

Similar to McMahon's spindle model, our final model incorporates a sensory region of each intrafusal fiber model as a pure elastic element with a spring constant KSR. The tension within this region is

Formula 1(1)

MECHANICS OF SENSORY AND POLAR REGIONS. With the intention of keeping our model as simple as possible, we initially assumed our intrafusal fiber model to have McMahon's intrafusal fiber model structure. By manually tuning all the parameters in the model, the model's ability to accurately account for experimentally observed spindle afferent activity during a variety of kinematic conditions and fusimotor activations was assessed and the potential need for additional terms was identified. Addition of every further term involved extensive testing of the model's performance to reduce the danger of overcomplicating the model with terms whose functions were unknown.

equal to

Formula 2(2)
where L is fascicle length, LPR is polar region length, and L0SR is the unloaded sensory region length, all in units of L0.

The polar region is modeled as a spring with a spring constant KPR and a parallel contractile element that consists of an active force generator and a damping element, both of whose properties are modulated by fusimotor input. The tension within this region is equal to

Formula 3(3)
where

Formula 4(4)
and

Formula 5(5)
M refers to the intrafusal fiber mass required for computational stability in a series-elastic system with velocity-dependent contractility (Brown et al. 1996bGo). L0PR is a polar region's rest length. beta represents the polar region's damping term; increases in beta result in increases in the velocity sensitivity of the primary afferent, which plays an important role in modulating the spindle's behavior during dynamic fusimotor stimulation of bag1 fiber (beta1) (the only intrafusal fiber receiving dynamic fusimotor drive) (Crowe and Matthews 1964Go). By contrast, static fusimotor activation produces a small decrease in beta (note that beta2 is negative) (Crowe and Matthews 1964Go). {Gamma} is defined as the active-state force generator term; increases in {Gamma} result in an increase in the bias activity of the dependent afferent. Static fusimotor input causes a sustained, strong contraction within the bag2 and chain polar regions ({Gamma}2 x fstatic), producing a stretch in the sensory region and a bias in the afferent activity. Dynamic fusimotor input produces a similar but much weaker effect ({Gamma}2 x fdynamic). The model incorporates the experimentally observed nonlinear velocity dependency of the spindle's afferent response and was modeled with the velocity power term (LFormula 5R)a (Houk et al. 1981Go; Prochazka and Gorassini 1998Go). C is a constant describing the experimentally observed asymmetric effect of velocity on force production during lengthening and shortening. Although this asymmetry has been observed directly only in extrafusal striated muscle (Scott et al. 1996Go), we have assumed its existence also in the case of the intrafusal fiber's contractile polar region. C was set to unity during polar region lengthening (C = CL = 1) and to CS during shortening (C = CS). Finally, the model incorporates the length dependency of the force–velocity relationship (term LPR – R), where an increase in fascicle length results in increased slope of the force–velocity relationship for slow to moderate velocities. This effect, observed in extrafusal fibers, is believed to result from the influence of myofilament lattice spacing on cross-bridge kinetics (Brown et al. 1999Go). Under most physiological conditions the sarcomere length of the intrafusal fiber's polar region tends to follow that of extrafusal fibers (extrafusal sarcomere: Scott et al. 1996Go; intrafusal sarcomere: Barker 1974Go; Bessou et al. 1975; Boyd 1976Go; Poppele and Quick 1981Go), so the extrafusal fiber measurements were used to estimate the length (R) of the polar region of intrafusal fibers for this effect (assuming that length changes in sensory region are minor comparing to those in the polar region; see Implementation of the spindle model and parameter determination).

INTRAFUSAL FIBER CONTRIBUTIONS TO AFFERENT FIRING. In our spindle model, the stretches in the sensory and polar regions are independently calculated for each intrafusal fiber to determine their potential contributions to afferent firing. Because tensions in Eqs. 2 and 3 are the same, the sensory region's equation for tension (Eq. 2) was rearranged to express polar region length (LPR) in terms of tension (T) and fascicle length (L). This polar region length was then substituted into Eq. 3 to obtain a second-order differential equation oftension in terms of fascicle length

Formula 6(6)
Because the primary afferent endings are located on the sensory regions of all three intrafusal fibers, the stretch in the sensory region of each intrafusal fiber is calculated (T/KSR). Once the afferent endings are stretched passed a certain sensory region length (sensory region threshold length: LNSR) the ion channels open and depolarization/impulse generation takes place. The stretch above the threshold length is scaled by a constant G (the term that relates the stretch of the intrafusal fiber's sensory region to primary afferent firing; see Implementation of the model and parameter determination) to obtain the intrafusal fiber's contribution to the activity of the primary afferent (before nonlinear intrafusal fiber firing summation between bag1 and combined bag2 and chain fiber models; for more details see next section)

Formula 7(7)
Contrary to the primary afferent endings, the secondary afferent transduction zones are located more eccentrically, straddling both sensory and polar regions of bag2 and chain intrafusal fibers. Therefore action potential generation reflects the stretch in both sensory and polar regions

Formula 8(8)
X represents the percentage of the secondary afferent located on the sensory region, which can vary among spindles. The stretch within the part of the secondary afferent that is located on the sensory region is obtained by multiplying X by sensory region stretch above the sensory region threshold length (LNSR) and normalizing it by the ratio of the secondary afferent rest length (Lsecondary) and sensory region rest length (L0SR). The stretch within the part of the secondary afferent located on the polar region is similarly obtained. The polar region length at and above which secondary afferent sensory endings are stretched is defined as the polar region threshold length (LNPR). Once the stretches of the secondary afferent portions that are located on both sensory and polar regions are obtained, they are summed together and multiplied by G to obtain the intrafusal fiber contribution to secondary afferent firing.

Afferent firing model

The output of the primary afferent is captured by nonlinear summation between the bag1 and combined bag2 plus chain intrafusal fiber outputs, to account for the effect of partial occlusion that has been observed in primary afferents during simultaneous static and dynamic fusimotor stimulation (Banks et al. 1997Go; Carr et al. 1998Go; Fallon et al. 2001Go). Such combined fusimotor stimulation produces an afferent response that is greater than the larger of the individual responses (during either static or dynamic fusimotor stimulation) but smaller than their sum. Although the mechanism responsible for the occlusion is still debated (Banks et al. 1997Go; Carr et al. 1998Go; Fallon et al. 2001Go; Proske and Gregory 2002Go), one likely explanation assumes the existence of two impulse generators in the spindle afferent. One impulse generator is believed to be located on the bag1 fiber, whereas the second combines generator potentials from receptors on both the bag2 and chain fibers (Banks et al. 1997Go; Carr et al. 1998Go; Celichowski et al. 1994Go; Fallon et al. 2001Go). The primary afferent firing results from competition between these two impulse generators in which the dominant generator wins and suppresses all activity in the weaker generator by resetting. Although that would produce total occlusion, the mechanism responsible for partial occlusion is believed to involve spread of transduction current between the suppressed and dominant generator, resulting in increased impulse generation at the dominant site. Regardless of the exact mechanism, the effect has been well described experimentally and was incorporated into the model. The driving potentials produced by bag1 and combined bag2 and chain intrafusal fibers are compared and the larger of the two plus a fraction (S) of the smaller are summed to obtain the primary afferent firing. The secondary afferent output (which is not influenced by bag1 receptors) is obtained from the simple summation of the outputs of bag2 and chain intrafusal fiber models.

Implementation of the spindle model and parameter determination

The anatomical and mathematical structure of the spindle model was embodied as a set of nested blocks in the MATLAB Simulink modeling environment. First, a number of the model parameters were estimated and set directly from experimental measurements. The remainder of the parameters were then fit using standard least-squares fitting algorithms.

DIRECT PARAMETER ESTIMATION. The values of a number of model parameters, described below, were estimated directly based on a variety of previously published data. These parameters were estimated directly rather than leaving them as free parameters because in each case there was enough experimental evidence to justify such a direct approach.

The term "a," representing the nonlinear velocity dependency of afferent firing, was set to value 0.3 (Houk et al. 1981Go). We found a velocity power term of 0.3 (Houk et al. 1981Go) to produce a better fit when capturing the increase in the spindle primary's dynamic response during three velocities of ramp stretch than 0.5–0.6 as suggested by Prochazka and Gorassini (1998)Go.

Based on Boyd's figures, we estimated "L0SR," "L0PR," and "Lsecondary," which represent the sensory and polar regions' spring rest lengths and secondary afferent rest length (Boyd 1976Go). By dissecting the individual cat's spindle and examining it under the microscope, Boyd found that at the spindle's rest length, some 5% of its length belongs to the sensory region (where primary endings are located) and 95% to the polar region. The length of the secondary afferent ending that spans both the sensory and polar regions is roughly 5% of the total spindle's rest length. To obtain these percentages in terms of optimal fascicle length we estimated the spindle's rest length from the soleus fascicle slack length (0.8L0) (Scott et al. 1996Go). Therefore "L0SR," "L0PR," and "Lsecondary" were estimated to be 0.04L0, 0.76L0, and 0.04L0, respectively.

"R" represents the polar region length above which the lattice spacing of the myofilaments has effects on the cross-bridge kinetics that drive the force–velocity relationship. We assumed that this basis length would be similar for intrafusal and extrafusal fibers. In the case of extrafusal fibers, the effect of myofilament lattice spacing on cross-bridge kinetics was observed for fascicle lengths between 0.8L0 and 1.2L0, although the data suggest that the effect exists for the shorter fascicle lengths as well (see Fig. 3C in Brown et al. 1999Go). By extracting the suggested force–fascicle length curve, we estimated that this effect exists for fascicle lengths {gtrsim}0.5L0. To determine the polar region length at 0.5L0 fascicle length, we first needed to estimate the sensory region length. Because the evidence suggests that length changes in the sensory region are almost negligible compared with the length changes of the polar region (Boyd 1976Go), we assumed that at 0.5L0 fascicle length, the sensory region is approximately equal to its rest length (0.04L0). Therefore the polar region length at 0.5L0 fascicle length was estimated to be 0.46L0 ("R").


Figure 3
View larger version (16K):
[in this window]
[in a new window]
 
FIG. 3. Spindle model performance during 6-mm whole muscle ramp stretches. Primary afferent response at 3 different velocities (whole muscle stretches at 5, 30, and 70 mm/s; fascicle stretches at 0.11, 0.66, and 1.55L0/s; fascicle length changes 0.95–1.08L0) were performed in the absence of fusimotor stimulation (A, B, C), during constant dynamic fusimotor stimulation at 70 pps (D, E, F), and during constant static fusimotor stimulation at 70 pps (G, H, I). Solid thin lines represent model output; experimental data are shown as dots (·). Secondary afferent responses at 3 different velocities (whole muscle stretches at 5, 30, and 50 mm/s; fascicle stretches at 0.11, 0.66, and 1.12L0/s) are compared with the secondary afferent activity from 2 different muscle spindles (shown as +, · ) originating from the same experimental preparation (J, K, L).

 
"G," the term that relates the stretch of the intrafusal fiber's sensory endings to primary afferent firing, was estimated based on experimental data on cat tenuissimus muscle (Boyd 1976Go; Boyd et al. 1977Go) in the presence of maximal dynamic fusimotor stimulation; the bag1 sensory endings stretch by some 2–8% (scaling it to units of L0: 5% of primary sensory ending rest length = 0.05 x 0.04L0 = 0.002L0) while generating about 40 pps of the primary afferent firing. By combining these two values, the "G" term for bag1 fiber was estimated [G(for bag1) = 40 pps/0.002L0 {approx} 20,000 pps/L0]. Similarly, we obtained the "G" value for the combined bag2 and chain intrafusal fiber model, where maximal observed stretch during maximal static fusimotor stimulation of bag2 and chain was 12–30 and 15–20%, respectively (average 19% of sensory region rest length = 0.19 x 0.04L0 = 0.0076L0) (Boyd 1976Go; Boyd et al. 1977Go), and primary afferent firing about 150 pps (Boyd 1986Go) [G(for bag2&chain, primary) = 150 pps/0.0076L0 {approx} 20,000 pps/L0]. In the case of the secondary afferent endings, the maximal observed stretch during maximal static fusimotor stimulation of bag2 and chain fibers was comparable to the primary afferent case (Boyd 1976Go; Boyd et al. 1977Go), whereas the secondary firing observed at this stretch was about 110–115 pps (Boyd 1986Go) [G(for bag2&chain, secondary) = 110 pps/0.0076L0 {approx} 14,500 pps/L0].

The term "S" that represents the amount of partial occlusion that occurs in primary afferent firing also originated from direct experimental measurements (Fallon et al. 2001Go). Although there is some experimental evidence suggesting a length dependency to S, at this point, the data are very limited and noisy. As such, we have assumed that S is constant (S = 0.156), with the acknowledgment that if further data support an actual length dependency to S, then the model will need to be so modified.

The parameters that are used in mapping the fusimotor stimulation frequency to the model's fusimotor activation ("p," "freqbag1," "freqbag2," and "freqchain") were estimated based on Boyd's measurements of the intrafusal fiber's polar region contraction in response to different fusimotor stimulation frequencies (Boyd 1976Go). Because the chain fiber's measurements were obtained at relatively high frequencies, we used these data to estimate the Hill-equation parameters for chain fiber and then modified them for other two fibers to account for their different saturation frequencies. We should mention that, although limited experimental data suggest that fusimotor stimulation above the saturation frequencies [100 Hz for bag1 (fdynamic = 0.735) and bag2 (fstatic = 0.735); 150 Hz for chain (fstatic = 0.735)] will have little or no additional effect beyond that observed at the saturation frequency, our model actually has a slowly increasing effect because of the nature of the Hill equation. Although it was possible to correct this by introducing some nonlinearity in the model we chose not to do so to avoid complications during the model inversion (see DISCUSSION). We believe that this does not represent a limitation of our model because fusimotor neurons usually fire at frequencies below the saturation level (see DISCUSSION).

DATA SELECTION FOR FREE PARAMETER OPTIMIZATION. In a series of experiments during the 1960s and 1970s, the activities of identified primary and secondary afferents were recorded from isolated or partly isolated spindles from soleus muscle of the cat. We used essentially the same published records that were used in previous spindle modeling studies (i.e., Lin and Crago 2002Go; Maltenfort and Burke 2003Go; Schaafsma et al. 1991Go) to represent afferent activity during a broad range of experimentally controlled kinematic and fusimotor conditions. Afferent activity records during different velocities of ramp, triangular, and sinusoidal stretches under various constant fusimotor drives were used to develop and tune the model (Crowe and Matthews 1964Go; Hulliger et al. 1977aGo,bGo; Lennerstrand 1968Go; Lennerstrand and Thoden 1968aGo,bGo; Matthews 1963Go; Prochazka 1996Go). Once the selected spindle activity patterns were digitized manually from the figures in the journal articles, MATLAB's Spline Toolbox was used to facilitate parameter optimization, as described below.

The independent data set used for validation of the model originated from the medial gastrocnemius (MG) muscle of the cat rather than the soleus muscle and was similarly digitized manually (Taylor et al. 2000Go).

An important modification to the above-described data was made to facilitate curve-fitting. The ramp stretch records that were used to develop the model (Crowe and Matthews 1964Go) were modified in several ways because of their noisiness during the dynamic part of the stretch (which the authors attributed to experimental artifact). This was accomplished by using a cleaner ramp stretch record collected during a single velocity of stretch (Prochazka 1996Go) and scaling it to match the record recorded by Crowe and Matthew. First we divided Prochazka's afferent record into static and dynamic phases. We assumed that during the stretch the static component increased linearly from initial firing before the stretch to the final level 2–2.5 s after the stretch. The dynamic response was obtained by subtracting this static component from the afferent firing record. Although the static response of Prochazka's data did not require scaling, the dynamic response was scaled to match the dynamic response of afferent activity recorded by Crowe and Matthew. The rest of the data used in tuning and validating the model did not require any such alterations or scaling.

Because the length records from the experiments used in developing the model represent whole muscle length, appropriate calculations were done using an extrafusal model of cat soleus muscle (Brown et al. 1999Go) to determine the relative contributions of tendon stretch and muscle fascicle (i.e., spindle) length (see APPENDIX 1, B and C).

FREE PARAMETER OPTIMIZATION. Free parameter optimization was accomplished using the Levenberg–Marquardt optimization method (Press et al. 1986Go). Data records of ramp, triangular, and sinusoidal stretches in the absence or during constant static and/or dynamic fusimotor stimulation were used to obtain a single set of model parameters for the simulations presented throughout this report.

A single intrafusal fiber model of bag2 and chain fibers was assumed during this optimization because the bag2 and chain intrafusal fiber contributions to the activity of the primary afferent could not be quantitatively distinguished during the experiments in which constant fusimotor activity was used. Afterward, the combined bag2 and chain model was separated based on the simplifying assumption that they share exactly the same structure and parameter values, while having different temporal properties in response to step changes in fusimotor drive (see Fusimotor activation). This approach required two parameters relating to static fusimotor stimulation (beta2 and {Gamma}2) to be scaled to accommodate their different fusimotor saturation points (see Table 1).

Initial parameter optimization runs produced very similar values of KSR and KPR for both the bag1 and combined bag2 and chain fibers. We therefore made the simplifying assumption that KSR and KPR were identical for all intrafusal fiber types. Likewise, CS was found initially to be very similar for the different fiber types and so it, too, was assumed to be a single constant. These simplifying assumptions left ten free parameters as necessary to capture fully the primary afferent activity and two additional parameters to capture secondary afferent activity (Table 1).

In the case of primary afferent parameter optimization, the optimization was broken down into two stages to reduce the number of free parameters optimized simultaneously. The first set of optimizations involved optimization of eight parameters [KSR, KPR, beta0(bag1), beta0(bag2&chain), beta2(bag2&chain), {Gamma}2(bag2&chain), CS, LNSR] for data recorded either in the absence of fusimotor stimulation or during constant static fusimotor stimulation. This was done to introduce just enough dynamic sensitivity into the combined model of the passive bag2 and chain fibers so that it could be reduced appropriately during the static fusimotor stimulation. Hundreds of optimization runs starting from randomly selected initial points were performed to avoid local minima. After this initial set of optimizations, the parameter values that generated the best fits were used in the next set of optimizations where data involving dynamic fusimotor stimulation were used and two more parameters were optimized [beta1(bag1), {Gamma}1(bag1)]. In addition to using the best-fit parameter values, randomly chosen parameter values in the proximity of the best-fit values were used as well to run a number of optimizations of all ten primary afferent parameters simultaneously (Table 1).

In the case of the secondary afferent model, two additional parameters ("X" and "LNPR") had to be introduced, explained by the fact that part of the afferent transduction region lies within the polar region of the intrafusal fiber. These two parameters were optimized on the secondary afferent record during 2- and 8-mm/s triangular stretches in the presence of constant static fusimotor stimulation ({gamma}static = 70 pps). The reason for not using other secondary afferent records (5-, 30-, and 50-mm/s ramp stretches and 2- and 8-mm/s triangular stretches in the absence of fusimotor stimulation) is attributed to the existence of several traces of secondary afferent activity during the same kinematic and fusimotor conditions. "LNPR", the term representing the polar region length above which extension of polar region's secondary afferent endings takes place (see Eq. 8), produced the best fit at the value of 0.89L0. Interestingly, the parameter "LNPR", when summed with sensory region rest length ("L0SR"), resembles approximately the intrafusal fiber length (0.89L0 + 0.04L0 = 0.93L0) at which the passive tension develops. Although this value was never previously measured in the case of intrafusal fibers, the extrafusal fiber data indicated that passive tension develops at similar sarcomere lengths (Brown et al. 1996aGo).

After the parameter optimization, the intrafusal fiber mass (M) was added to individual intrafusal fiber models to increase the model's stability; it had no significant influence on the model's properties.


    RESULTS
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX 1
 ACKNOWLEDGMENTS
 REFERENCES
 
In the following subsections we show the model's prediction compared with data for each type of movement used in the parameter optimization. This is then followed by validation of the model against novel data that were not used for parameter optimization.

Ramp-and-hold stretches

The model's ability to reproduce primary afferent activity during ramp-and-hold stretches is demonstrated and compared with the experimental data in Fig. 3 (Crowe and Matthews 1964Go). Whole muscle stretches of 6 mm length (about 0.95L0 to 1.08L0 fascicle length stretches; see APPENDIX 1, B and C) at three different velocities (5, 30, and 70 mm/s whole muscle or about 0.11, 0.66, and 1.55L0/s fascicle length stretches) were used to run the spindle model simulations in which the fusimotor drive was kept either at zero (Fig. 3, A, B, and C) or at a constant value (Fig. 3, D, E, and F: {gamma}dynamic = 70 pps; Fig. 3, G, H, and I: {gamma}static = 70 pps). In the absence of fusimotor stimulation, the model (solid lines) closely agreed with the experimental records of instantaneous primary afferent frequency (dots) for 5 and 30 mm/s, whereas it underestimated the peak dynamic response of the 70-mm/s record by some 20 pps (about 18%). Because the researchers reported having problems with 70-mm/s records and because their published records would suggest an almost linear relationship between dynamic response and velocity that was subsequently refuted (Houk 1981), all 70-mm/s records (with and without any fusimotor stimulation) were excluded from the free parameter optimization. The presence of dynamic fusimotor stimulation during three velocities of stretch increased the dynamic response of the biological spindle, whereas the presence of constant static fusimotor drive had the opposite effect while introducing a strong static bias (Fig. 3, G, H, and I). The primary afferent's dynamic response for 5 and 30 mm/s during static fusimotor stimulation was decreased, whereas the 70-mm/s record remained unaffected. Other published records of primary afferent activity show a consistent decrease in dynamic behavior in response to static fusimotor drive at all velocities, providing further justification to exclude these 70-mm/s records from parameter optimization. Finally, small oscillations in the behavior of the model (e.g., at the end of stretch in Fig. 3, F and H) were artifacts of the abrupt velocity changes; these can be significantly reduced or eliminated by making such transitions smoother, as they would be in the actual spindle.

Predictions of the secondary afferent model during ramp-and-hold stretches were compared with the secondary afferent data records of two different muscle spindles (Matthews 1963Go) (Fig. 3, J, K, and L). As seen from the figure, there is great variability in the original data among spindles. Instead of optimizing the X and LNPR for these records, we used the values (0.7 and 0.89L0) that were optimized on the 2- and 8-mm/s triangular stretches in the presence of constant static fusimotor drive where only a single secondary afferent trace existed. In the case of three velocities tested (whole muscle velocities: 5, 30, and 50 mm/s; fascicle velocities: 0.11, 0.66, and 1.12L0/s), predictions of the spindle model are contained within the range of the observed variability among secondary afferent firing (see DISCUSSION).

Triangular stretches

The model's ability to capture primary afferent activity during 8-mm whole muscle triangular stretches (fascicle length: 0.90L0–1.08L0) at ±8 mm/s (fascicle velocity: ±0.18L0/s) was evaluated using data reported by Lennerstrand and Thoden (Lennerstrand 1968Go; Lennerstrand and Thoden 1968aGo,bGo). Lennerstrand and Thoden did not identify from which extensor muscle (soleus or lateral gastrocnemius) the afferent activity originated, but they suggested that the two muscles produced quantitatively similar afferent activity. Triangular stretches in the presence of either dynamic or static fusimotor drives at 35, 70, and 200 pps were available (Fig. 4). The primary afferent model accurately accounted for the spindle's behavior during lengthening (similar to the ramp stretches in Fig. 3), as well as during the shortening case in which the afferent firing was silenced when influenced by dynamic fusimotor stimulation or was maintained when under static fusimotor stimulation. The model performed well during the static fusimotor stimulation. In the case of dynamic fusimotor stimulation, it overestimated the initial dynamic response for the cases of 70- and 200-pps stimulations, and failed to reproduce the recovery of the low-level afferent activity during the second half of shortening in the presence of 200-pps dynamic fusimotor stimulation.


Figure 4
View larger version (12K):
[in this window]
[in a new window]
 
FIG. 4. Model's ability to capture primary afferent activity during triangular stretches. Whole muscle (8 mm) stretches (fascicle length changes 0.90–1.08L0) at 8 mm/s (fascicle stretches at 0.18L0/s) were performed during constant dynamic (A, B, C) or static (D, E, F) fusimotor stimulations at 3 different frequencies (35, 70, and 200 pps). Solid thin lines represent model output; experimental data are shown as dots.

 
The model's ability to reproduce secondary afferent firing during triangular stretches (Lennerstand and Thoden 1968aGo,bGo) is shown in Fig. 5. Whole muscle triangular stretches (2 and 8 mm/s, corresponding to fascicle velocities 0.05L0/s and 0.18L0/s) were analyzed in the absence of fusimotor drive (Fig. 5, A and B), and in the presence of static fusimotor stimulation (Fig. 5, C and D). The model accurately accounted for the maintained secondary afferent firing during shortening, in contrast to silencing of the primary afferent during this phase of triangular stretch. Because X and LNPR values were optimized on the triangular stretches in the presence of static fusimotor drive, it is not surprising that the model was very accurate in capturing those records; in the absence of fusimotor stimulation its predictions were contained within the range of the observed variability among secondary afferent firing. Finally, it should be mentioned that this particular set of secondary afferent activity most probably originated from the cat's ankle flexor muscles rather than the extensor muscles (soleus and medial gastrocnemius) that were used in other triangular stretches involving primary afferent activity.


Figure 5
View larger version (11K):
[in this window]
[in a new window]
 
FIG. 5. Model's ability to capture secondary afferent activity during triangular stretches. Whole muscle (8 mm) stretches (fascicle length changes 0.90–1.08L0) at 2 and 8 mm/s (fascicle stretches at 0.05 and 0.18L0/s) were performed in the absence of fusimotor stimulation (A, B) and during static fusimotor stimulation at 70 pps (C, D). Solid thin lines represent model output; experimental data are shown as dots.

 
Sinusoidal stretches

Primary afferent activity and the model's performance during whole muscle sinusoidal stretches (1,400 µm peak-to-peak, corresponding to fascicle length 0.995 ± 0.012L0) at a frequency of 1 Hz are shown in Fig. 6 (Hulliger et al. 1977aGo,bGo). In the first set of experiments, such stretches were applied during constant static or dynamic fusimotor drives with stimulation frequencies of 50, 75, and 125 pps (Fig. 6, A and B). The published database for sinusoidal stretches also includes simultaneous stimulation of dynamic and static fusimotor efferents. In one set of simulations, dynamic fusimotor drive was held constant at 125 pps, whereas static fusimotor drive was either 50, 75, or 125 pps (Fig. 6C) and in the other case static was held at 70 pps, whereas dynamic fusimotor activity was either 50, 75, or 125 pps (Fig. 6D).


Figure 6
View larger version (21K):
[in this window]
[in a new window]
 
FIG. 6. Model's ability to capture primary afferent activity during sinusoidal stretches and constant fusimotor stimulation. Whole muscle stretches of 1.4 mm peak-to-peak (fascicle length changes 0.995 ± 0.012L0) at 1 Hz. A: dynamic fusimotor stimulation at 50, 75, or 125 pps. B: static fusimotor stimulation at 50, 75, or 125 pps. C: dynamic fusimotor stimulation at 125 pps plus static fusimotor stimulation at 50, 75, or 125 pps. D: static fusimotor stimulation at 70 pps, plus dynamic fusimotor stimulation at 50, 75, or 125 pps. Solid thin lines represent model output; experimental data are shown as +, ·, and *.

 
The records in Fig. 6 are particularly important to demonstrate the effects of occlusion between primary transduction zones. A model with no occlusion (i.e., simple summation) predicts much higher activity in response to combined static and dynamic drive. For example, peak afferent activity during stretch with either pure dynamic drive (Fig. 6A) or pure static drive (Fig. 6B) was around 100 pps; with combined drive at similar rates (Fig. 6, C and D) it was only slightly higher. A model with complete occlusion (i.e., winner-take-all) would not have demonstrated the complex transition that occurs at the beginning of the stretch phase in the presence of combined static and dynamic drive (Fig. 6, C and D). Activity during the shortening phase requires static drive because the viscosity of the bag1 fiber causes its transduction zone to be slack. When the spindle starts to lengthen in the second half of this sinusoidal motion, this same viscosity, amplified by any concurrent dynamic drive, accentuates the stretch of the bag1 transduction zone, which then modulates but does not yet dominate the activity arising from the bag2 and chain transduction zones. When stretch velocity peaks, the afferent activity tends to be dominated by the dynamic drive, although residual contributions from different levels of static drive are still apparent. In Fig. 6, A and C, it is apparent that the model has overestimated the effects of dynamic drive at the highest frequency of 125 pps, perhaps because this particular bag1 fiber had a somewhat lower saturation frequency than that used in our model (100 Hz). It is important to remember that all of the experimental data from all of the experimental preparations and paradigms in the cited literature were fit with a single set of model parameters.

Model validation

During the early stages of this model's development, the bag2 and chain fibers were combined into one system because of unavailability of data that would allow quantitative distinction between their properties and because of their similar responses during constant gamma static drive. The bag2 and chain fibers are distinguished by their different dynamic responses to the onset and offset of fusimotor drive (modeled as low-pass filter intrafusal fiber properties), but these effects had not been studied systematically in recordings of afferent activity during controlled fusimotor stimulation. To test the accuracy of the low-pass filter property of the bag2 intrafusal fiber included in the model, we used the secondary afferent activity recorded from MG muscle where direct recordings of two types of static fusimotor activity during decerebrate locomotion in the cat and secondary afferent activity were successfully obtained (Taylor et al. 2000Go). Although the existence of two types of static fusimotor drives is still not universally accepted, the direct fusimotor recordings of Taylor et al. provided convincing new evidence suggesting that type-1 static fusimotor drive innervates chain or bag2 and chain fibers together, whereas type-2 static fusimotor drive innervates solely the bag2 fiber.

In the experiments used for our model's validation, Taylor et al. (2000)Go used a decerebrate cat preparation where locomotor movements of three legs occurred in contact with the treadmill. The experimental leg was denervated except for MG and tibialis anterior (TA) muscles, kept clear of the belt by fixation at midfemur and at the lower end of tibia and allowed to freely rotate about ankle through the phasic contractions of MG and TA. Single-unit recordings were obtained simultaneously from type identified spindle afferents and fusimotor efferents. To reveal the effects of fusimotor drive, Taylor et al. (2000)Go plotted the secondary difference signal. The secondary difference signal is defined as the difference between the secondary afferent activity recorded during locomotor cycling of the intact MG (where two types of static fusimotor drives were present) and the activity recorded with the leg deefferented (no fusimotor or extrafusal drive), while the ankle was moved through the same trajectory as recorded initially (secondary difference in Fig. 7). We used this experimentally obtained secondary difference signal and compared it to the predicted secondary difference of our model under the same conditions (meaning: secondary activity during same kinematics and simultaneous stimulation of two static fusimotor drives minus the secondary activity during same kinematics and no fusimotor drive), both with (Fig. 7B) and without (Fig. 7A) the inclusion of the low-pass filter for bag2 intrafusal fiber model (X and LNPR were set to 0.7L0 and 0.89L0, respectively).


Figure 7
View larger version (20K):
[in this window]
[in a new window]
 
FIG. 7. Model's prediction of the secondary difference signal of secondary afferent during the cat's locomotor step. Secondary difference signal is the difference between the secondary afferent activity recorded during locomotor cycling of the intact MG and the activity recorded with the leg deefferented (no fusimotor or extrafusal drive) while the ankle was moved through the same trajectory as recorded (thin solid lines). In this figure the model's prediction of secondary difference signal during simultaneous stimulation of 2 types of static fusimotor drives [type 1 ({triangleup}) and type 2 ({square}), left frequency scale] is compared with experimental data (Taylor et al. 2000Go) (+, left frequency scale). Although the amplitude stretch ({alpha}) was about 7° of ankle rotation, the mean cycle period lasted for about 0.83 s. Performance of the model lacking low-pass filter property of bag2 fiber is shown in A (thick interrupted line, right frequency scale), compared with the performance of the model when low-pass filter was introduced in bag2 intrafusal fiber model (B, thick solid line, right frequency scale). C: magnified portion of the cycle (enclosed areas in A and B) where the greatest difference in the performance of the model without and with low-pass filters occurred.

 
Based on the ankle angle kinematics provided during decerebrate locomotion it was not possible to calculate the MG fascicle length, but instead a reasonable guess about these values was made. Records suggest that during normal walking, the MG usually operates between the lengths of 0.94L0 and 1.14L0, which typically corresponds to some 45° of ankle extension (Goslow et al. 1973Go). Based on the suggested limb's position we assumed that the MG was somewhere in the vicinity of the optimal fascicle length and chose the fascicle length range for the smaller recorded movements to be between 1.05L0 and 1.085L0 because it produced the best fit for the purpose of modeling the secondary difference signal. Furthermore, a typical spindle contains four to ten chain fibers rather than a single chain fiber (as our model does) which are driven similarly but asynchronously by separate type-1 static fusimotor neurons. Thus our model's secondary difference prediction had to be scaled to fit the experimental data (experimentally obtained secondary difference signal (+ +): left frequency scale; modeled secondary difference signal [without low-pass filter (---); with low-pass filter (—): right frequency scale)]). The inclusion of the low-pass filter in the bag2 fiber model improved the model's performance (Fig. 7C), although a small discrepancy remains at the onset of muscle shortening.


    DISCUSSION
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX 1
 ACKNOWLEDGMENTS
 REFERENCES
 
Historically the muscle spindle has received much attention because of its importance in sensorimotor control. We chose to build our model of the spindle based on the extensive information available about the internal components and their properties as measured directly in intrafusal fibers, as inferred indirectly from spindle activity during controlled experiments, or as extrapolated from general properties of striated muscle fibers. We are pleased that a single model with one set of coefficients could account for spindle activity from many different preparations under widely varying conditions of motion and fusimotor drive.

Our analysis of the model features and results presented has focused on their ability to capture qualitative behavior rather than quantitative measurements of goodness of fit (e.g., Figs. 6 and 7). Qualitative comparisons are often more instructive than overall quantitative comparisons, which tend to be dominated by portions of the records that happen to have relatively little modulation. Many estimates were required to pool data from various preparations and their disparate descriptions in the literature. Errors in these estimates or heterogeneity in the spindles themselves would be expected to produce offsets in afferent activity levels that would shed little light on the ability of the model to capture the physiological properties of muscle spindles.

An important feature of afferent activity that our model describes accurately is the partial occlusion effect that has been identified in the recent experimental literature (Banks et al. 1997Go; Carr et al. 1998Go; Fallon et al. 2001Go). Modeling of this effect between bag1 and combined bag2 and chain intrafusal fiber models reconciled the primary afferent activity during combined static and dynamic fusimotor drive. A somewhat weaker partial occlusion effect has also been reported in the biological spindle during simultaneous stimulation of multiple gamma static motoneurons (Carr et al. 1998Go; Fallon et al. 2001Go) but the interpretations of these results are less clear. One suggestion is that the bag2 fiber drives one impulse generator site, whereas all the other chain fibers (four to 11 of them) collectively drive the other, or that every single chain fiber has a separate impulse-generating site as well. This seems unlikely given their proximity. The nonlinear summation observed during activation of multiple static fusimotor efferents may instead reflect nonlinear summation of the generator potentials themselves. Because of the uncertainty of interpretation of partial occlusion observations during stimulation of multiple static fusimotor efferents and because many researchers still recognize only two impulse-generating sites (one on the bag1 and the other on the bag2 and chain fibers), we excluded other occlusion effects in our model.

Our model accounted well for the distinctive activation dynamics of bag2 versus chain fibers, but there are other differences in the mechanical properties of these intrafusal fibers that were omitted from modeling because of insufficient data. Evidence suggests that the bag2 fiber is more viscous than the chain fiber [the bag2 fiber contains slow twitch myosinATPase, whereas the chain fiber has the fast twitch form (Ovalle and Smith 1972Go)], although the afferent records demonstrating this distinction are very noisy. Therefore instead of making arbitrary guesses regarding viscosities of the fibers, we assumed the same viscosities and structures for two intrafusal fiber models but different fusimotor saturation points and temporal properties. Note, however, that the basic structure of the model (Fig. 1) assumes only a single source of static fusimotor drive controlling the bag2 and chain fibers together, whose composite mechanical properties appear to be well captured by our model.

Our model accounted well for the reported secondary afferent behavior. Because its parameters were optimized on the secondary afferent activity record during two ramp stretches (2 and 8 mm/s) in the presence of static fusimotor activity, it accurately captured those records; in the case where more than a single afferent record was available for the same kinematic and fusimotor conditions, the model predicted the values within the range of observed variability. The reason for the existence of large variability among secondary afferent records is not well understood. Initially we thought that such variability might reflect the variations in the secondary afferent location along different intrafusal fibers. If so, then it could be eliminated by adjusting the X parameter value. However, we realized that this was not sufficient and that it was also necessary to adjust the threshold length of the polar region for transmitting stretch to the secondary ending (LNPR) to capture the experimental variability. One anatomical feature related to the LNPR term could be the presence of kinks within the chain fibers, especially in the polar region areas where secondary afferent endings reside (Boyd 1976Go). A kinked chain fiber might need to be stretched more (i.e., larger LNPR) and pulled straight before its secondary endings are stretched. To account for variability among individual secondary afferent data, both X and LNPR parameters probably need to be adjusted.

Finally, it should be mentioned that the model differs from the biological spindle by omitting two effects of movement history: the response to noncycling cross-bridges and the effects of prior movement on the response of the receptors that are unrelated to noncycling cross-bridges (Nichols and Cope 2004Go). The response to noncycling cross-bridges is believed to give rise to the effect of "stiction" or the initial burst (Proske and Morgan 1999Go). This mechanism is believed to account for the experimentally observed phenomenon of very high viscosity in the quiescent bag1 fiber for a very small range of motion, purportedly as a result of the existence of a small number of residual cross-bridges between myofilaments, particularly in the absence of background fusimotor drive (Hill 1968Go). The consequences of omitting this effect are apparent during very slow and small-amplitude sinusoidal stretches in the absence of fusimotor stimulation, where the model noticeably underestimated the primary afferent activity (not shown), as well as during initial moments of stretch, where an initial burst in afferent firing was not captured by the model (Fig. 3, B and C). We chose not to model the "stiction" because it seems unlikely to be a factor during most natural motor behavior, in which spindles usually operate over longer stretches and with continuously modulated fusimotor input, both of which eliminate stiction.

Haftel et al. (2004)Go recently described the phenomena of reduced dynamic response (RDR), in which immediately successive trials of triangular-stretch release produced a systematic reduction in the dynamic response in deefferented spindles in rat triceps surae muscle. In their experiments, the RDR and initial burst properties of the rat's spindle afferents were studied and it was concluded that the initial burst and RDR exhibit very different behavior during successive trials of muscle-stretch release. For example, the initial burst is expressed over a much shorter time and smaller magnitude of muscle stretch than RDR. Additionally, the initial burst recovers more rapidly than RDR and its magnitude is variable and occasionally absent in repeated sets of stretch-release trials, whereas the greater dynamic response is always present in the first trial compared with subsequent trials. At the moment, there are various explanations for the differences between the initial burst and RDR response properties (Haftel et al. 2004Go). Interestingly, using a similar triangular-stretch–release paradigm in the decerebrate cat where fusimotor activity was present, Houk et al. (1992)Go found no evidence of change in the dynamic response. If this history-dependent effect occurs only in the absence of fusimotor tone, then its omission should not give rise to significant errors during most conditions of normal use, which typically include substantial fusimotor modulation.

Comparison with previous modeling attempts

Several models of spindle afferent activity have been developed to account for the growing base of physiological data (Chen and Poppele 1978Go; Crowe 1968Go, 1970Go; Hasan 1983Go; Houk et al. 1981Go; Lin and Crago 2002Go; Maltenfort and Burke 2003Go; Matthews and Stein 1969Go; Rudjord 1970Go; Schaafsma et al. 1991Go). Prochazka and Gorassini (1998)Go compared the performance of some of these models to experimentally recorded primary and secondary afferent activity from cat hamstring muscles during locomotion. Because the afferent records used for this purpose incorporated intact fusimotor drive and the spindle models did not model these effects, certain assumptions regarding the fusimotor activity had to be made. We chose not to conduct a similar comparison analysis for the purpose of testing the performance of our model because accurate modeling of the fusimotor effects was our major goal. Instead, we chose to concentrate on comparisons only with those models incorporating fusimotor effects (Hasan 1983Go; Lin and Crago 2002Go; Maltenfort and Burke 2001; Schaafsma et al. 1991Go).

Hasan (1983)Go presented one of the first models consisting of mathematical elements representing the anatomical components found in the biological spindle. The model included a sensory region, which was modeled as a spring, and a polar region, which included nonlinear velocity dependency and a property akin to friction. It addressed issues such as the initial burst response (which we ignored) and fusimotor modulation of afferent activity. Because of both its innovative nature and its easy implementation, Hasan's model attracted the attention of many researchers and its performance has been tested extensively. A major drawback of this model, however, arises from its rigid structure, which requires changing its parameters to account for each fusimotor state. Furthermore, during a series of stretches, the modeled initial burst in afferent activity persists at the start of each stretch, whereas in the biological spindle it is visible only during the first stretch after a period of rest (Matthews 1972Go).

The next major step in spindle modeling came from Schaafsma et al. (1991)Go, who designed a complex structural model of only the primary afferent that incorporated fusimotor effects. The model consisted of two intrafusal fiber submodels (bag1 model and combined bag2 and chain model) and included thixotropic effects, although a mechanism for recovery from these effects was not provided. Subsequently the model was extended to incorporate a simulated chain fiber (Scheepstra et al. 1995Go). Overall, the model performed well, especially during the ramp-and-hold stretches for which the model parameters were optimized. Limitations in its performance were apparent for sinusoidal stretches (particularly small-amplitude stretches), for which the model underestimated afferent activity. Because the model predates the elucidation of the effect of partial occlusion during simultaneous stimulation of dynamic and static fusimotor efferents (Banks et al. 1997Go; Carr et al. 1998Go; Fallon et al. 2001Go), the model assumes unrealistic, extreme occlusion. In other words, it assumes that the afferent activity will be controlled completely by the larger of the two fusimotor inputs. Although the developers of this model do not report on its performance during conditions of simultaneous static and dynamic fusimotor activation, it is likely that the model would incorrectly predict the afferent firing during those conditions as a consequence of its incorrect occlusion model. Fina