# Proprioceptors and Models of Transduction

Post-publication activity

Curator: Gerald E. Loeb

Even when deprived of exteroceptive sensory information such as vision and touch, we are aware of the posture and motion of our bodies (kinaesthesia) and the amount of effort being exerted by our muscles. Most of this information comes from large diameter, myelinated, afferent nerve fibers innervating specialized mechanoreceptors called muscle spindles and Golgi tendon organs, whose structure and function have been modelled in detail as reviewed herein. Muscles, tendons, ligaments and joint capsules contain a host of other sensory endings serviced by smaller diameter and more slowly conducting nerve fibers, some of which are also sensitive to mechanical stimuli such as strain (Paintal, 1960, Rymer et al., 1979, Grigg et al., 1986, Martin et al., 2006). Their structure is more heterogeneous and their roles in sensorimotor control are less well understood.

## Muscle spindle

The muscle spindle is a sense organ found in most vertebrate skeletal muscles. In a typical mammalian lower limb muscle, 20-500 muscle spindles can be found lying in parallel with extrafusal muscle fibers and experiencing length changes representative of muscle length changes (Chin et al., 1962, Voss, 1971). Their sensory transducers (primary and secondary afferents) provide the central nervous system with information about the length and velocity of the muscle in which the spindle is embedded (Figure 1
Figure 1: A typical 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, 1980)).
). Their sensitivity is dynamically modulated by a complex fusimotor system of tiny intrafusal muscle fibers under separate control by the central nervous system. The spindles provide the main source of proprioceptive feedback for kinesthetic perception and spinal coordination.

Early spindle literature described very complex firing patterns from spindle afferents when subject to various kinematic and fusimotor conditions. Various stretches (ramp, triangular, sinusoidal) at different speeds, starting from different fiber lengths or during different fusimotor activations were utilized to help elucidate the mechanisms responsible for the observed afferent behaviour (transducer properties, mechanical properties, etc). Early models of muscle spindles usually addressed only specific afferent behavior such as small amplitude sinusoidal stretches that were hypothesized to be adequately described by linear system (Matthews and Stein, 1969, Poppele and Bowman, 1970, Chen and Poppele, 1978). They proposed transfer functions to describe the system and in some cases also discussed potential ways to integrate fusimotor activity. While accurately capturing the specific range of behaviour, such models failed to capture the nonlinear afferent behaviour during larger amplitudes and higher velocities of stretch (Houk et al., 1981). With time, the understanding of spindle structure and function became more complete and the models became more detailed and complex. We first review models of spindle transduction that are based on their anatomical structure, followed by a black box model by Maltenfort and Burke (Maltenfort and Burke, 2003).

## Spindle components

The structural models rely on anatomy and physiology to explain spindle afferent behaviour during different kinematic and fusimotor conditions. The structural spindle models usually use the extensive information available about the internal spindle components and their properties as measured directly in intrafusal fibers, or as inferred indirectly from spindle activity during controlled experiments, or as extrapolated from general properties of striated muscle fibers. Similarly to the biological spindle, the spindle models are driven by a combination of length changes and fusimotor activation (where modelled) in order to generate primary and/or secondary afferent activity. Most models use the extrafusal fiber length as the kinematic input, assuming that it is representative of the length changes being experienced by a spindle lying in parallel to extrafusal fibers. Because of series elasticity in tendons, the velocities of the spindles and extrafusal muscle fibers may differ transiently from the length of the overall musculotendon unit from origin to insertion, which more directly reflects the skeletal posture (Hoffer et al., 1992). The models address either both the primary and the secondary afferents activity (Lin and Crago, 2002, Mileusnic et al., 2006) or concentrate only on the primary afferent (Rudjord, 1970, Hasan, 1983, Schaafsma et al., 1991).

### Intrafusal fiber models

The structural spindle models include explicit models of intrafusal fibers. All the spindle models include a model representing a bag1 intrafusal fiber, while a bag2 and chain fibers are either lumped together in a single model (Schaafsma et al. 1991; Rudjord et al. 1970; and partly Mileusnic et al 2006) or modelled separately (Lin and Crago 2002; Hasan, 1983; partly Mileusnic et al., 2006). In a particular model system, all the intrafusal fiber models typically share the same structure, but their parameters and fusimotor innervation (when included) depend on fiber type. The intrafusal fiber is generally modelled as consisting of a polar and sensory region in series.

In a biological spindle, the centrally located sensory (transduction) region contains afferent endings wrapped around it. This region is commonly modeled as an elastic element (spring). The stretch in this region is related to the distortion of afferent endings and therefore depolarization of the membranes and generation or contribution to generation of the action potential firing (depending on the model).

On each side of the sensory region are the polar regions which share many properties with well-studied striated skeletal (extrafusal) muscle fibers. When modelling the polar region, most of the models take advantage of this similarity and utilize Hill-type muscle models or typical force-length and force-velocity curves obtained for extrafusal fibers to model the polar region’s behavior (Hill et al. 1968). Therefore, the polar region is usually modeled as a passive spring and a contractile element consisting of an active force generator and a damping element. In models where fusimotor effects are included, the polar region receives fusimotor activation, dynamically changing its mechanical properties. Dynamic fusimotor activation innervating the bag1 fibers has its major effect on increasing the stiffness of the myofilaments, thereby increasing the velocity sensitivity of the primary afferent. Static fusimotor activation of the bag2 and chain fibers induces polar region contraction, effectively introducing a positive bias to the activity of primary and secondary afferents.

Several models also use the polar region to model the experimentally observed effect of an initial burst (‘stiction’) in afferent firing at the start of small stretching movements from rest (Schaafsma et al. 1991; Lin and Crago 2002; Rudjord 1970; Hasan, 1983; Matthews and Stein, 1969; Chen and Poppele, 1978). This seems likely to be caused by a very high viscosity in the quiescent bag1 fiber that exists for a very small range of motion, as a consequence of a small number of residual cross-bridges between myofilaments that form in the absence of background fusimotor drive (Hill, 1968, Chen and Poppele, 1978). In addition to being responsible for the initial burst, this phenomenon could be also a factor responsible for the observed elevated afferent activity during very slow and small-amplitude sinusoidal stretches. The effect is modelled differently depending on the model. The model by Schaafsma et al. (1991), for example, assumes the existence of one hundred cross-bridge regions that are initially fixed and thereby make the muscular part rigid. When a stretch is introduced, all elongation is initially transmitted to the sensory region resulting in intrafusal fiber force. As the force increases, the cross-bridges gradually release when their individual breaking forces are reached. The model by Lin and Crago (2002) represented this effect as a potentiation effect (tx) that is a negative (decreasing) exponential function of the activation level. Although several models omit this “stiction” effect, it is debatable whether that represents a significant limitation because the effect 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.

