An Assessment of Six Muscle Spindle Models for Predicting Sensory Information during Human Wrist Movements

Background

The realization of restoring movement using brain-machine interfaces has begun and there is a great deal more work needed to refine and develop this technology (Wessberg et al., 2000; Nicolelis, 2003; Donoghue et al., 2004; Schwartz, 2004). An area of current development is the incorporation of sensory feedback from proprioceptors, or artificial proprioceptor-like sensors (Clark et al., 2014; Fisher et al., 2014, 2015; McGee et al., 2014; Niu et al., 2014). One of the afferents involved in proprioception is the primary (Ia) afferent innervating the specialized sensory structures in skeletal muscle: muscle spindles (Gandevia, 1996; Prochazka, 1996). It is our belief that studying the sensory coding schemes that have evolved in muscle spindles will be important for developing biomimetic prosthetics.

Information theory has been a useful tool to estimate the capacity of communication channels in engineering and design codes that take advantage of that capacity (Shannon, 1948). Engineering closed-loop neuroprosthetics with proprioceptive sensory encoding would benefit from measuring the information about movement encoded in the firing rate of the muscle spindle afferents. These data can be used to estimate the bandwidth of proprioceptive feedback in the intact sensorimotor system that is being replaced with a neuroprosthetic device. Information theory also provides an estimate of the precision of the neural code, informing the choice of resolution of synthetic sensors (McGee et al., 2014) or stimulus parameters for activating intact sensory systems (Clark et al., 2014; Fisher et al., 2014). To characterize the potential information content it is important to capture the full entropy of the sensory source that would be encountered in ecological conditions. In the visual system investigations of neural coding moved from synthetic input, e.g., oriented bars and sinusoidal gratings, to natural stimuli improving the decoding algorithms (Simoncelli and Olshausen, 2001; Olshausen and Reinagel, 2003). We propose that the same approach should be used to define the information content of the sensory systems underlying proprioception. Prior to collecting experimental data, our intention was to predict the sensory information encoded in muscle spindle firing during movements of the wrist joint in humans. To make these predictions, a plausible model of human muscle spindle firing in response to movement needed to be determined.

Quantitative models based on cat muscle spindle physiology started in the late 1960s (Matthews and Stein, 1969) and continued slowly over the next three decades (Houk et al., 1981; Hasan, 1983; Scott and Loeb, 1994; Prochazka and Gorassini, 1998b; Mileusnic et al., 2006; Niu et al., 2014). To date, the mathematical models of muscle spindle primary afferents have only been tested with data sets acquired in acute and chronic cat experiments. It is imperative that these models be tested against responses from human muscle spindles, especially given suspected differences between species (Prochazka, 1999). Thus a secondary aim of this study is to determine whether mathematical models based on the cat, have predictive value for human muscle spindle responses. We chose to follow the approach of Prochazka and Gorassini (1998b) who compared the performance of six models to muscle spindle data from cat hamstring muscles during locomotion: a behavioral task with natural sensory statistics. In that study there was only a modest improvement from adding an EMG signal to mimic coactivation of gamma and beta motor neurons. Therefore, we have excluded this parameter in our initial study of these models. The six models were assessed in relation to the ensemble firing rate profile of eight human extensor carpi radialis (ECR) muscle spindes (presumed group Ia) recorded during a step-tracking, center-out task. Our emphasis is on a first order assessment of the utility of these models for capturing basic features of the step-tracking task and then using the models to predict the response during a more complex task: random two-dimensional tracking (Paninski et al., 2004).

Methods

Overview

The simulations and models developed for this study are schematically illustrated in Figure 1. The first step was generation of hypothetical wrist movements for two different tasks: center-out and random movements. The kinematics of the wrist movements were simplified to rotations about two orthogonal and intersecting axes. The movements were measured in degrees with respect to a neutral position defined as the origin of a co-ordinate system with flexion/extension movements along the x-axis (extension being positive) and radial/ulnar deviation movements along the y-axis (radial deviation being positive). This Cartesian system was transformed into a polar co-ordinate system. In the polar system pure extension movements pointed toward the right at 0°, pure flexion movements to the left at 180°, radial deviation movements pointed upwards at 90°, and ulnar deviation movements pointed downwards at 270° (Figure 1, left).

FIGURE 1

Figure 1. Schematic of the steps involved in the simulations. All simulations are hypothetical movements of the wrist in right-hand coordinates. The wrist joint coordinate system is illustrated to orient the reader to subsequent polar plot figures where angle is determined by the rotation of the wrist about its two axes. Two types of task are examined: center-out and random movements. fem, flexion/extension movement; rud, radial/ulnar deviation movement; ECRB, extensor carpi radialis brevis muscle.

The resulting rotations of the wrist were input to a model of tendon displacement for the extensor carpi radialis brevis (ECRb) muscle (Loren et al., 1996). The tendon displacement output was then used as input to six different muscle spindle models to generate predicted Ia afferent firing rate. These models have been evaluated in relation to the natural movements during cat locomotion and exhibited good fits to firing rate profiles of ensembles of cat Ia afferents (Prochazka and Gorassini, 1998a,b). For our goal of predicting sensory feedback from Ia afferents during natural human wrist movements, these models seemed a clear choice. A more detailed and comprehensive model was not included in our assessment (Mileusnic et al., 2006). A major goal of this model is the accurate modeling of fusimotor effects, which we did not investigate in this initial assessment. In addition, the length input to this model is fascicle length (in units of optimal muscle fascicle length). Converting from our approximation of tendon displacement to fascicle length was dubious given the sparse experimental data from the human ECR muscle.

All simulations and analysis were performed in Matlab ver 7.0.4 and Simulink ver 6.2 with a time step of 1 ms.

Center-Out Movement Simulations