Like all muscle fibers, intrafusal fibers exhibit a substantial delay between fusimotor excitation and the generation of contractile output in the myofilaments. This is particularly relevant given that some theories of fusimotor function propose rapid modulation of fusimotor activation such as during a specific phase of the locomotor cycle. In the two models addressing this aspect (Lin and Crago 2002; Mileusnic et al. 2006), different types of fusimotor fibers were given different time-constants for a low-pass filter between excitation and mechanical output, as reported by (Boyd, 1976) and (Boyd et al., 1977). While both models assume short time constant for the rapidly contracting chain fibers and longer for bag1 and bag2 fibers, they disagree regarding the time constants for the bag fibers. Boyd et al. (1977) reported a very long time-constant for bag1 fibers based on video observations of visible contractions, but this may under-estimate the rise-time for the viscosity increase in these tonic-like muscle fibers that underlies the velocity sensitivity of the primary afferents.

### Primary and secondary afferent

Because the primary afferent endings spiral around all three intrafusal fibers, most of the models assume that the primary afferent ending reflects stretches in the sensory zones of all intrafusal fibers, except for a model by Lin and Crago (2002) where only the contributions from bag1 and bag2 fibers are recognized.

The few spindle models that address the secondary afferent model its activity as resulting from the bag2 and chain fibers where its spiral endings are located (Hasan, 1983; Lin and Crago 2002; Mileusnic et al. 2006). In these models, the fiber’s potential to contribute to the afferent firing is obtained either by calculating the stretch in the sensory region (Lin and Crago 2002; Hasan, 1983) or the stretch in both sensory and polar regions (Mileusnic et al. 2006) to account for the juxtaequatorial location of the secondary afferent endings straddling both sensory and polar regions of bag2 and chain intrafusal fibers. While assuming that the sensory endings are located only on the central sensory region might be acceptable during muscle lengthening where the primary and secondary afferents behave similarly, it might introduce errors during shortening where very different behaviors have been observed in primary and secondary afferents (especially during low or absent static fusimotor stimulation, when the primary afferent becomes silent, whereas the secondary afferent continues to fire)(Hulliger et al. 1977).

Early spindle models computed the resulting primary or secondary afferent activity as a simple summation of the relevant intrafusal fiber contributions (Rudjord, 1970a Rudjord, 1970b, Hasan, 1983), while the more recent models assume the existence of multiple spike-generating sites (Schaafsma et al. 1991; Lin and Crago 2002; Mileusnic et al. 2006). Multiple, competing spike-generation sites appear to be necessary to account for the experimentally observed effect of occlusion, in which primary afferent firing that results from the simultaneous activation of static and dynamic fusimotor input innervating the same spindle is less than the sum of the rates produced by either fusimotor input individually (Banks et al., 1997, Carr et al., 1998, Fallon et al., 2001). The number of the sites and the exact interaction between them, however, varies across the models. Some suggest two spike-generating sites, one on the bag1 and the other on the bag2 (and chain) (Mileusnic et al. 2006; Schaafsma et al. 1991), while the others suggest three sites, each related to individual intrafusal fiber (Lin and Crago, 2002). The secondary afferent activity results from either pure summation of bag2 and chain activities (Mileusnic et al. 2006) or the competition between them (Lin and Crago, 2002). Some models permit complete occlusion (Lin and Crago, 2002; Schaafsma et al. 1991), while a recent model assumes that the spread of transduction current takes place between the impulse generating sites prior to the competition and suppression, therefore resulting in increased impulse generation at the dominant site (partial occlusion) (Mileusnic et al. 2006). Partial occlusion between two or more spike initiation regions may give rise to irregular firing patterns as a result of collisions between action potentials (Loeb and Marks, 1985) but this has not been modeled explicitly; see Ensemble spindle activity below.

### Spindle motoneuron innervation

All the spindle models addressing the fusimotor innervation assume the existence of single dynamic and single static fusimotor input. While not being widely accepted, there exists limited evidence supporting the existence of two different sources of static fusimotor drives. The direct fusimotor recordings of two static fusimotor efferents firing during decerebrate locomotion in the cat suggest that the two efferent firing profiles are somewhat out of phase; type 2 drive (presumed to be innervating bag2 fiber) leads the type 1 drive (innervating chain fibers) by about 0.17s (Taylor et al., 2000). A recent spindle model (Mileusnic et al. 2006) incorporating the temporal properties of bag2 and chain fibers divided the static fusimotor input into two different static drives. The model more accurately predicted experimentally recorded secondary afferent activity when incorporating the different temporal properties of the different intrafusal fibers. Furthermore, the model findings suggested that the experimentally observed phase differences between the two static fusimotor drives would compensate for the differences in contractile dynamics of their respective intrafusal fibers, resulting in synchronous static fusimotor modulation of the receptors. This provides support for the simplification of including only one type of static fusimotor fiber and drive in spindle models, but it seems likely that they are actually used differentially in some behaviors.

In addition to incorporating fusimotor effects, two spindle models also addressed the effect of beta (β) motoneurons (Maltenfort and Burke 2001; Lin and Crago 2002). The beta (or also called skeleto-fusimotor) motoneurons sometimes innervate both intra- and extrafusal muscle fibers. Various evidence suggests that beta motoneurons receive monosynaptic group Ia excitation comparable to that in motoneurons, which would then result in a positive feedback loop (Burke et al., 1973). While the model of Lin and Crago (2002) suggests that the beta input could be received by the model, this is somewhat unclear from their description. In contrast, the development of Maltenfort’s black box spindle model was largely inspired by the author’s desire to study the effects of beta motoneurons. After designing a spindle model that responds to fusimotor input, the β feedback was simulated as a spindle model receiving a fusimotor drive that is proportional to the instantaneous firing rate of that afferent. The model assumes that the effect of β motoneuron activity on spindle afferent firing has the same strength and time course as the effect for a fusimotor neuron of the same type (static or dynamic), but it is not clear whether or how to account for the general tendency of alpha (extrafusal) motoneurons to fire much slower than gamma (fusimotor) neurons.

### Model parameters and performance

The number of model parameters varies across the models depending on their complexity. When possible, many parameters were derived directly from spindle or extrafusal fiber literature. The remaining free parameters were typically optimized to account for various afferent records under experimentally controlled conditions of kinematics and fusimotor stimulation, mostly from feline hindlimb muscles. Usually no independent data sets were available to validate the models. Their strengths and weaknesses in fitting the available data are discussed for each model below.

## Spindle models

The structural models rely on anatomy and physiology to explain spindle afferent behaviour during different kinematic and fusimotor conditions. The structural spindle models usually use the extensive information available about the internal spindle components and their properties as measured directly in intrafusal fibers, or as inferred indirectly from spindle activity during controlled experiments, or as extrapolated from general properties of striated muscle fibers. Similarly to the biological spindle, the spindle models are driven by a combination of length changes and fusimotor activation (where modelled) in order to generate primary and/or secondary afferent activity. Most models use the extrafusal fiber length as the kinematic input, assuming that it is representative of the length changes being experienced by a spindle lying in parallel to extrafusal fibers. Because of series elasticity in tendons, the velocities of the spindles and extrafusal muscle fibers may differ transiently from the length of the overall musculotendon unit from origin to insertion, which more directly reflects the skeletal posture (Hoffer et al., 1988). Some models address the primary and secondary afferents separately (Mileusnic et al. 2006; Lin and Crago, 2002); others concentrate only on the primary afferent (Rudjord, 1970b; Schaafsma et al. 1991; Hasan, 1983).

### Structural model of Rudjord 1970

Contrary to the initial modelling attempts that used linear transfer functions to describe spindle afferent behaviour, Rudjord pointed out various non-linear transfer properties and developed a second-order model of primary afferents (Figure 2). The model was created to describe the behaviour of the de-efferented spindle, but it included two types of intrafusal muscle fibers (bag and chain) having different mechanical properties. The mechanical model consisted of two parallel mechano-elastic transducers (that simulate the effect of the terminations on the bag and chain fibers) that are in series with the polar region. Therefore, when the length changes are imposed on the entire system, the length changes of the two elastic elements (kλ – chain and kα- bag) are considered to activate two mechano-electric transducers. An amplitude dependent gain of the α-branch transducer is included to account for nonlinear transfer properties prior to superimposing two output components. The two transducers were assumed to have different gains.

The equations of movement for the entire system were developed and a transfer function derived that completely characterized the static and dynamic behaviour of the model. Rather than optimizing the free parameters for specific experimental data, the parameters were adjusted to capture qualitatively various types of published data. By introducing a much higher gain for bag than for chain transducers, the model demonstrated also the ability to capture the stiction overshoot at the initiation of ramp stretches. For sinusoidal stretches of varying frequency, both the amplitude ratio and the phase advance showed only moderate agreement with experimental data (Matthews and Stein, 1969).

Figure 2: Primary afferent model from (Rudjord, 1970).

### Structural model of Hasan 1983

Hasan (1983) presented one of the first models consisting of mathematical elements representing the anatomical components found in the biological spindle. Hasan’s model separates the mechanical filtering process from the later transducing and encoding processes. It describes the mechanical filtering that converts the change in muscular length to change in extension of the nerve endings as a nonlinear process. This is supported by the observation in isolated spindles that increasing amplitudes of sinusoidal stretch produce less than proportional increases in the amplitudes of modulation of spindle tension or receptor potentials of both primary and secondary endings. Transduction and spike encoding were linear, on the assumption that any nonlinearities could be incorporated into the mechanical filtering.

The model differentiates between the sensory and polar region of the intrafusal fibers. The sensory region is modelled as a spring. The polar region is assumed to share many similarities with mechanical features observed in extrafusal fibers and included features of nonlinear velocity dependency and a property akin to friction. Each of the regions (sensory, polar) is described by an equation and the two equations are then summarized in a differential equation. By solving the equation, the length of the sensory region was obtained. Once the equation for length in the sensory region was calculated, the instantaneous discharge rate (afferent activity) was assumed to be a function of the rate and extent of sensory region elongation.

The model also addresses the fusimotor activation but in a qualitative rather than quantitative way. Particularly when integrating fusimotor activation, the model demonstrates a very rigid structure that requires changing its parameters to account for each fusimotor state. It cannot account for the modulated fusimotor activation usually found in cyclic tasks such as locomotion. Like Rudjord’s model, it can account for the initial burst in afferent activity at the start of each stretch, whereas in the biological spindle this occurs only during the first stretch after a period of rest (Matthews, 1972). Nevertheless, the model captures qualitatively many aspects of afferent activity. Because of both its innovative nature and its easy implementation, Hasan’s model attracted the attention of many researchers.

### Structural model of Schaafsma et al. 1991

This structural model of the primary afferent incorporates fusimotor effects that arise from two intrafusal fiber submodels (Figure 3
Figure 3: The structural model of primary afferent from Schaafsma et al, 1991). The figure presents primary afferent response during passive stretch without fusimotor stimulation, with constant dynamic, constant static and combined static and dynamic fusimotor stimulation.
). One submodel represents the bag1 intrafusal fiber and receives dynamic fusimotor input, while the other represents bag2 and chain fibers together and receives static fusimotor input (called bag2 submodel). Both submodels contain a sensory part (the central nuclear bag around which primary afferents terminate) and a muscular part (the contractile portion of a nuclear bag, containing sarcomeres). The sensory parts are modeled as identical linear elastic elements. The muscular part of the bag1 submodel is modelled as a slow twitch muscle fiber, and that of the bag2 as a fast twitch fiber. The forces generated by this part depend on fiber activation, fiber length and contraction velocity (similar to extrafusal muscle fibers). The two submodels differ in terms of the parameters that control the maximal isometric force, the passive elasticity, the maximal force for an eccentric contraction relative to the maximal isometric force, and the passive damping coefficient. The fusimotor input is linearly related to the active force, and scaled with respect to the fusimotor excitation for maximal effect. The model furthermore includes the force enhancement in response to stretch that has been observed in extrafusal muscle fibers.

Each bag submodel generates a depolarization (transducer process) that depends linearly on the length of its sensory region. The firing rate of the bag submodel is linearly dependent on the receptor potential but also on the rate of change of the receptor potential. The final primary afferent firing is obtained by competition between the two action potential trains. The model assumes complete occlusion – the final primary afferent discharge is driven by the intrafusal fiber providing the highest discharge rate.

The bag1 submodel’s muscular part also includes the cross-bridge fixation (or stiction) in order to account for the high sensitivity of muscle spindles to small amplitudes of stretching (e.g. initial burst at the onset of ramp-and-hold stretches). The model assumes one hundred cross-bridge regions that are initially fixed and thereby make the muscular part rigid. If a stretch is introduced, all elongation is initially transmitted to the sensory region resulting in intrafusal fiber force. As the force increases, the cross-bridges gradually release one after the other when their individual breaking forces are reached.

While many parameters used in the model are based on spindle or extrafusal fiber literature, it was necessary to define the remaining ten parameters. The parameters were optimized on experimentally recorded primary afferent activity during 36 ramp-and-hold stretches under 6 different conditions of fusimotor activity (constant static or dynamic fusimotor input, not simultaneous) and 6 different velocities of stretching. The resulting model accounted well for the data used for the fitting process.

The output of the model during small amplitude, 1 Hz sinusoidal stretches was compared to experimental data in terms of amplitude, modulation depth and phase during constant static or dynamic fusimotor input. The model underestimated the afferent activity, generating about one-half the discharge modulation reported in experimental studies. The phase advances predicted by the model tended to be slightly lower than those observed for real spindles.

### Structural model of Lin and Crago 2002

A combined model of extrafusal muscle and spindles was designed in order to study their interaction (Figure 4
Figure 4: The structural spindle model by Lin and Crago (2002). A. A cobmined model of extrafusal muscle and spindles assumes intrafusal length to be proportional to extrafusal fiber length. B. Spindle model consists of three intrafusal fiber models. C. Intrafusal fiber models share the same structure while having different parameters. Each intrafusal fiber has a sensory region in series with a polar region consisting of intrafusal contractile element (CEi), intrafusal fiber mass (Mi), parallel viscous element (Ri) and nonlinear spring corresponding to series elastic element (Kt,i) and passive elastic element (Kp,i).
). A Hill-type muscle model was used to obtain the extrafusal fiber length, which was assumed to be proportional to the intrafusal fiber length. The spindle model is a structural one consisting of three intrafusal fiber models (bag1, bag 2 and chain). The three intrafusal fiber models share the same structure (normalized Hill-type muscle model) while having different parameters (maximal muscle force, optimal muscle fiber length, tendon slack length, muscle mass, etc). Each intrafusal fiber has a sensory region in series with a polar region consisting of intrafusal contractile element (CEi), intrafusal fiber mass (Mi), parallel viscous element (Ri) and nonlinear spring corresponding to series elastic element (Kt,i) and passive elastic element (Kp,i).