The center-out task is a staple of sensorimotor neurophysiology and has been used extensively for 2D and 3D reaching experiments in monkeys and humans (Georgopoulos et al., 1982; Gordon et al., 1994). This task, also referred to as step-tracking, has been intensively studied during wrist movements by Hoffman and Strick in both man and monkey (Hoffman and Strick, 1986a,b, 1990, 1993, 1999; Kakei et al., 1999). Our simulations of afferent responses during the center-out task build on the wide-spread use of these movements for studying sensorimotor systems.

The kinematics of the center-out movements were calculated to eight equally spaced targets around a circle in wrist joint space (Figure 1). Minimum jerk trajectories were used as they are descriptive of wrist movements (Stein et al., 1988) and similar to the kinematics predicted by a more accurate causative optimal control model for the wrist (Haruno and Wolpert, 2005). The amplitude of the movements was 7° of joint rotation to match with a subset of data extracted from a previous human muscle spindle study (Jones et al., 2001). The simulations included a short period of 500 ms in the neutral wrist position followed by a movement phase of 650 ms that followed the minimum jerk trajectory (Flash and Hogan, 1985). The simulations continued for a period of 500 ms after reaching the target. The average velocity of movement to the targets was 10.8 deg/s with a peak velocity equal to 20.7 deg/s. These movement velocities correspond to those previously reported for the human data: average 9–10 deg/s with peak velocities of 20–30 deg/s (Jones et al., 2001).

Center-Out Movement—Human Data

Eight ECR muscle spindle afferents were selected from a larger sample of data to compare with the predicted firing rate of the spindle model. This subset of data were selected because: the data were from putative Ia afferents in the ECR muscle; the responses were obtained during a center-out task with kinematics that approximate minimum jerk trajectories; and responses were obtained during movements to all eight targets. In a previous study of the ensemble response of spindle primary afferents in the cat hamstring, the authors reported that after five or six had been averaged to estimate a population response, the addition of more units did not change the population response significantly (Prochazka and Gorassini, 1998b). On this precedent, we felt confident that our data sample was sufficient to estimate general ensemble firing rate behavior of human ECR spindle primary afferents. The spike trains were aligned to movement onset and the average movement duration was 634 ± 19 ms (mean ± SD). To calculate the mean firing rate and 95% confidence intervals for the ensemble of ECR afferents, the individual single-trial spike trains were convolved with a Gaussian kernel then averaged. The width of the kernel function (60 ms) was chosen iteratively considering the phasic response and background discharge rate while comparing the ensemble to the predicted firing rates from the six muscle spindle models.

Random Movement Simulations

While the center-out task has been widely used in sensorimotor neurophysiology to investigate neural coding, there are a number of shortcomings that have been noted (Paninski et al., 2004). First, these movements, as with all point-to-point movements, are relatively straight and exhibit bell-shaped velocity profiles where the peak velocity is proportional to the target distance. These invariant features of movement result in the coupling of position and velocity, two variables that can have independent effects on responses of muscle spindles and proprioception. Another difficulty with the center-out task is that the data are statistically non-stationary, i.e., the mean and variance change over time during the task. The non-stationarity precludes the use of any analysis methods that require a stationary stochastic time series; information-theoretic analysis for example. In addition to the inevitable coupling of position and velocity and the non-stationary data, the center-out task results in movements that occupy a limited amount of the available joint space, and while natural, do not approximate the statistics of ergonomic movements of the wrist. For these reasons, we decided to simulate a random pursuit-tracking task introduced to sensorimotor neural coding studies (Paninski et al., 1999, 2004).

The goal of the random movement simulations was to reproduce the statistics of human wrist movements during typing (Serina et al., 1999). The random movements were generated by a Gaussian random number generator so that the standard deviation of position matched the ergonomic data (Table 1). In order to match the statistics of wrist velocity during typing, a 6th-order low-pass Butterworth filter with cut-off frequency of 1.5 Hz was applied.

TABLE 1

Table 1. Wrist movement statistics.

Given the importance of the simulations of the random movements, we have illustrated the key features in Figure 2. Gaussian curves overlying the histograms (Figures 2A,B) illustrate that a constant position distribution was conserved in the final random wrist movement signals. The curves have the same parameters for a given axis at fast and slow speeds. Power spectral analysis was carried out for position and velocity along each wrist axes to check for a flat spectrum below eleven cut-off frequencies: 0.5–1.5 Hz in 0.1 Hz increments. The flat nature of the power spectral density plots indicated the signals were approximately white within the band-pass of interest and therefore considered random (Figures 2E,F).

FIGURE 2

Figure 2. Features of the random wrist movement simulations. (A,B) The trajectory of the random wrist movements in 2D wrist joint space is illustrated for the fastest (1.5 Hz cut-off) and slowest (0.5 Hz cut-off) movement speeds. The histograms illustrate the distribution of angles about the two wrist axis and a Gaussian curve is overlaid to illustrate the distribution of position is the same at the two different speeds. (C,D) Thirty seconds of the wrist position time series for each wrist axis to illustrate the random nature of the simulated movements. (E,F) Power spectral densities of wrist position in both wrist dimensions for the 1.5 and 0.5 Hz signals respectively.

Tendon Displacement Model

The rotations about the two axes of the wrist were used to calculate the tendon displacement of the ECRb muscle resulting from the movements. These calculations were based on a human cadaveric forearm study that reported equations for instantaneous moment arms with respect to joint rotation (Loren et al., 1996). The equation for ECRb, in the neutral wrist posture, was integrated to generate an equation for instantaneous tendon displacement (in mm) with respect to joint rotation. Tendon displacement and velocity, computed with the transfer function 500 s/(s+500) where s is the Laplace operator, were then used as input to the six muscle spindle models (Figure 1). Corrections for compliance or muscle fiber pennation were not features of the spindle models tested, so were ignored for the present study.