Dynamic fusimotor activity excites the bag1 fiber’s polar region, while static fusimotor activity excites bag2 and chain fibers’ polar regions. The model calculates the intrafusal muscle force by multiplication of the activation level, the output of the force–length curve and the output of the force–velocity curve, similarly to extrafusal fibers. The initial burst is modeled as a potentiation of stretch sensitivity (tx) that is a negative (decreasing) exponential function of the activation level. The model assumes that each intrafusal fiber has an independent impulse generator site. The stretch in each intrafusal fiber’s sensory region is calculated to obtain its activity. The primary afferent firing rate is modelled as resulting from the complete occlusion between the bag1 and bag2 fiber activities.

The secondary afferent endings are assumed to be located only on the sensory region rather than juxtaequatorially. The secondary afferent activity results from the complete occlusion among the activities in bag2 and chain fibers. Individual intrafusal fibers also include time constants of activation dynamics (low-pass filter properties). The chain fiber is modelled as having a significantly shorter time constant than slow bag fibers. Overall, 89 parameters, out of which 24 were free parameters, were used in this model.

The model correctly predicted the fractional power dependence on stretch velocity of the primary and secondary ending responses. In the case of the primary afferent, the model performed well during ramp-and-hold and sinusoidal stretches, especially during the stimulation of individual fusimotor inputs, while demonstrating certain limitations in the absence of such inputs. In the case of the ramp-and-hold stretch, the model underestimated the afferent decay time after completion of the stretch, especially in the absence of any fusimotor stimulation (0.2 vs. 0.5s). At the completion of the stretching phase of a ramp-and hold stretch, the model predicted an abrupt and brief cessation of firing, which might be a troubling artefact if the model were to be incorporated into a model of segmental regulation of motor output. Finally, the performance of the model during combined fusimotor stimulation is not described, except for a single ramp-and-hold example presented without matching experimental data. Because the model uses complete rather than partial occlusion, it seems likely that the model will underestimate afferent activity during conditions of simultaneous fusimotor stimulation.

### Structural model of Mileusnic et al. 2006

This model consists of mathematical elements closely related to the anatomical components found in the biological spindle (Figure 5
Figure 5: The structural spindle model by Mileusnic et al. (2006). (a) Spindle model consisting of 3 intrafusal fiber models; (b) intrafusal fiber model consisting of sensory and polar region. Three intrafusal fibers share same structure while having different parameter values.
). It is composed of three nonlinear intrafusal fiber models (bag1, bag2, and chain). Each intrafusal fiber model responds to two inputs: fascicle length and the relevant fusimotor drive (dyanmic fusimotor drive to the bag1 and static fusimotor drive to the bag2 and chain fibers). 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.

All the intrafusal fiber models share the same general structure, a modified version of McMahon’s spindle model (McMahon, 1984). The relative importance of model parameters differs for three intrafusal fiber models to account for the different properties of three fiber types. The intrafusal fiber model consists of two regions, sensory and polar. 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, 1984). The sensory region of each intrafusal fiber model is modelled as a linear elastic element. The model calculates the stretch in the sensory region to obtain the intrafusal fiber’s contribution to the activity of the afferent (prior to taking into account the partial occlusion described below). The polar region, constituting essentially a striated muscle fiber innervated by fusimotor endings, is modelled as a spring 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 polar region’s damping term is primarily modified by dynamic fusimotor stimulation and to a smaller extent by static fusimotor stimulation. Increases in damping result in increases in the velocity sensitivity of the primary afferent, which plays an important role in modulating the spindle’s behaviour during dynamic fusimotor excitation of the bag1 fiber. By contrast, static fusimotor activation produces a small decrease in damping. Active force generation is strongly modulated by the static fusimotor input, which causes a sustained, strong contraction within the bag2 and chain polar regions, producing a stretch in the sensory region and a bias in the afferent activity. Dynamic fusimotor input produces a similar but much weaker effect in the bag1 fiber. The model also incorporates the experimentally observed nonlinear velocity dependency of the spindle’s afferent response (Houk et al. 1981). Furthermore, it incorporates the experimentally observed asymmetric effect of velocity on force production during lengthening and shortening observed in extrafusal striated muscle. The intrafusal fiber models include a biochemical Hill-type equation to capture the saturation effects that take place at high fusimotor stimulation frequencies. Additionally, the models incorporate low-pass filters when converting the fusimotor input into model activation to account for the experimentally measured temporal properties of individual intrafusal fibers, particularly for slow bag intrafusal fibers.

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., 1997, Carr et al., 1998, Fallon et al., 2001) . 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.

The model behaviour is described by a second order differential equation which uses muscle fascicle length and fusimotor inputs (static and dynamic) to calculate the primary and secondary afferent firing. Numerous model parameters were based on the experimental literature for the intrafusal or extrafusal muscle fibers while the remaining parameters were fit using standard least-squares fitting algorithms. 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. The primary and secondary afferent firing during ramp and triangual stretches at different velocities and during different fusimotor stimulation conditions on which model data were optimized is accurately captured (Figure 6
Figure 6: Performance of a spindle model during whole mucle ramp stretches (Mileusnic et al. 2006). Primary afferent response is shown for 3 different whole muscle stretch velocities and 3 different fusimotor conditions (no fusimotor stimulation, constant dynamic fusimotor stimulation, constant static fusimotor stimulation). Secondary afferent response is shown for 3 different whole muscle stretch velocities in absence of fusimotor stimulation.
and 7)
Figure 7: Performance of a spindle model during whole muscle triangular stretches (8mm/s) (Mileusnic et al. 2006). Primary afferent activity was recorded during constant static or dynamic fusimotor stimulation.
. During simultaneous static and dynamic fusimotor stimulation the model correctly captured the primary afferent activity as influenced by the partial occlusion effect (Figure 8)
Figure 8: Performance of a spindle model during sinusoidal stretched and fusimotor stimulation (Mileusnic et al. 2006). Two upper graphs show the primary afferent activity during constant dynamic or static fusimotor stimulation. Two lower graphs show the activity during constant, stimultaneous static and dynamic fusimotor stimulation.
.

The model was validated on an independent data set recorded in the decerebrate cat preparation where two types of naturally occurring static fusimotor drives were recorded (Taylor et al. 2000). Although the existence of two types of static fusimotor drives is still not universally accepted, the researchers suggested 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. Although originally planned to receive one type of static fusimotor drive, the spindle model was modified to account for two static drives. The experimentally recorded kinematic and fusimotor profiles were used to run the modified model and calculate secondary afferent activity. The model simulations were repeated twice, with and without the inclusion of the low-pass filter for bag2 intrafusal fiber model. When subject to two static fusimotor drives, the model more accurately predicted experimentally recorded secondary afferent activity when incorporating the known different temporal properties of the two fiber types. The model suggests that the experimentally observed phase differences between two static fusimotor drives compensate for the differences in contractile dynamics of their respective intrafusal fibers, resulting in synchronous static fusimotor modulation of the receptors, at least in this preparation.

### Black box model of Maltenfort and Burke 2001

The development of the model was largely inspired by the desire to study the effects of beta (β) (or also called skeleto-fusimotor) motoneurons which innervate both intra- and extrafusal muscle fibers (Laporte and Emonet-Dénand, 1976). Some evidence suggests that β motoneurons receive monosynaptic group primary afferent excitation comparable to that in motoneurons, potentially resulting in a positive feedback loop (Burke and Tsairis, 1977).

As proposed previously by several researchers, the model describes the cat’s primary afferent (Ia) as a sum of four components: a pure velocity sensitivity, a pure position sensitivity, a mixed velocity and position sensitivity, and baseline afferent activity at the initial length of the muscle (Lennerstrand and Thoden, 1968b, Lennerstrand and Thoden, 1968a, Hasan, 1983, Prochazka and Gorassini, 1998) (Figure 9
Figure 9: The primary afferent spindle model from Maltenfort and Burke, 2001.
). Spindle afferent response R to muscle length (x) and velocity (v) is summarized by the following equation

$$R = [S_v(\gamma, v) + S_{SS}(\gamma)]x + Q(\gamma, v) + B(\gamma)$$

where $\gamma$ denotes fusimotor drive, $R$ instantaneous primary afferent firing rate, $S_{ss}(\gamma)$ a velocity-independent position sensitivity, $S_v(\gamma,v)$ a velocity-dependent position sensitivity, $Q(\gamma,v)$ a pure velocity sensitivity, and $B(\gamma)$ the baseline firing of the spindle at the initial length.

In their primary afferent spindle model Maltenfort and Burke extended the above mentioned theory and subdivided $B$, $S_{ss}$, $Q$ and $S_v$ into the response of the passive spindle to stretch ($R_{passive}$; no $\gamma$ input) and additive responses to stretch modulated by dynamic ($γ_d$) and static ($\gamma_s$) fusimotor inputs ($ΔR_{dynamic}$, $ΔR_{static}$):

$$R(t) = R_{passive}(t) + f_{occlusion}(\Delta R_{dynamic}, \Delta R_{static}).$$

The passive response is given by

$$R(t) = R_{passive}(t) + f_{occlusion}(\Delta R_{dynamic}, \Delta R_{static}).$$

The portion of spindle response to stretch that is modulated by fusimotor drive is

$$\Delta R_g(t) = [S_{v,g}(\gamma_g, v_f) + S_{ss, g}(\gamma_g)]*x + Q_g(\gamma_g, v_f) + B_g(\gamma_g),$$

where subscript $g$ represents static or dynamic drive, $x$ is muscle mass, $v_f$ estimate of velocity.

The model also includes occlusion:

$$f_{occlusion}(\Delta R_{g1}, \Delta R_{g2}) = \Delta R_{g1} + \Delta R_{g2} - \frac{1}{\frac{1}{\Delta R_{g1}} + \frac{1}{\Delta R_{g2}}},$$

where $ΔR_{g1}$ and $ΔR_{g2}$ are the increments in primary firing rate (over passive) caused by fusimotor inputs $\gamma_1$ and $\gamma_2$ individually.

Finally, baseline firing ($B$), position sensitivity ($S_{ss}$) and velocity-dependent terms ($Q$ and $S_v$) during lengthening and shortening were estimated based on the published experimental records of spindle activity.

The model fairly accurately captured primary afferent activity during ramp stretches at different velocities without fusimotor stimulation as well as with either constant static or dynamic fusimotor stimulation present. During small sinusoidal stretches without fusimotor stimulation as well as during individual or simultaneous static and/or dynamic fusimotor stimulations, the model also performed relatively well for most of the conditions.

In terms of the role of β feedback that inspired the development of the spindle model in the first place, this was simulated as a single spindle receiving a fusimotor drive that was proportional to the instantaneous firing rate of that afferent. The model also assumes that β input has the same influence on primary afferent firing as $\gamma$ motoneurons. In other words, the effect of β motoneuron activity on spindle afferent firing had the same strength and time course as the effect of a $\gamma$ motoneuron of the same type (static or dynamic). During triangular stretching, the β feedback enhanced the primary afferent firing during shortening, while the closure of the β dynamic loop increased primary afferent firing during lengthening. In general, the loop gain increased with velocity and amplitude of stretch but decreased with increased superimposed $\gamma$ rates. The effect of β feedback was most noticeable when β loop and $\gamma$ bias were of different type.

## Ensemble spindle activity

Computing kinaesthesia from spindle afferent activity is not trivial. The firing rates of individual afferents tend to be noisy (Prochazka et al., 1977, Loeb and Duysens, 1979) as a result of transduction noise and mechanical ripple from fusimotor activity. Combining signals from many asynchronously firing afferents results in a Poisson distribution for which signal-to-noise ratio improves with the square root of spike events (Loeb and Marks, 1985). Accurate information about the length and velocity of each muscle can only be extracted by combining activity from all available spindle afferents and only after accounting for fusimotor drive and tendon stretch, as noted above.

Most joints have more than one axis of motion and most muscles cross more than one joint, so information about the length and velocity of many muscles must be combined to build up an accurate sense of posture and motion. The numbers of spindles in almost all human skeletal muscles have been counted (Voss, 1971) and have been used in large-scale musculoskeletal models to test hypotheses about how these signals are combined to account for psychophysical measurements of kinaesthetic resolution at each joint (Scott and Loeb, 1994).

## Golgi tendon organs

Golgi tendon organs (GTO) are tension-sensitive mechanoreceptors found in mammalian skeletal muscles that supply the central nervous system with information regarding active muscle tension by their Ib afferents. The GTOs are most commonly located at the junctions between the muscle fibers and the collagen strands composing tendons and aponeuroses (Figure 10
Figure 10: The structure of Golgi tendon organ. The GTO receptor is located in series between the tendon and muscle fibers that insert into it. A single GTO afferent axon innervates GTO’s complexly interwoven collagen strands. (Figure adapted from Mileusnic and Loeb, 2006).
) (Golgi, 1878, Golgi, 1880). The GTO receptor consists of bundles of collagen fibers that connect small fascicles of muscle to the whole muscle tendon or aponeurosis (Schoultz and Swett, 1972). In other words, the GTO is placed in series between muscle fibers (“muscle end”) and tendon and aponeurosis (“tendon end”), the oppostie of the muscle spindle which lies in parallel with extrafusal muscle fibers. Each GTO receptor is usually innervated by a single large, myelinated Ib afferent that enters the GTO capsule near its equator. The activation of a muscle fiber that inserts onto tendon strands that pass through the GTO straightens some of the collagen strands, compressing and depolarizing the pressure-sensitive afferent endings and eventually resulting in membrane depolarization and initiation of action potentials in the Ib axon (Fukami and Wilkinson, 1977).