The tendon displacement and velocity during the random wrist movements are illustrated in Figure 3 for the fastest and slowest speeds. The key feature of this figure is that it illustrates the range and distribution of the two inputs to muscle spindles during statistically natural human wrist movements. These data allow comparison of the statistics of tendon movements in humans and that in other species, e.g., cats during locomotion. In the middle column the velocity distribution has been rectified prior to calculating the median value at each of the 11 different filter rates. A Kurskal-Wallis test showed a significant difference in the median tendon velocities as a function of filter rate (p < 0.05). The last column illustrates polar plots where the angle is determined by the position in wrist joint space and the distance to a point is the value of tendon displacement or velocity. The polar plots of displacement illustrate the ECRb tendon is stretched (positive displacement, black) for wrist positions combining flexion and ulnar deviation and shortened for wrist positions in the opposite direction. Tendon velocity is independent of wrist position and can have positive or negative values anywhere in wrist joint space. This illustrates a key objective of the random movement simulations, uncoupling of wrist joint position, and tendon velocity that are unavoidably linked in the center-out task.

FIGURE 3

Figure 3. Features of ECRb tendon displacement and velocity during random human wrist movements. The simulated random movements result in Gaussian distributed inputs to the muscle spindle models. The distribution of tendon displacement is constant across the slow to fast movement speeds, whereas the velocity increases. Since the simulated movements are constrained to the statistics gleaned from ergonomic studies of wrist movements, these stimuli can be considered to have a statistically natural distribution. (A,B) tendon displacement (above) and velocity (below), for fast (1.5 Hz cut-off) and slow (0.5 Hz cut-off) random wrist movements. (C,D) Histograms of tendon displacement and rectified velocity. At the faster movement, the tendon velocity distribution spreads out to higher values. (E,F) Polar plots of tendon displacement and rectified velocity (below) for 1.5 Hz (E) and 0.5 Hz (F). Positive displacement or velocity is black, negative is gray.

Predicted Ia Afferent Firing Rates

Six muscle spindle models were used to generate predicted firing rate time series resulting from the simulated movements (Matthews and Stein, 1969; Chen and Poppele, 1978; Houk et al., 1981; Hasan, 1983; Prochazka and Gorassini, 1998b). All of these models were shown to have some predictive value for estimating ensemble Ia afferent firing rate during normal cat locomotion, though when compared with slow ramp-&-hold stretches one emerged as the most general and accurate. The muscle spindle models were implemented in Simulink, using a previous published block-diagram (Figure 7 in Prochazka and Gorassini, 1998b). As the models were no longer available from the cited internet site, correct implementation was verified by comparing model behavior to the published responses during cat locomotion (Figure 3 in Prochazka and Gorassini, 1998b). The solver used in our simulations was the Euler algorithm with a step-size of 1 ms. For comparison to the human data set, the baseline firing rate of the models was set to 10 imp/s from its original 82 imp/s. This decision was made after averaging the mean firing rate for each of the eight human muscle spindles while holding in the central starting position: 10.2 (SD 1.4) imp/s. Our goal was to make the minimal adjustment to the models rather than leaving parameters free and running a traditional optimization algorithm. This was done to evaluate if the models had predictive value with minimal change, to test the hypothesis that there is little difference between cat and human muscle spindle responses.

Data Analysis

Center-Out

The most wide spread approach to analysis of neurophysiological data during the center-out, or step-tracking, task is the “mean vector” analysis method by which the directional tuning is determined (Georgopoulos et al., 1982; Gribble and Scott, 2002). This approach has been extensively used for analysis of muscle spindle coding during passive movements of the ankle (Bergenheim et al., 2000; Roll et al., 2000, 2004; Ribot-Ciscar et al., 2002, 2003) as well as our previous study of active and passive movements of the wrist (Jones et al., 2001). The length of the mean vector was normalized to values between 0 and 1. To test for significance of directional tuning, a bootstrap test was used where the mean rate and target angle are resampled from the original data (4000 resampling trials) and the resulting resampled vector length is compared to the original. If fewer than 200 resampled vectors are longer than the original, the directional tuning is considered significant with 95% confidence. Note, resampling is not done if the original mean vector length is <0.001, as this already indicates a non-directional distribution in circular statistics (Batschelet, 1981).

Quantitative Comparison of Human Data and Model

Both qualitative and quantitative comparisons were done to evaluate the ability of the six muscle spindle models to capture the dynamics of the human experimental data during the center-out task. Root mean square (RMS) error was calculated between each model and the ensemble data on the following measures: mean firing rate, directional tuning vector length, static index, dynamic index, and temporal similarity. For the first two measures, an RMS error was calculated during two phases of the task: (1) movement to the target, and (2) holding on the target. The RMS error for a particular measure was normalized by dividing by the median RMS error for the six models. The normalized errors were summed over all the measures to give a total error score. The error score has a lower limit of zero, which corresponds to perfect prediction, and an unbounded upper limit. The consolidated error score was used to rank order the six models in their ability to predict the human ensemble data. Most of these measures are self-explanatory with full details given in Malik (2006). The measure of temporal similarity evaluated the difference between the ensemble instantaneous firing rates and the firing rates predicted by the models over intervals of 100 ms. The squared difference in firing rate at each interval was summed over the duration of a trial to each of the eight targets. The final RMS error for each model was computed by dividing the sum of squared differences by the number of 100 ms bins and number of targets, then taking the square root.

Random Movements

The random movement data was analyzed using entropy and mutual information for a continuous channel (Shannon, 1948). This approach to analysis of sensory stimuli in neural systems has proven to be a powerful framework for nonparametric and nonlinear analysis (e.g., Bialek et al., 1991; Bialek and Rieke, 1992). Following Shannon we adopt the logic that information = uncertainty = entropy and, using base 2 for our logarithms, interpret the entropy as the smallest number of bits needed to communicate our signal of interest, which is tendon displacement or velocity. The ability to transmit these signals that convey the mechanical state of the muscle using muscle spindle firing rate was assessed by mutual information. This quantity measures the average number of bits an observer receives about the tendon displacement or velocity by observing firing rate. The probability density function for a signal was estimated from a histogram of the signal after which entropy was calculated. Mutual information is a linear addition of the respective entropies, I_xy = H_x + H_xy, where H is entropy, x is the input signal (tendon displacement or velocity) and y is the firing rate of one of the six muscle spindle models. The relevant equations can be found in Shannon (1948), and have been reproduced in other sources (e.g., Cover and Thomas, 1991; Rieke et al., 1999). Calculations were done in Matlab using (Moddemeijer, 1989).