## GTO models

To date, only two models of GTO activity have been published (Houk and Henneman, 1967, Mileusnic and Loeb, 2006). The later model was subsequently used to study the ability of the ensemble firing of all GTO receptors in a muscle to accurately encode whole muscle force (Mileusnic and Loeb, 2009).

### Model of Houk and Henneman 1967

Although the main objective of the published work by (Houk and Henneman, 1967) was to test whether active contraction rather than muscle stretching was more effective in activating the GTO, they also offered a linear model of GTO response. In the studies, a laminectomy was performed on adult cats and the dorsal and ventral roots subdivided to record individual GTO activity during stimulation of individual motor units. The firing rate of the Ib afferent $r(t)$ to any input whose time course is specified by $f(t)$ is calculated from the convolution or superposition integral

$$r(t) = \int_0^t f(t-t') h(t') dt'$$

where $h$ represents a model response to a unit impulse. In other words, the convolution integral computes the output by adding up the responses to an infinite series of infinitesimal impulses. The dynamic properties are defined by the means of $h$, where $h$ is defined as

$$h(t) = K(1 + B + C]u_0(t) - K[bBe^{-bt} + c C e^{-ct}]u_{-1}(t)$$

and characterizes the tendon organ. $B$, $b$, $C$, $c$, $K$ are parameters, $u_0(t)$ is the unit impulse function, and $u_{-1}(t-t’)$ is the unit step function the effect of exponential parameters can be best appreciated by looking at the response to a step input

$$r(t) = K[1 + Be^{-bt} + Ce^{-ct}], \;\; t \geq 0$$

The dynamic properties of tendon organs to stimulations of single motor units at different velocities were recorded. Subsequently, these different responses were fitted by the simulated responses of the mathematical model of the receptor where model parameters were adjusted for each individual GTO response. The exponential parameters of the model indicate the magnitude of overshoot ($B$ and $C$) and the decay rate ($b$ and $c$), while $K$ is a gain factor whose physiological correlate is the sensitivity of the receptor (Figure 11
Figure 11: GTO response to a tetanic contraction of a single motor unit (upper curve-slow contraction, lower curve-fast contraction). The input force is represented by solid line while modelled GTO rate of discharge ‘r’ in pulses per second as +. (From Houk and Henneman, 1967)
).

In addition to GTO responses to individual motor unit stimulations, the researchers recorded also GTO activity during multiple motor unit stimulation, which results in non-linear summation. If two motor units are stimulated simultaneously, the resultant response is greater than that which could be obtained from either individually but smaller than that expected from simple addition. While providing an experimental example of non-linear summation, the researchers do not report on model parameters that would capture it.

### Model of Mileusnic and Loeb 2006

A recently developed, anatomically based GTO models attempts to explain the experimentally recorded GTO firing by interactions of various histological components of the GTO. It is based on feline medial gastrocnemius (MG) muscle in which the GTO has been most extensively studied.

A biological GTO consists of bundles of collagen fibers that connect small fascicles of muscle fibers to the whole muscle tendon or aponeurosis (Figure 10). Similarly to a biological GTO, the GTO model receives the tension input from multiple muscle fibers inserting into its capsule. A typical GTO in feline MG is attached to 20 muscle fibers belonging to 13 different motor units (Gregory, 1990).The structure of the GTO model is also strongly influenced by the published literature which suggests that the collagen within the GTO capsule is unevenly packed. The marginal areas of the GTO are typically occupied by densely packed collagen whose fibers run in parallel with one another and only rarely make contact with the GTO afferent axon. The central area of the capsule is occupied largely by loosely packed collagen containing many afferent branches among the complexly interwoven collagen strands (Schoultz and Swett, 1972, Nitatori, 1988).

The GTO model assumes that a single muscle fiber inserting into the capsule interacts with both types of collagen (Figure 12
Figure 12: Structural GTO model (Mileusnic and Loeb, 2006). a. The structure of the GTO model. For simplicity the figure represents only two of the many muscle fibers that insert into the GTO capsule. Each fiber has its collagen divided among bypassing and innervated collagen which is then divided between two common collagen networks having separate impulse generating sites. b. The flower-shaped cross sectional area of collagen within the GTO. The muscle fibers’ collagen was assumed to be arranged in the manner of the flower petals where the fiber’s collagen that occupies the central area of flower-shaped structure belongs to the innervated collagen type (yellow) while the more distally located collagen is of bypassing type (dark grey).
). Because at least two large myelinated branches are typically found within the GTO capsule, each innervating separate portions of the GTO’s loosely packed collagen network, the model assumes the existence of two separate transduction zones, each having its own impulse-generating site. The bypassing collagen attached to an inserting muscle fiber is modelled as uninnervated nonlinear spring, while the remaining innervated collagen is distributed between two common transduction zones. The loosely packed collagen within each transduction zone consists of a nonlinear sensory spring (collagen in direct contact with sensory endings) in series with an element consisting of a nonlinear spring in parallel with a damper that reflects the tendency of collagen strands within the network to gradually rearrange in response to stretch.

The amount of innervated and bypassing collagen that is attached to a single muscle fiber is determined by using a conceptual geometrical arrangement of fibers inserting into the GTO where the GTO’s cross-sectional area is assumed to be flower shaped and the fibers are arranged in the manner of flower petals. The cumulative innervated collagen belonging to all the muscle fibers inserting into the GTO is assumed to occupy an inner zone (10% of the total GTO collagen) of the flower-shaped cross-sectional area. Therefore 90% of the cross-sectional area is assumed to be of the bypassing type. The various extrafusal fiber types (S, FF, FR) have different cross-sectional areas and therefore different amounts of the collagen attached to them, accounting for their different abilities to deform innervated collagen and excite the nerve endings.

The sensory region is modelled as a spring in series with loosely packed collagen and represents the collagen that is in direct contact with the sensory endings. The model assumes that transduction in the sensory endings reflects the amount of stretch of the innervated collagen attached to a given muscle fiber weighted by the cross-sectional area of that collagen. The two separate transduction zones each have their own impulse-generating sites and the net activity of the Ib afferent reflects the dominant site (complete occlusion).

The equations of tensions in various model regions were derived and combined. The model parameters were manually tuned rather than optimized because the available experimental data are sparse and highly variable, particularly for the more complex aspects of GTO behaviour.

The model captures GTO activity well as reported in the literature, including the following complexities:

Dynamic and static responses
After a sudden step activation of the MU, the GTO response consists of a burst (dynamic response) that gradually decays to a constant afferent firing (static response) as a result of the damping term as well as the time-dependencies of force generation in various fiber types. The Figure 13 demonstrates the model’s response to slow, non-fatiguable, and fast fatiguable motor unit activation.
This refers to a phenomenon whereby a GTO’s dynamic response during motor unit activation is decreased after prior activation of the same motor unit (self-adaptation) or a different motor unit (cross-adaptation) (Gregory and Proske, 1979, Gregory et al., 1985). In Figure 13
Figure 13: The GTO model’s ability to capture self- and cross-adaptation (Mileusnic and Loeb 2006).The self-adaptation is given for three types of MUs (S: ··, FR: ― and FF:--). Different amounts of cross-adaptation were obtained by varying the amount of innervated collagen that 2 MUs share. (Mileusnic and Loeb 2006).
, the second responses of the modeled GTO to the same tetanic activation of same (or different) motor unit shows a typical decrease confined to the dynamic response. The experimental literature describes large variability in the amount of cross-adaptation that exists between different MUs (Gregory et al., 1985). Two FF muscle fibers were modeled as sharing either 90% or 10% of the innervated collagen. Figure 13 shows a large or small amount of cross-adaptation, respectively, when the units are activated in succession.
Nonlinear summation
During the activation of multiple MUs the GTO demonstrates nonlinear summation, where two MUs when stimulated simultaneously produce afferent activity that is smaller than the linear summation of the GTO firing rates that individual MUs produce when stimulated independently (Gregory and Proske, 1979) (Figure 14
Figure 14: The GTO model’s ability to capture nonlinear summation (Mileusnic and Loeb, 2006). The GTO afferent activity is obtained in response to progressive recruitment of 8 MUs.
).

### Ensemble GTO activtity

Typically, tens of GTOs are distributed unevenly across the myotendinous junction, where each of them responds vigorously to active tension produced by those muscle fibers that insert into its capsules. The ensemble firing of all GTO receptors in the muscle has been hypothesized to represent a reliable measure of the whole muscle force during normal, size-ordered recruitment of motor units (Jami, 1992) but the precision and accuracy of that information are largely unknown because it is impossible to record activity simultaneously from all GTOs in a muscle.

The GTO model developed by Mileusnic and Loeb (2006) was used subsequently in a larger model of ensemble GTO firing to test Jami’s hypothesis. The model of ensemble GTO firing in a muscle was designed to examine the reliability of the ensemble sensory information (to study the relationship between the aggregate GTO activity and total muscle force). To derive the ensemble GTO activity, the model utilizes several other models. The previously described GTO model (Mileusnic and Loeb, 2006) was employed to obtain afferent activity of each individual GTO receptor in response to tensions of muscle fibers inserting into its capsule. A model of the cross-sectional area of the muscle-tendon junction (including the distributions of various types of muscle fibers, motor units and GTOs) was designed in order to obtain a realistic sample of muscles fibers that insert into each GTO receptor. Finally, the realistic muscle model was used to calculate the tensions of the fibers inserting into individual GTOs.

The relationship between the aggregate GTO activity and whole muscle tension was investigated under conditions of natural muscle recruitment and natural distribution of GTOs and muscle fibers as well as in pathological muscles (for example, reinnervated muscles, artificial recruitment by functional electrical stimulation). The ensemble activity was studied for all conditions during very slow ramp muscle activations as well as during phasic excitation similar to that associated with cyclical behavior such as locomotion.

The model suggests that in the intact muscle where MUs are recruited based on the size principle (Henneman and Mendell, 1981), the ensemble GTO activity accurately encodes force information according to a nonlinear, monotonic relationship that has its steepest slope for low force levels and tends to saturate at the highest force levels. Various types of neuromuscular pathology and electrical recruitment produced different and sometimes non-monotonic input-output relationships, as illustrated in Figure 15
Figure 15: The relationship between the ensemble Golgi tendon organ activity and whole muscle tension (as predicted by ensemble GTO model from Mileusnic and Loeb, 2009). The relationship is shown for four conditions: intact muscle, self-reinnervated muscle, intramuscular electrode stimulation and nerve-cuff electrode stimulation.
.

## References