Results

Center-Out Movements

In this section we assess whether a group of six muscle spindle models, which have been developed and used to interpret data from cat studies, give a general and accurate prediction of the responses of human muscle spindles during a step-tracking task.

Comparing Time Series Data

In the human data set that was analyzed for this study, the subjects were asked to move quickly to the targets with an emphasis on accuracy. The main movement phase to the target was in a straight line with a bell-shaped velocity. Movements to some targets, for which the ECR muscle is the primary agonist, resulted in a decrease in firing rate of the human muscle spindles as illustrated in Figure 4. During movement to target 2 (negative tendon displacement, Figure 3E), this muscle spindle was temporarily unloaded resulting in a pause in the spike train (Figure 4A, left upper trace). A continuous estimate of firing rate generated from the ensemble firing rates of all eight spindles in the data set showed a transient period of decreased firing rate during movement followed by recovery of firing rate during the static hold-target phase (Figure 4A, left middle trace). Movement in the opposite direction to target 6 (positive tendon displacement, Figure 3E) resulted in a burst of action potentials from the muscle spindle afferent with instantaneous firing rates reaching almost 90 imp/s followed by adaptation to a lower rate during the static phase. The continuous firing rate estimate from the ensemble showed a clear separation of dynamic and static firing rate responses during “movement to” and subsequent “holding on” the target.

FIGURE 4

Figure 4. Comparing movements to opposite targets. Upper trace shows an example spike train from a human ECR muscle spindle during movements to two targets. Middle trace shows the continuous firing rate generated from the ensemble data of 8 ECR spindle afferents. The first dashed line indicates the onset of movement to the targets and the second is 650 ms later, the end of the minimum jerk movement phase.

The changes in firing rate predicted by the models during movements to these same targets had similar dynamic and static phases. During movements to target 2, five of the six models showed a transient decrease in firing rate during movement with recovery of firing rate during the static phase of holding on the target. Simulations of movements to target 6 showed a transient increase in predicted firing rate that was higher than the rate after reaching the target. The firing rates predicted by four of the models during the dynamic phase were of similar magnitude to the human data, but only three of these predicted rates similar to the human example in the static phase.

While the individual movements to these two targets were nearly 100 ms shorter than the simulated movements (650 ms), the grouped data for all eight spindles during movements to all targets was 634 ± 19 ms (mean ± SD). Therefore, the durations in the experimental data and simulations are of a similar magnitude to allow qualitative comparison. To check how well the models predicted the temporal behavior and variability during this task, the responses of a single spindle during four repeated movements to the same target were averaged (Figure 5). The data were aligned to movement onset and the average firing rate with 95% confidence intervals calculated (Figure 5, gray band). The human spindle showed a clear dynamic phase where the average maximum was 18 imp/s, up from a baseline of 8.5 imp/s, and a static component where average rate was 12 imp/s.

FIGURE 5

Figure 5. Comparing the temporal dynamics of firing rate and variability for a human muscle spindle and six models. Continuous firing rate responses from a single human ECR spindle, recorded during four repeat movements to the target at 270°. A smoothed response was generated from the spike trains using a Gaussian kernel (60 ms bandwidth) and the five repeated movements were averaged to find the mean response and 95% confidence intervals (gray band). The movement duration was 662 ± 95 ms (mean, SD). Overlaid are predicted firing rates from the six models for a minimum jerk trajectory (650 ms duration) to the same target. Baseline firing rate in the models was set to 10 imp/s and the average for the human data was 8.5 imp/s. Records are aligned at movement onset (0.5 s).

Since these movements are well described by a minimum jerk trajectory, we calculated the minimum jerk trajectory to the same target, with movement duration equal to 650 ms (Figure 5). Three of the models showed qualitative features similar to the human data (Figure 5). The remaining three models did not make good predictions of the static phase activity: one was too strong (orange) while the others were too weak (yellow and green) compared to the human data.

These results suggest that some of the models do a better job at predicting the temporal profile of human muscle spindle firing rates during this task. While we have not attempted to match the human kinematics in detail, the predicted minimum jerk trajectories result in predicted firing rates that have clear dynamic and static phases that are similar in magnitude to the average human data. Next we examine whether the models can predict firing rates during movements to all targets in the center-out task.

Comparing Directional Tuning

The human ECR spindle data were analyzed to measure the mean firing rate during the dynamic and static phases of movement to each of the eight targets. Six of the eight afferents were directionally tuned during the dynamic (move) phase (bootstrap, p < 0.05) while only three were significantly tuned during the static (hold) phase. The data were then used to estimate the directional tuning for the population. The mean vector for the dynamic phase had a normalized length of 0.28 and an angle of 239°. The mean vector for the static phase was 0.11 long with an angle of 240°. A bootstrap test showed that both vectors were significantly tuned (p < 0.05). This result differs from the previous study in which these units were part of a larger sample of ECR spindle afferents. In the previous study, the data were not significantly tuned during the static phase (Jones et al., 2001).

All the models showed directional tuning during the dynamic phase with a mean vector angle of 225°. The lengths of the mean vectors are given in the legend for Figure 6 where the directional tuning is illustrated in polar coordinates. The central gray area illustrates the distribution of the human data and the black arrow the direction of the mean vector. The accuracy of the mean vector predicted by the models (Figure 6 red arrow) was surprising given the variability in the human data and the simplifying assumptions for the biomechanical wrist model and simulations. The human data shows higher firing rates than those predicted by the models in the direction opposite the preferred direction, which could be due to co-activation of gamma or beta motor neurons not accounted for in these simulations.

FIGURE 6