• Bakker, G J (1980). Histochemical characteristics of muscle spindles in cat dorsal neck muscles. In: Physiology, Vol. M.Sc. Kingston, Ontario, Canada: Queen's University.
• Banks, R; Hulliger, M; Scheepstra, K and Otten, E (1997). Pacemaker activity in a sensory ending with multiple encoding sites: The cat muscle spindle primary ending. The Journal of Physiology 498: 177-199.
• Boyd, I; Gladden, M H; McWilliam, P and Ward, J (1977). Control of dynamic and static nuclear bag fibres and nuclear chain fibres by gamma and beta axons in isolated cat muscle spindels. The Journal of Physiology 265: 133-162.
• Boyd, I A (1976). The response of fast and slow nuclear bag fibres and nuclear chain fibres in isolated cat muscle spindles to fusimotor stimulation, and the effect of intrafusal contraction on the sensory endings. Quarterly Journal of Experimental Physiology and Cognate Medical Sciences 61: 203-254.
• Boyd, I A and Smith, R S (1984). The muscle spindle. In: P J Dyck et al. (Eds.), Peripheral Neuropathy, Vol. 2 (pp. 171-202). W.B. Saunders company.
• Burke, R; Levine, D; Tsairis, P and Zajac, F, III (1973). Physiological types and histochemical profiles in motor units of the cat gastrocnemius. The Journal of Physiology 234: 723.
• Burke, R E and Tsairis, P (1977). Histochemical and physiological profile of a skeletofusimotor (β) unit in cat soleus muscle. Brain Research 129: 341-345.
• Carr, R W; Gregory, J E and Proske, U (1998). Summation of responses of cat muscle spindles to combined static and dynamic fusimotor stimulation. Brain Research 800: 97-104.
• Chen, W and Poppele, R (1978). Small-signal analysis of response of mammalian muscle spindles with fusimotor stimulation and a comparison with large-signal responses. Journal of Neurophysiology 41: 15-27.
• Chin, N K; Cope, M and Pang, M (1962). Number and distribution of spindle capsules in seven hindlimb muscles of the cat. In: D Barker (Ed.), Symposium on Muscle Receptors (pp. 241-248). Hong Kong: Hong Kong University Press.
• Fallon, J B; Carr, R W; Gregory, J E and Proske, U (2001). Summing responses of cat soleus muscle spindles to combined static and dynamic fusimotor stimulation. Brain Research 888: 348-355.
• Fukami, Y and Wilkinson, R S (1977). Responses of isolated golgi tendon organs of the cat. The Journal of Physiology 265: 673-689.
• Golgi, C (1878). Intorno alla distribuzione e terminazione dei nervi nei tendini dell’uomo e di altri vertebrati. Rend R Ist Lomb Sci Lett B11: 445-453.
• Golgi, C (1880). Sui nervi dei tendini dell’uomo et di altri vertebrati e di un nuovo organo nervosa terminale muscolo-tendineo. Mem R Acad Sci Torino 32: 359-385.
• Gregory, J E (1990). Relations between identified tendon organs and motor units in the medial gastrocnemius muscle of the cat. Experimental Brain Research 81: 602-608.
• Gregory, J E; Morgan, D and Proske, U (1985). Site of impulse initiation in tendon organs of cat soleus muscle. Journal of Neurophysiology 54: 1383-1395.
• Gregory, J E and Proske, U (1979). The responses of Golgi tendon organs to stimulation of different combinations of motor units. The Journal of Physiology 295: 251-262.
• Grigg, P; Schaible, H G and Schmidt, R F (1986). Mechanical sensitivity of group III and IV afferents from posterior articular nerve in normal and inflamed cat knee. Journal of Neurophysiology 55: 635-643.
• Hasan, Z (1983). A model of spindle afferent response to muscle stretch. Journal of Neurophysiology 49: 989-1006.
• Henneman, E and Mendell, L M (1981). Functional organization of the motoneuron pool and its inputs. In: V B Brooks (Ed.), Handbook of Physiology Sect I The Nervous System, Vol II, Part 1 (pp. 423-507). Washington, DC: American Physiological Society.
• Hill, D (1968). Tension due to interaction between the sliding filaments in resting striated muscle. The effect of stimulation. The Journal of Physiology 199: 637-684.
• Hoffer, J A; Caputi, A A and Pose, I E (1988). Role of muscle activity on the relationship between muscle spindle length and whole muscle length in the freely walking cat. Proceedings of the Conference on Afferent Control of Posture and Locomotion 1-4: 32.
• Hoffer, J A; Caputi, A A and Pose, I E (1992). Activity of muscle proprioceptors in cat posture and locomotion: Relation to EMG, tendon force and the movement of fibres and aponeurotic segments. IBRO Symposium Series 1-6.
• Houk, J and Henneman, E (1967). Responses of Golgi tendon organs to active contractions of the soleus muscle of the cat. Journal of Neurophysiology 30: 466-481.
• Houk, J C; Rymer, W Z and Crago, P E (1981). Dependence of dynamic response of spindle receptors on muscle length and velocity. Journal of Neurophysiology 46: 143-166.
• Hulliger, M; Matthews P B C and North J (1977). Static and dynamic fusimotor action on the response of Ia fibers to low frequency sinusoidal stretching of widely ranging amplitude. Journal of Physiology 267: 811-838.<(span>
• Jami, L (1992). Golgi tendon organs in mammalian skeletal muscle: Functional properties and central actions. Physiological Reviews 72: 623-666.
• Laporte, Y and Emonet-Dénand, F (1976). The skeleto-fusimotor innervation of cat muscle spindle. In: Progress in Brain Research Understanding the Stretch Reflex, Vol. 44 (pp. 99-105).
• Lennerstrand, G and Thoden, U (1968a). Position and velocity sensitivity of muscle spindles in the cat. II. Dynamic fusimotor single-fibre activation of primary endings. Acta Physiologica Scandinavica 74: 16-29.
• Lennerstrand, G and Thoden, U (1968b). Position and velocity sensitivity of muscle spindles in the cat. III. Static fusimotor single-fibre activation of primary and secondary endings. Acta Physiologica Scandinavica 74: 30-49.
• Lin, C-C K and Crago, P E (2002). Structural Model of the Muscle Spindle. Annals of Biomedical Engineering 30: 68-83.
• Loeb, G E and Duysens, J (1979). Activity patterns in individual hindlimb primary and secondary muscle spindle afferents during normal movements in unrestrained cats. Journal of Neurophysiology 42: 420-440.
• Loeb, G E and Marks, W B (1985). Optimal control principles for sensory transducers. In: I A Boyd and M H Gladden (Eds.), Proceedings of the International Symposium: The Muscle Spindle (pp. 409-415). London: MacMillan Ltd.
• Maltenfort, M G and Burke, R E (2003). Spindle model responsive to mixed fusimotor inputs and testable predictions of β feedback effects. Journal of Neurophysiology 89: 2797-2809.
• Martin, P G; Smith, J L; Butler, J E; Gandevia, S C and Taylor, J L (2006). Fatigue-sensitive afferents inhibit extensor but not flexor motoneurons in humans. The Journal of Neuroscience 26: 4796-4802.
• Matthews, P B C (1972). Mammalian Muscle Receptors and Their Central Actions. London: Edward Arnold Publishing Ltd.
• Matthews, P B C and Stein, R B (1969). The sensitivity of muscle spindle afferents to small sinusoidal changes of length. Journal of Physiology 200: 723-743.
• McMahon, T A (1984). Muscles, Reflexes and Locomotion. Princeton, NJ: Princeton University Press.
• Mileusnic, M and Loeb, G (2009). Force estimation from ensembles of Golgi tendon organs. Journal of Neural Engineering 6: 036001.
• Mileusnic, M P; Brown, I E; Lan, N and Loeb, G E (2006). Mathematical models of proprioceptors. I. Control and transduction in the muscle spindle. Journal of Neurophysiology 96: 1772-1788.
• Mileusnic, M P and Loeb, G E (2006). Mathematical Models of Proprioceptors. II. Structure and Function of the Golgi Tendon Organ. Journal of Neurophysiology 96: 1789-1802.
• Nitatori, T (1988). The fine structure of human Golgi tendon organs as studied by three-dimensional reconstruction. Journal of Neurocytology 17: 27-41.
• Paintal, A S (1960). Functional analysis of group III, afferent fibres of mammalian muscles. Journal of Physiology 152: 250-270.
• Poppele, R E and Bowman, R J (1970). Quantitative description of linear behavior of mammalian muscle spindles. Journal of Neurophysiology 33: 59-72.
• Prochazka, A and Gorassini, M (1998). Models of ensemble firing of muscle spindle afferents recorded during normal locomotion in cats. Journal of Physiology (London) 507(Pt. 1): 277-291.
• Prochazka, A; Westerman, R A and Ziccone, S P (1977). Ia afferent activity during a variety of voluntary movements in the cat. Journal of Physiology (London) 268: 423-448.
• Rudjord, T (1970a). A mechanical model of the secondary endings of mammalian muscle spindles. Kybernetik 7: 122-128.
• Rudjord, T (1970b). A second order mechanical model of muscle spindle primary endings. Kybernetik 6: 205-213.
• Rymer, W Z; Houk, J C and Crago, P E (1979). Mechanisms of the clasp-knife reflex studied in an animal model. Experimental Brain Research 37: 93-113.
• Schaafsma, A; Otten, E and Van Willigen, J D (1991). A muscle spindle model for primary afferent firing based on a simulation of intrafusal mechanical events. Journal of Neurophysiology 65: 1297-1312.
• Schoultz, T W and Swett, J E (1972). The fine structure of the Golgi tendon organ. Journal of Neurocytology 1: 1-25.
• Scott, S H and Loeb, G E (1994). The computation of position sense from spindles in mono-and multiarticular muscles. The Journal of Neuroscience 14: 7529-7540.
• Taylor, A; Durbaba, R; Ellaway, P H and Rawlinson, S (2000). Patterns of fusimotor activity during locomotion in the decerebrate cat deduced from recordings from hindlimb muscle spindles. Journal of Physiology 522: 515-532.
• Voss, H (1971). Tabelle der absoluten und relativen muskelspindelzahlen der menschlichen skelettmuskulatur. Anatomischer Anzeiger 129: 562-572.