Figure 6. Comparing directional tuning in polar plots. The top panel (A) shows mean firing rate for the models (see legend) and averaged human data (gray area). The preferred direction of the models (computed by mean vector) was 225° (red arrow) and the preferred direction of the human data was 239° (black arrow). The length of the normalized mean vectors (not illustrated) were: 0.28 (human), 0.28 (red), 0.46 (orange), 0.62 (yellow), 0.45 (green), 0.52 (blue), and 0.46 (purple). The bottom panel (B) shows mean firing rates during the static phase of holding on the target. Two models (yellow and green) were not significantly tuned during the hold phase (mean vector length = 0.0), while the others had the same preferred direction (225°) and normalized mean vector lengths of: 0.25 (red), 0.64 (orange), 0.28 (blue), and 0.28 (purple). The mean vector for the human data had an angle of 240° and a length of 0.11, which was significant (p < 0.05).

The data in Figures 4–6 provide qualitative evaluation of the ability of the models to match the human ensemble data. RMS errors were calculated for seven different measures during the center-out task during movements to, and holding on, the eight targets to provide a quantitative comparison of the models. The RMS errors were normalized and summed to produce a total error score, which is a relative measure of the models against each other. These data are presented in Table 2 and indicate that the three models with the lowest error scores in rank order are: red, Chen and Poppele (1978); purple, Hasan (1983); and blue, Prochazka and Gorassini 2 (Prochazka and Gorassini, 1998b). Therefore, these models provide the best overall prediction of the ensemble human muscle spindle data.

TABLE 2

Table 2. Overall error for model ranking.

Random Movements

In the previous section we found that previous models of cat muscle spindles could approximate firing rate and directional tuning of a sample of human ECR muscle spindles. The three models with the lowest overall error scores were used to estimate how much information muscle spindle firing rate carries about movement. This was achieved by simulating random movements that dissociate joint position and velocity to maximize the entropy of the stimulus while remaining an ecologically valid input.

Prediction of the Firing Rate Time Series

The time series plots of firing rate predicted by the three spindle models during simulated random movements at two speeds are shown in Figure 7A. The predicted firing rates during the fast movements had a striking “peaky” structure that showed sharp rises to a peak firing rate followed by periods of low to zero firing rates. The primary differences predicted by the three models are the dynamic gain and frequency response with the Chen and Popele (red) model having the lowest dynamic gain and frequency response. The mean predicted firing rates, over 5 min of simulated movement, for the three models ranged from 10.5 to 13.1 imp/s with peak firing rates <50 imp/s (Table 3). Using the mean firing rate to characterize the different models was complicated by periods of zero firing rate. The models predicted periods of unloading, or silencing, that ranged between 12.5 and 35.0% of movement duration. This had a significant effect on the probability distributions for the firing rate signals (Figure 7B). The probability distribution of firing for the Chen and Popele model is unimodal (if the peak at zero is ignored) while the other two models predict a bimodal distribution of firing rates. The differences in distribution are associated with differences in the dynamic gain and frequency response illustrated in the time series plots.

FIGURE 7

Figure 7. The time series and histograms of firing rates predicted for random movements. The random 2D movements in wrist joint space resulted in changes in ECRb tendon length that were input to the muscle spindle models to predict firing rates during movements of the human wrist encountered in normal movements. (A). During fast movements, the predicted firing rates are <40 imp/s. The main distinguishing feature when comparing the three models is the peak rate and periods of silence (i.e., firing rate of zero). Even during slow movements the Prochazka and Gorassini 2 model predicts periods of zero firing rate. (B). The histograms all have the same scale starting at zero firing rate and the dotted line indicates 20 imp/s. The means of these distributions, with and without using zero firing rate in the calculations, are given in Table 3. The shapes of firing rate histograms predicted by the three models are different: unimodal and bimodal. Experimental data during a similar task should be plotted in a similar fashion to evaluate the model predictions.

TABLE 3

Table 3. Descriptive statistics of firing rates during fast random movements.

Assessing the Representation of Tendon Displacement and Velocity Using Entropy and Mutual Information

The final goal of this investigation was to characterize the muscle spindle responses to sensory stimuli using information-theoretic analysis. The amount of information in the sensory stimuli was quantified by the entropy of tendon displacement and velocity at each of the eleven random movement speeds. The average entropy and standard error for tendon displacement was 3.05 ± 0.01 bits per symbol across all eleven speeds while the entropy of the tendon velocity signal increased linearly from 3.88 ± 0.01 at the slowest speed to 5.49 ± 0.01 bits per symbol at the fastest speed.

The muscle spindle is a physiological transducer and the mathematical representation of transduction states that the entropy of the firing rates should be less than or equal to the entropy of the source they are encoding (Theorem 7 in Shannon, 1948). The entropies of the firing rates predicted by the models at the fastest speed (Table 3) are all greater than the entropy for tendon displacement and less than the entropy for the tendon velocity signal. Theoretically, if the models were perfect inverse transducers of displacement and velocity then the entropy of the firing rates would be equal to the sum of their source entropies, in this case 3.05+5.49 = 8.54. Instead the firing rate entropies are about half this value. But these measures of firing rate entropy do not quantify how much information firing rates encode about the stimulus.

To measure the transmission of tendon displacement or velocity information in the firing rates of the models we calculated Shannon's mutual information. The mutual information values were very low between firing rate and tendon displacement (I_{rate, disp} in Table 3). On average about 3.3% of the information about tendon displacement is preserved in the firing rates of the models. In comparison about 22.8% of the information about tendon velocity is preserved in the firing rates of the models during the fastest movements (I_{rate, vel} in Table 3).

Discussion

Muscle spindles are an important source of sensory feedback that signal information about the position and movement of joints (Stein et al., 2004). Over the past three decades recordings have been made from single muscle spindle afferents in behaving cats, monkeys and humans. The experiments are technically difficult so valid models are helpful to consolidate understanding and guide further experiments. Until now, none of the six models compared here had been tested against a human data set. Our goal in this report was to compare six models against human data and then predict responses expected from human muscle spindles during a novel experimental paradigm: random pursuit tracking.

In simulations of center-out movements, we found that the ensemble firing rate and directional tuning of human ECR Ia afferents was adequately captured by three of the six models: Chen and Poppele (red), Prochazka and Gorassini 2 (P&G2, blue) and Hasan (purple). The three remaining models were considered a poorer representation of the ensemble given the error scores (Table 2).

Limitations

We have used minimum jerk trajectories rather than fitting to the actual kinematics; we have excluded any terms that would capture phasic gamma or beta motor neuron activity; we have excluded the finer details of the muscle geometry and biomechanics by using tendon displacement as the main stimulus input to the models; and, we have used models that output firing rate as a continuous function rather than a series of action potentials. Despite these simplifications, as a first order approximation these models have clear predictive value for estimating human ECR muscle spindle response during wrist movements. One of the differences we noted between the models and human data was the prediction of zero firing rates during movements opposite to the preferred direction. The human spindle response during movement to target 2 (Figure 4A) was atypical for the sample of eight spindles reported here. The majority of human ECR spindles (75%) did not fall silent during movement to targets opposite to the preferred direction. In comparison, the majority of the models predict a period of zero firing rate. This difference between model predictions and human data could be improved by adding an EMG-linked gamma motor neuron activation term. This should be done in conjunction with human experiments that include a load at the wrist since the “non-forceful” contractions in these step-tracking data do not recruit the ECR muscle strongly. We hypothesize that by creating a load that would recruit the ECR muscle during this task, responses would switch from “hamstring-like” to “triceps surae-like.” The triceps surae muscles are more strongly recruited during locomotion and required an EMG-linked term to improve the prediction of the models (Prochazka and Gorassini, 1998b).

After convincing ourselves that the models had predictive value for human spindles, we simulated the random movement task. We had a number of objectives for pursuing these simulations. First, we wanted to address the issue of whether firing rates approaching those reported during cat locomotion would be predicted during human movements that were ecologically valid. The human experimental data has been criticized for being unnaturally slow, and this has been used as a possible explanation for the lower overall firing rates (Prochazka, 1999). Second, we wanted to use a task that could predict tuning properties of muscle spindles when joint position and velocity were dissociated. Finally, we wanted to use a task that generated data amenable to analysis by statistical methods other than the “mean vector” directional tuning approach. These objectives could have been achieved in a number of ways, but we chose to match the statistics of the movements to the ergonomics of typing.

Importance of Natural Movement Statistics?

What is the appropriate stimulus ensemble for evaluating the sensory information encoded in muscle spindle responses? If you subscribe to the notion that sensory neurons have adapted and evolved in an environment where some stimuli are more likely than others, then the answer is found in the statistics of that natural environment. Over the past decade, it has become clear that sensory encoding in the visual and auditory systems is specifically adapted to the statistics of natural environments (Simoncelli and Olshausen, 2001; Olshausen and Reinagel, 2003). In the primary visual system for example, it has been shown that models based on synthetic stimuli do not generalize to natural stimulus statistics (David et al., 2004). Based on these findings from other areas of sensory neurophysiology, we wanted to test the muscle spindle models with natural stimulus ensembles—in a statistical sense.

The random movement simulations resulted in a number of predictions that are amenable to testing. While we have ranked the models according to similarity of responses to human data in the center-out task, the firing rates predicted during random movements will allow more definitive ranking of the models against human data. The information-theoretic analysis was useful for estimating a lower boundary of information transmission by muscle spindles. Since the models output continuous functions, that are filters of the original spike train, the entropy in the firing rates is less than that of the original spike trains. Less entropy means a lower capacity for transmitting information about the stimulus.

Conclusions

We have found that models based on cat muscle spindle primary afferents have some predictive value for predicting the dynamic and static features of human muscle spindle firing during a simple task. The remaining discrepancy between cat and human data is the lower overall mean firing rate in humans. These models were then used to predict the information content of natural wrist movements and the mutual information between kinematics and spindle firing. This analysis predicts that the information that is represented in spindle firing rate is more strongly weighted to velocity when compared to length. This prediction, which needs to be tested empirically, suggests that biomimetic neuroprosthetic systems should concentrate more on providing velocity feedback to central structures responsible for decoding muscle spindle signals.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The primary support for this work was a grant from The Whitaker Foundation with additional support from the Alberta Heritage Foundation for Medical Research (AHFMR) and the Canadian Institutes of Health Research (CIHR). PM received a scholarship from CIHR-INMHA. We are grateful to Lora Major and Neil Tyreman who provided assistance with preparation of figures. We are also indebted to Professor Prochazka for the implicit and explicit mentoring.

References

Batschelet, E. (1981). Circular Statistics in Biology. London; New York, NY: Academic Press.

Google Scholar

Bergenheim, M., Ribot-Ciscar, E., and Roll, J. P. (2000). Proprioceptive population coding of two-dimensional limb movements in humans: I. muscle spindle feedback during spatially oriented movements. Exp. Brain Res. 134, 301–310. doi: 10.1007/s002210000471

PubMed Abstract | CrossRef Full Text | Google Scholar

Bialek, W., and Rieke, F. (1992). Reliability and information transmission in spiking neurons. Trends Neurosci. 15, 428–434. doi: 10.1016/0166-2236(92)90005-S

PubMed Abstract | CrossRef Full Text | Google Scholar

Bialek, W., Rieke, F., de Ruyter van Steveninck, R. R., and Warland, D. (1991). Reading a neural code. Science 252, 1854–1857. doi: 10.1126/science.2063199

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W. J., and Poppele, R. E. (1978). Small-signal analysis of response of mammalian muscle spindles with fusimotor stimulation and a comparison with large-signal responses. J. Neurophysiol. 41, 15–27.

PubMed Abstract | Google Scholar

Clark, G. A., Wendelken, S., Page, D. M., Davis, T., Wark, H. C., Normann, R. A., et al. (2014). “Using multiple high-count electrode arrays in human median and ulnar nerves to restore sensorimotor function after previous transradial amputation of the hand,” in 36Th Annual International Conference Of The IEEE Engineering In Medicine & Biology Society. doi: 10.1109/EMBC.2014.6944001

PubMed Abstract | CrossRef Full Text | Google Scholar

Cover, T. M., and Thomas, J. A. (1991). Elements of Information Theory. New York, NY: Wiley.

Google Scholar

David, S. V., Vinje, W. E., and Gallant, J. L. (2004). Natural stimulus statistics alter the receptive field structure of v1 neurons. J. Neurosci. 24, 6991–7006. doi: 10.1523/JNEUROSCI.1422-04.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Donoghue, J. P., Nurmikko, A., Friehs, G., and Black, M. (2004). Development of neuromotor prostheses for humans. Suppl. Clin. Neurophysiol. 57, 592–606. doi: 10.1016/S1567-424X(09)70399-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Fisher, L. E., Ayers, C. A., Ciollaro, M., Ventura, V., Weber, D. J., and Gaunt, R. A. (2014). Chronic recruitment of primary afferent neurons by microstimulation in the feline dorsal root ganglia. J. Neural Eng. 11:036007 doi: 10.1088/1741-2560/11/3/036007

PubMed Abstract | CrossRef Full Text | Google Scholar

Fisher, L. E., Hokanson, J. A., and Weber, D. J. (2015). “Neuroprostheses for somatosensory function,” in Implantable Neuroprostheses for Restoring Function, ed K. Kilgore (Cambridge, MA: Elsevier), 127–151.

Flash, T., and Hogan, N. (1985). The coordination of arm movements: an experimentally confirmed mathematical model. J. Neurosci. 5, 1688–1703.

PubMed Abstract | Google Scholar

Gandevia, S. C. (1996). “Kinesthesia: roles for afferent signals and motor commands,” in Handbook of Physiology, Exercise: Regulation and Integration of Multiple Systems, eds L. B. Rowell and J. T. Shepherd (New York, NY: Oxford University Press), 89–127.

Georgopoulos, A. P., Kalaska, J. F., Caminiti, R., and Massey, J. T. (1982). On the relations between the direction of two-dimensional arm movements and cell discharge in primate motor cortex. J. Neurosci. 2, 1527–1537.

PubMed Abstract | Google Scholar

Gordon, J., Ghilardi, M. F., and Ghez, C. (1994). Accuracy of planar reaching movements. I. Independence of direction and extent variability. Exp. Brain Res. 99, 97–111. doi: 10.1007/BF00241415

PubMed Abstract | CrossRef Full Text | Google Scholar

Gribble, P. L., and Scott, S. H. (2002). Method for assessing directional characteristics of non-uniformly sampled neural activity. J. Neurosci. Methods 113, 187–197. doi: 10.1016/S0165-0270(01)00492-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Haruno, M., and Wolpert, D. M. (2005). Optimal control of redundant muscles in step-tracking wrist movements. J. Neurophysiol. 94, 4244–4255. doi: 10.1152/jn.00404.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Hasan, Z. (1983). A model of spindle afferent response to muscle stretch. J. Neurophysiol. 49, 989–1006.

PubMed Abstract | Google Scholar

Hoffman, D. S., and Strick, P. L. (1986a). Activity of wrist muscles during step-tracking movements in different directions. Brain Res. 367, 287–291. doi: 10.1016/0006-8993(86)91604-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoffman, D. S., and Strick, P. L. (1986b). Step-tracking movements of the wrist in humans. I. Kinematic analysis. J. Neurosci. 6, 3309–3318.

PubMed Abstract | Google Scholar

Hoffman, D. S., and Strick, P. L. (1990). Step-tracking movements of the wrist in humans. II. EMG analysis. J. Neurosci. 10, 142–152.

PubMed Abstract | Google Scholar

Hoffman, D. S., and Strick, P. L. (1993). Step-tracking movements of the wrist. III. Influence of changes in load on patterns of muscle activity. J. Neurosci. 13, 5212–5227.

PubMed Abstract | Google Scholar

Hoffman, D. S., and Strick, P. L. (1999). Step-tracking movements of the wrist. IV. Muscle activity associated with movements in different directions. J. Neurophysiol. 81, 319–333.

PubMed Abstract | Google Scholar

Houk, J. C., Rymer, W. Z., and Crago, P. E. (1981). Dependence of dynamic response of spindle receptors on muscle length and velocity. J. Neurophysiol. 46, 143–166.

PubMed Abstract | Google Scholar

Jones, K. E., Wessberg, J., and Vallbo, A. B. (2001). Directional tuning of human forearm muscle afferents during voluntary wrist movements. J. Physiol. 536, 635–647. doi: 10.1111/j.1469-7793.2001.0635c.xd

PubMed Abstract | CrossRef Full Text | Google Scholar

Kakei, S., Hoffman, D. S., and Strick, P. L. (1999). Muscle and movement representations in the primary motor cortex. Science 285, 2136–2139. doi: 10.1126/science.285.5436.2136

PubMed Abstract | CrossRef Full Text | Google Scholar

Loren, G. J., Shoemaker, S. D., Burkholder, T. J., Jacobson, M. D., Fridén, J., and Lieber, R. L. (1996). Human wrist motors: biomechanical design and application to tendon transfers. J. Biomech. 29, 331–342. doi: 10.1016/0021-9290(95)00055-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Malik, P. (2006). Muscle Spindle Ia Afferents: Cat Models Applied to Human Data. Department of Biomedical Engineering, University of Alberta Library, University of Alberta.

Matthews, P. B., and Stein, R. B. (1969). The regularity of primary and secondary muscle spindle afferent discharges. J. Physiol. 202, 59–82. doi: 10.1113/jphysiol.1969.sp008795

PubMed Abstract | CrossRef Full Text | Google Scholar

McGee, T. G., Para, M. P., Katyal, K. D., and Johannes, M. S. (2014). “Demonstration of force feedback control on the Modular Prosthetic Limb,” in 2014 IEEE International Conference on Systems, Man & Cybernetics (SMC) (San Diego, CA), 2833–2836. doi: 10.1109/SMC.2014.6974358

CrossRef Full Text | Google Scholar

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. J. Neurophysiol. 96, 1772–1788. doi: 10.1152/jn.00868.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Moddemeijer, R. (1989). On estimation of entropy and mutual information of continuous distributions. Signal Process. 16, 233–248. doi: 10.1016/0165-1684(89)90132-1

CrossRef Full Text | Google Scholar

Nicolelis, M. A. (2003). Brain-machine interfaces to restore motor function and probe neural circuits. Nat. Rev. Neurosci. 4, 417–422. doi: 10.1038/nrn1105

PubMed Abstract | CrossRef Full Text | Google Scholar

Niu, C. M., Nandyala, S. K., and Sanger, T. D. (2014). Emulated muscle spindle and spiking afferents validates VLSI neuromorphic hardware as a testbed for sensorimotor function and disease. Front. Comput. Neurosci. 8:141. doi: 10.3389/fncom.2014.00141

PubMed Abstract | CrossRef Full Text | Google Scholar

Olshausen, B., and Reinagel, P. (2003). Sensory coding in the natural environment. Network 14, 3:i–3:ii. doi: 10.1088/0954-898X_14_3_001

CrossRef Full Text | Google Scholar

Paninski, L., Fellows, M. R., Hatsopoulos, N. G., and Donoghue, J. P. (2004). Spatiotemporal tuning of motor cortical neurons for hand position and velocity. J. Neurophysiol. 91, 515–532. doi: 10.1152/jn.00587.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Paninski, L., Fellows, M. R., Hatsopoulos, N. G., and Donoghue, J. P. (1999). Coding Dynamic Variables in Populations of Motor Cortex Neurons. Miami Beach, FL: Society for Neuroscience. Program No. 665.669. Neuroscience Abstracts.

Prochazka, A. (1996). “Proprioceptive feedback and movement regulation,” in Handbook of Physiology: Exercise: Regulation and Integration of Multiple Systems eds L. B. Rowell and J. T. Shepherd (New York, NY: Oxford University Press), 89–127.

Prochazka, A. (1999). Quantifying proprioception. Prog. Brain Res. 123, 133–142. doi: 10.1016/S0079-6123(08)62850-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Prochazka, A., and Gorassini, M. (1998a). Ensemble firing of muscle afferents recorded during normal locomotion in cats. J. Physiol. 507(Pt 1), 293–304. doi: 10.1111/j.1469-7793.1998.293bu.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Prochazka, A., and Gorassini, M. (1998b). Models of ensemble firing of muscle spindle afferents recorded during normal locomotion in cats. J. Physiol. 507(Pt 1), 277–291. doi: 10.1111/j.1469-7793.1998.277bu.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ribot-Ciscar, E., Bergenheim, M., Albert, F., and Roll, J. P. (2003). Proprioceptive population coding of limb position in humans. Exp. Brain. Res. 149, 512–519. doi: 10.1007/s00221-003-1384-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ribot-Ciscar, E., Bergenheim, M., and Roll, J. P. (2002). The preferred sensory direction of muscle spindle primary endings influences the velocity coding of two-dimensional limb movements in humans. Exp. Brain Res. 145, 429–436. doi: 10.1007/s00221-002-1135-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Rieke, F., Bialek, W., Warland, D., and de Ruyter van Steveninck, R. (1999). Spikes: Exploring the Neural Code. Cambridge; Mass: MIT Press.

Roll, J. P., Albert, F., Ribot-Ciscar, E., and Bergenheim, M. (2004). Proprioceptive signature of cursive writing in humans: a multi-population coding. Exp. Brain Res. 157, 359–368. doi: 10.1007/s00221-004-1853-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Roll, J. P., Bergenheim, M., and Ribot-Ciscar, E. (2000). Proprioceptive population coding of two-dimensional limb movements in humans: II. Muscle-spindle feedback during drawing-like movements. Exp. Brain Res. 134, 311–321. doi: 10.1007/s002210000472

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, A. B. (2004). Cortical neural prosthetics. Annu. Rev. Neurosci. 27, 487–507. doi: 10.1146/annurev.neuro.27.070203.144233

PubMed Abstract | CrossRef Full Text | Google Scholar

Scott, S. H., and Loeb, G. E. (1994). The computation of position sense from spindles in mono- and multiarticular muscles. J. Neurosci. 14, 7529–7540.

PubMed Abstract | Google Scholar

Serina, E. R., Tal, R., and Rempel, D. (1999). Wrist and forearm postures and motions during typing. Ergonomics 42, 938–951. doi: 10.1080/001401399185225

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, C. E. (1948). A mathematical theory of communication. Bell Syst. Tech. J. 27, 379–423. doi: 10.1002/j.1538-7305.1948.tb01338.x

CrossRef Full Text | Google Scholar

Simoncelli, E. P., and Olshausen, B. A. (2001). Natural image statistics and neural representation. Annu. Rev. Neurosci. 24, 1193–1216. doi: 10.1146/annurev.neuro.24.1.1193

PubMed Abstract | CrossRef Full Text | Google Scholar

Stein, R. B., Cody, F. W., and Capaday, C. (1988). The trajectory of human wrist movements. J. Neurophysiol. 59, 1814–1830.

PubMed Abstract | Google Scholar

Stein, R. B., Weber, D. J., Aoyagi, Y., Prochazka, A., Wagenaar, J. B., Shoham, S., et al. (2004). Coding of position by simultaneously recorded sensory neurones in the cat dorsal root ganglion. J. Physiol. 560, 883–896. doi: 10.1113/jphysiol.2004.068668

PubMed Abstract | CrossRef Full Text | Google Scholar

Wessberg, J., Stambaugh, C. R., Kralik, J. D., Beck, P. D., Laubach, M., Chapin, J. K., et al. (2000). Real-time prediction of hand trajectory by ensembles of cortical neurons in primates. Nature 408, 361–365. doi: 10.1038/35042582

PubMed Abstract | CrossRef Full Text | Google Scholar