Encoding mechanisms for sensory neurons studied with a multielectrode array in the cat dorsal root ganglion

Recent advances in microelectrode array technology now permit a direct examination of the way populations of sensory neurons encode information about a limb's position in space. To address this issue, we recorded nerve impulses from about 100 single units simultaneously in the L6 and L7 dorsal root ganglia (DRG) of the anesthetized cat. Movement sensors, placed near the hip, knee, ankle, and foot, recorded passive movements of the cat's limb while it was moved pseudo-randomly. The firing rate of the neurons was correlated with the position of the limb in various coordinate systems. The firing rates were less correlated to the position of the foot in Cartesian coordinates (x, y) than in joint angular coordinates (hip, knee, ankle), or in polar coordinates. A model was developed in which position and its derivatives are encoded linearly, followed by a nonlinear spike-generating process. Adding the nonlinear portion significantly increased the correlations in all coordinate systems, and the full models were able to accurately predict the firing rates of various types of sensory neurons. The observed residual variability is captured by a simple stochastic model. Our results suggest that compact encoding models for primary afferents recorded at the DRG are well represented in polar coordinates, as has previously been suggested for the cortical and spinal representation of movement. This study illustrates how sensory receptors encode a sense of limb position, and it provides a general framework for modeling sensory encoding by populations of neurons.

757 ARTICLE / ARTICLE Encoding mechanisms for sensory neurons studied with a multielectrode array in the cat dorsal root ganglion1 R.B. Stein, Y. Aoyagi, D.J. Weber, S. Shoham, and R.A. Normann Abstract: Rcccnt advances in microclcctrodc array technology now permit a direct examination of the way populations of sensory neurons encode information about a limb's position in space. To address this issue, we recorded nerve im­pulses from about 100 single units simultaneously in the L6 and L7 dorsal root ganglia (DRG) of the anesthetized cat. Movement sensors, placed near the hip, knee, ankle, and foot, recorded passive movements of the cat's limb while it was moved pseudo-randomly. The firing rate of the neurons was correlated with the position of the limb in various co­ordinate systems. The firing rates were less correlated to the position of the foot in Cartesian coordinates (.v, y) than in joint angular coordinates (hip, knee, ankle), or in polar coordinates. A model was developed in which position and its derivatives arc encoded linearly, followed by a nonlinear spike-generating process. Adding the nonlinear portion signifi­cantly increased the correlations in all coordinate systems, and the full models were able to accurately predict the fir­ing rates of various types of sensory neurons. The observed residual variability is captured by a simple stochastic model. Our results suggest that compact encoding models for primary afferents recorded at the DRG arc well repre­sented in polar coordinates, as has previously been suggested for the cortical and spinal representation of movement. This study illustrates how sensory receptors encode a sense of limb position, and it provides a general framework for modeling sensory encoding by populations of neurons. Key words: sensory, encoding, multielectrode, dorsal root ganglion, cutaneous, muscle. Resume : Grace aux progrcs recents dc la technologic des reseaux dc microelectrodes, nous pouvons maintenant obser­ver directement comment les populations dc neurones sensoriels encodent l'information relative a la position d'un mcmbre dans l'espace. Pour etudier cette question, nous avons enregistre les impulsions nerveuses d'une centaine d'unites discriminables, simultanement dans les ganglions dc la racine dorsale (GRD) L6 et L7 d'un chat sous anes- thesie. Des detecteurs du mouvement, places pres dc la hanche, du genou, dc la cheville et du pied, ont enregistre les mouvements passifs du mcmbre du chat pendant que celui-ci etait deplace dc fa?on pseudo-aleatoire. La cadence des neurones a etc mise en correlation avec la position du mcmbre dans divers systcmes dc coordonnees. La correlation entre la cadence et la position du pied etait moins grande dans les coordonnees cartesienncs (.v, y) que dans les coor- donnees articulaires (hanche, genou, cheville) ou dans les coordonnees angulaires. Un modclc a etc con^u, dans lequel la position et ses derives sont encodes dc fa^on lineaire, suivi par un processus non lineaire generateur dc pointes. L'ajout dc la portion non lineaire augmente dc fa?on significative les correlations dans tous les systcmes dc coordon- nees et les modclcs complcts ont pu predire avec precision la cadence des divers types dc neurones sensoriels. La va- riabilite residuclle observee est bien saisie par un simple modclc stochastique. En conclusion, nos resultats semblent indiquer que les modclcs d'encodage compact pour les afferents primaires cnregistres au GRD sont bien reprcscntes dans les coordonnees angulaires, commc ccla a deja etc indique pour la representation corticale et spinale du mouve­ment. Cette etude illustre comment les recepteurs sensoriels encodent la position d'un mcmbre et fournissent un cadre general pour la modelisation dc l'encodage sensoriel par des populations dc neurones. Received 6 July 2004. Accepted 6 July 2004. Published on the NRC Research Press Web site at http://cjpp.nrc.ca on 1 October 2004. R.B. Stein,2 Y. Aoyagi, and D.J. Weber. Centre for Neuroscience, University of Alberta, 513 Heritage Medical Research Institute, Edmonton, AB T6G 2S2, Canada. S. Shoham and R.A. Normann. Department of Bioengineering, University of Utah, Salt Lake City, UT, USA. JThis paper is one of a selection of papers published in this Special Issue, entitled Nerve, muscle, and beyond: the R.B. Stein Symposium. Corresponding author: (email: [e-mail]). Can. J . Physiol. Pharmacol. 8 2 : 7 5 7^ 7 6 8 (2 0 0 4 ) doi: 1 0 .1139/Y0 4 -0 7 5 © 2 0 0 4 NRC Canada Mots cles : scnsoricl, codagc. multi-clcctrodcs, ganglion dc la rarinc dorsalc, cutanc, musdc. [Traduit par la Redaction] 758 Can. J. Physiol. Pharmacol. Vol. 82, 2004 Introduction The senior author of this article (R.B.S.) has had a life-long interest in the mechanisms involved in encoding information in sensory neurons (Matthews and Stein 1969; Stein 1967). Encoding models describe the transformation that takes place from a sensory stimulus to a neural response, as well as the variability of this response (Stein 1965). Early work was based on a simple linear summation of length and velocity in­puts to muscle receptors, and the variability of the discharge was shown to limit the amount of information that could be transmitted about these variables. Much more complex non­linear models have since been developed for muscle spindles and other sensory receptors (Hasan and Houk 1975; Otten et al. 1995; Prochazka and Gorassini 1998i>). Similarly, it has been suggested that the variability in sensory discharge con­veys information, rather than simply being neural "noise" (Perkel and Bullock 1968; Rieke et al. 1997; Stein 1970). The nervous system is a massively parallel computing sys­tem that deals with a vast amount of activity in real time. Can we begin to understand the general processing of activity in sensory systems, over and above our detailed understand­ing of single receptors? With the development of multi­electrode recording techniques (Loeb et al. 1977; Prochazka et al. 1976; Serruya et al. 2002; Taylor et al. 2002; Wessberg et al. 2000), we can now record from many neurons and study properties of the ensemble of signals passing through a structure, such as the dorsal root ganglion. Previously, sin­gle neurons were recorded sequentially over many hours or many preparations to build up a population, with the hope that conditions remained constant over time and between an­imals. Now, by recording from a population of neurons si­multaneously, that assumption is no longer required; the nature of the processing can be addressed in single experiments. To deal with a large number of simultaneous signals, we have developed a variety of automated methods. This ongo­ing development will be described in this paper. We also found that the activity of the sensory neurons can be de­scribed in a general framework. The responses still contain length and velocity components, but they can be described in terms of extrinsic variables, such as Cartesian (x, y) or polar coordinates, and intrinsic variables, such as joint angles. In addition to linear correlations, we found systematic non­linearities that are well described by a Wiener cascade. A Wiener cascade has a linear filter followed by a static non­linearity (Kearney and Hunter 1990; Korenberg and Hunter 1999). This has an intuitive appeal, because the ionic cur­rents generated by sensory inputs or synaptic potentials in a neuron may sum more or less linearly, but there are inherent nonlinearities in spike generation (see Chander and Chichilnisky 2001; Shoham 2001). For example, the firing rate saturates at a maximum determined by the neuron's re­fractory period and at a minimum of zero (no firing). Even with the best linear/nonlinear models, residual variability still remains between the actual firing and predicted patterns. This variability can be treated as a stochastic process, the properties of which are well described by our model. The major goal of this paper is to describe the methods and structure that we have developed to understand the encoding mechanisms in sensory neurons. Methods A series of cats were studied, during which methods grad­ually improved. Three acute experiments were done at the University of Utah (Salt Lake City, Utah). A manipulandum with 2 degrees of freedom was attached to the cat's paw, and the foot was moved manually in the sagittal plane. From the 2-joint angles of the manipulandum, the position of the cat's paw was calculated and compared with the neural signals re­corded. Three acute experiments were also done at the Uni­versity of Alberta (Edmonton, Alta.). The local animal welfare committees at the respective universities approved all procedures, which were conducted in accordance with the Helsinki Convention. The results from all experiments were essentially identical, and we have selected examples to illus­trate the common findings. Movement tracking In the first 2 experiments at the University of Alberta, movements of the limb were produced manually (Fig. 1), and were recorded by electromagnetic motion-tracking sen­sors (6D-Research, Skill Technologies Inc., Phoenix Ariz.). Four sensors were used: sensor 1 was placed on the skin near the hip joint, sensor 2 was placed on the lateral epicondyle of the femur near the knee, sensor 3 was placed on the lateral malleolus of the tibia near the ankle, and sen­sor 4 was placed on the surface of a foot holder near the lat­eral metatarsal joint of the foot. To avoid skin slippage or displacement during movement, sensors 2 and 3 were rigidly fixed to the femur and tibia with surgical sutures through holes drilled in the respective bones. No instruments near the sensors, including sections of the spinal frame, contained metal, so that there was no electromagnetic interference with the signals recorded from the motion sensors. The distance of each sensor from neighboring joints was measured. In the most recent experiment, a robotic manipulandum with 2 DC servo motors (Parker BE233DJ) was programmed to deliver repeatable movements. To generate random move­ments, the manipulandum moved sequentially from one po­sition to the next. The positions were selected at random from a rectangular grid of points, and the velocity of each movement was selected at random over a range of speeds up to 0.6 m/s. The movements continued until all points on the grid had been reached, so there was a uniform coverage of the workspace. Movements were tracked in this experiment with a digital video camera (JVC DVL9800) at 120 frames/s, using 5-mm reflective markers placed over the hip, knee, ankle, foot, and toes. No invasive procedures were used to track motion, so that the results could be compared with results from awake behaving animals (Prochazka et al. 2003; Weber et al. 2002). Because the skin over the knee moved to some extent, the knee position was calculated using the hip and ankle mark- © 2 0 0 4 N RC Canada Stein et al 759 Fig. 1. Movement sensors (O) were placed at the hip, just above the knee, at the ankle, and on the paw near the metatarsal- phalangeal joint. From the positions of the sensors, a stick figure of the cat's hindlimb in the sagittal plane was calculated. The random movement of the paw over its passive range of motion during a 45-s trial is shown as a dotted line. The position of the paw can be represented in terms of the forward (.v) and vertical (v) position with respect to the hip (Cartesian coordinates). It can also be represented in polar coordinates as the distance (r) and the angle (p) of the paw with respect to the hip. Finally, the position can be represented in terms of the forward flexion of the hip (h) (with respect to the vertical), the flexion of the knee (k) (with respect to a straight knee), and the flexion of the ankle (a) with respect to a right angle (dorsiflexion). For simplicity, the ankle angle is not shown. The angles a and p have negative values in this figure. Forward (x) position (cm) ers and the known lengths of the femoral and tibial seg­ments. In some cases, the data were smoothed, using the decimate function in Matlab, to reduce the effective sample rate and, hence, the number of data points. Surgical procedure The animals were anesthetized with sodium pentobarbital (University of Alberta) or Halothane (University of Utah). A tracheal cannula was inserted to maintain respiration, and a jugular catheter was used to add fluids and anesthetic. The back was shaved and a skin incision was made along the midline of the back. Paraspinal muscles overlying the trans­verse processes of L5-S1 were removed, and a wide lamin­ectomy was performed to expose the spinal cord and dorsal roots. Two 5 x 10 arrays of penetrating microelectrodes (Utah Electrode Arrays, Bionic Technologies, Salt Lake City, Utah) were implanted in the L6 and L7 dorsal root ganglia (DRGs) on one side through the dura with a pneu­matic inserter (Rousche and Normann 1992). Reference wires were placed in the fluids surrounding the DRGs. The animals were suspended from a spinal frame and the skin was closed. Radiant heat was used to maintain the body tem­perature near 37 °C. Unit identification To identify individual muscle and joint receptors, the hip, knee, ankle, and toes were moved manually. Vibration over the tendon or belly was used to identify primary muscle spindle afferents, where possible. Palpation (touch, pressure, and pinch) and vibration were used to identify cutaneous re­ceptors, and gentle blowing or touch were used to identify hair receptors. The activity was recorded for 10 s from each electrode channel during the "best" manipulation for a given unit. Once all the units had been identified, the series of movements described above was applied. Data analysis The electrode arrays were connected to a 100-channel am­plifier. The gain of the amplifiers was 5000 (bandwidth from 250-7500 Hz), and each electrode was sampled at 30 kHz. A Pentium class personal computer saved the responses with a data acquisition system (Neural Signal Acquisition System (NSAS), Bionic Technologies). In this system, after thresh­olds are set, only short waveforms are saved around the time that the threshold is crossed (Guillory and Normann 1999). Offline, single units were discriminated from the set of re­corded waveforms on each electrode, using a Matlab-based algorithm developed by Shoham et al. (2003). In the Alberta experiments, the NSAS recorded the units, as described above, and the movements were recorded separately. A marker pulse was recorded on both systems for synchroniza­tion. Statistical tests were done using Student's t tests. Cor­relation coefficients were calculated using the "corrcoef' Matlab function. Results Figure 2 shows the responses of 64 neurons, the discharge of which varied during random movements of the limb ap­proximately in the sagittal plane. These were selected from about 100 neurons that were recorded simultaneously. Elec­trodes 1-50 recorded in the L7 DRG and electrodes 51-100 recorded in the L6 DRG. Figure 2 also shows the position of the cat's paw at the level of the metatarsal-phalangeal joints with respect to the hip position, which was fixed. Forward positions (x-direction) are positive and backward positions negative. Similarly, upward positions (v-direction) are posi­tive and downward positions negative. Because the paw was always below the hip, the y-values are all negative. The medio-lateral movements of the limb were also recorded in this experiment, but the deviations from the sagittal plane were small (<2 cm in either direction). With the large number of channels being recorded, we tried to automate the discrimination and preprocessing of the neural data as much as possible. Several spike waveforms with distinct firing patterns were automatically recognized on a single channel (see Methods and Shoham et al. 2003). However, the superposition of waveforms and the presence of noise in the recordings can lead to spikes being missed or misclassified. Figure 3 shows the effect on the instantaneous frequency. If a spike is missed, the instantaneous frequency falls to about half that of its neighbors (Fig. 3a). A log2 scale is used for instantaneous frequency because the fre­quency will drop about 1 log2 unit on such a scale and fall outside the standard error bars of neighboring intervals. In­serting an extra spike (Fig. 3b) will produce 2 intervals that lie within the error bounds. Similarly, an extra misclassified spike will create an instantaneous frequency at least 1 log2 unit above its neighbors. Eliminating the spike brings the in- © 2 0 0 4 N RC Canada 760 Can. J. Physiol. Pharmacol. Vol. 82, 2004 Fig. 2. The neural responses of 64 neurons varied with random movements of the cat's hindlimb. Each tick mark represents a single spike for 10 s out of a 45-s recording. The forward and vertical positions of the paw arc shown for comparison. The electrode num­bers in terms of the Utah Electrode Arrays arc given on the right side. Values after the decimal point indicate different units recorded from some electrodes. 60 50 40 CO H-» ‘c 30 iiiiiiiiiiiiiiiiiiiiiiili.!....In I in I in in ini 89.0 79.3 !!I!I]]IJ!..!!!!J J!!!.!! 79.1 70.4 70 2 .................. 69.0 INH I11 JillilI! I 1 I1" III1 67.0 63.0 llllllllllll 59.1 "I I I I I I 55.3 55.1 .........-~.... .....s „„„„„ ............ ml iniill" 1.....i......1............i......... .......... 95.0 86.1 79.2 75.0 70.3 70.1 68.1 65.0 60.3 59.3 58.0 55.2 54.0 45.0 43.2 37.2 36.2 35.0 31.3 31.1 28.4 28.2 27.2 25.2 24.1 23.0 21.1 15.2 14.2 13.3 13.1 11.2 30 32 34 36 Time (s) 38 40 tervals within bounds. This method of spike editing has been tested on simulated and actual neural data and has been found to significantly improve the accuracy of the records (Stein and Weber 2004). Next, the spikes had to be aligned with respect to the times at which the kinematic data were sampled. Typically, the number of spikes is counted in a given sample interval [s,- si+ J, where At = si+1 - Sj is the sample period for the ki­nematic data. When this is done, the timing within the sam­ple period is lost and the times are delayed on average by half a sample period. This can be avoided by "partially binning" the data (Richmond et al. 1987; Schwartz 1992) and making the weight at sample times s, and s, + j propor­tional to the magnitude of the time difference tj - s, , where tj is the time of the jth spike. The incremental nature of the firing rate, rh is described by © 2 0 0 4 N RC Canada Stein et al. 761 Fig. 3. (A) A missed spike (O) produces a lower frequency in­terval than the standard deviation (SD) of the four neighboring intervals (long bars represent two SD units from the mean). If a spike is missing, the frequency falls below the short horizontal bar. (B) Inserting a spike in the middle of the interval produces two intervals that arc within the error bounds. (C) An extra spike inserted in the train can produce 1 or 2 instantaneous frequencies (O) that arc above the expected SD (short error bar). (D) Re­moving 1 spike reduces the error bounds (long horizontal bars) and leaves one interval, which lies within the new error bounds. From Stein and Weber (2004). (A) Missed spike 6 o> 4 o ^ (B) After inserting spike © O) o 5 10 (C) Extra spike © © * . ©© 0 5 10 (D) After deleting spike © [1] if t; 0 5 10 Interval number Arf 0 5 10 Interval number At < At and the summation is over all spike occur­rences (j). In other words, if the spike occurs at time sh it will contribute a total weight of 1/Ar to the firing rate. The partial weights at the adjacent sample times drop off linearly until reaching zero if ts - s,- = At. One can easily show that the average time of the contribution is t-r the time when the spike actually occurred, so no time delays are introduced. The factor 1/Ar converts the weights to a firing rate in im­pulses/ s. Figure 4 illustrates the process. Finally, the rate function was filtered with a critically damped second-order low-pass filter. This corresponds to re­placing the rate occurring at s, with an excitatory postsynaptic potential (EPSP)-like waveform. The resulting function fi t) is described by [2] fit) = E rp.Ht - tj)exp(.-akt - r,)) The rate constant (a) used was 15 or 30 rad/s (approx. 2.4 or 4.8 Hz), corresponding to a time constant of 67 or 33 ms, depending on the frequency components in the input. The resultant waveform is analogous to the summed EPSPs oc­curring in a postsynaptic neuron. Before examining the relationship between the firing rate and kinematic variables, we should discuss the statistical properties of the kinematics. The same filtering procedures were applied to the kinematics, so that no time shifts were produced between filing rates and kinematics. Table 1 shows Fig. 4. Methods of sampling and filtering spike trains. A scries of spikes in a sensory neuron (top) must be sampled at discrete times (O). Each spike is convolved with a triangular kernel (equation 1), which assures that each spike is sampled with the same weight and without any net time shift. The rate is then smoothed with an excitatory postsynaptic potential (EPSP)-likc filter (equation 2), which introduces delays. The kinematics arc also filtered in the same way to prevent any net time shift. 1 I 0.5 % (A) -e-6-e- -e-e-e-e- -e-6-e- -e-6-e-s 0.1 0.2 0.3 0.4 0.5 -* 0.6 Table 1. Correlation coefficients between kinematic variables. (A) Cartesian coordinate system: forward (.v) vs. vertical (v) position and their derivatives (d). dx dx x dx x -0.03 -0.29 -0.21 0.11 -0.09 -0.10 (B) Polar coordinate system: distance from hip (r) vs. angle from vertical (p) and their derivatives (d). dr dp r dr P -0.10 0.53 -0.11 -0.05 0.59 -0.14 (C) Angular coordinate system: ankle (a), knee (k), and hip (/1) angles and their derivatives (d). da k dk h dh a -0.08 0.4 0.01 0.79 -0.23 da -0.14 0.75 0.12 0.64 k -0.08 0.1 -0.13 dk 0.02 0.61 h -0.08 the correlation coefficients between the x and y positions (Cartesian coordinates) and their derivatives, dx and dy. The movements were quite random, even though they were ap­plied manually in this experiment. The small correlation be­tween x and y arises, at least in part, because the passive range of motion is not rectangular. Backward movements tend to be associated to some extent with upward move- © 2 0 0 4 N RC Canada 7 62 Can. J. Physiol. Pharmacol. Vol. 82, 2004 ments. Table 1 also shows the correlation for other coordi­nate systems. The position can be described in polar coordinates as the distance (/•) of the paw from the hip, and the angle (/?) with respect to the vertical. The correlation be­tween r and p is greater than x and v, because the anatomy of the limb allows the limb to be fully extended (maximum r) only when the limb is forward. Finally, the position of the limb can be described in terms of the angles of the ankle (a), the knee (k), and hip (h), and the length of the limb seg­ments. In contrast to the relatively small correlation for Car­tesian coordinates, the correlation reaches a value of 0.79 between the ankle and hip. This is because of the bio­mechanical linkages between the limb segments and the pas­sive stiffness of the joints. Bosco et al. (1996) noted this after plotting the joint angles in a 3-dimensional figure; they found that the points formed a plane. In other words, al­though there are 3 separate joints, the number of independ­ent combinations during passive movements is effectively reduced to two dimensions. The same procedures were applied to all neurons. The va­riety of neurons recorded has been described in detail else­where (Aoyagi et al. 2003), so this topic will not be discussed further. The vertical axis in the "stem" plots of Fig. 5 is the firing rate of a neuron (unit 14.2 in Fig. 1) mea­sured in 0.1-s intervals. The other axes are the kinematic variables measured, for example the x and y positions of the paw with respect to the hip in Fig. 5a. This neuron re­sponded best to flexion of the ankle and toe extension, and could be a muscle spindle afferent in the flexor hallucis longus or flexor digitorum longus muscles. Figures 5b and 5c show the corresponding plots, using polar coordinates (distance and orientation with respect to the hip) and angular coordinates (ankle and knee joint angles), respectively. Be­cause of the high correlation between the hip and ankle, the hip angle is not shown. Two simple points emerge from Fig. 5. First, the firing of the neuron can be described in any of the 3 coordinate sys­tems. In Cartesian coordinates, the response is highest when the limb is forward and up. In polar coordinates, it is highest when the limb is oriented forward and the distance is small (i.e., when the limb is raised). In angular coordinates, it re­sponds best to flexion of the ankle, as expected from its identification as a spindle in an ankle extensor. Second, the response is not fully determined in any of the coordinate systems; the firing rates at adjacent positions can vary sub­stantially. The correlation coefficient between firing rate and a linear combination of x and y values plus a constant was 0.64. The correlation in polar coordinates was 0.78 and in angular coordinates was 0.77. All angles (ankle, knee, and hip) were used in calculating the correlation. The correlation coefficient for Cartesian coordinates is significantly less than the others; there were 445 data points in this sample (for the sake of clarity, not all of which are shown in Fig. 5). As described in the Introduction, our goal was to see whether a relatively simple model can describe the coding properties of the wide variety of large sensory neurons re­corded here. For example, the fluctuations in rate at a given position could arise if the neurons respond, not only to posi­tion, but also to derivatives of position, such as velocity and acceleration. If these variables are important, adding them to Fig. 5. Three-dimensional "stem" plots showing the firing rate (vertical axis), over periods of 0.1 s, as a function of position for the three coordinate systems (see Fig. 1). The unit shown here (unit 14.2) responded best to flexion of the ankle and extension of the toes, and may have been a muscle spindle afferent in the flexor hallucis longus or flexor digitorum longus muscles. (A) unit 14.2 Knee angle (k) Ankle angle (a) the linear model should increase the correlation. In fact, add­ing velocity terms increased the correlation coefficients to 0.67 (Cartesian coordinates), 0.79 (polar coordinates), and 0.78 (angular coordinates) for the cell in Fig. 5. Figure 6 shows the results after analyzing all 64 cells in this experiment. Similar results were seen in other experi­ments. Correlating the firing rate to the current and previous positions (shifted by 0.1 s) is equivalent to correlating the response to position and velocity, because the velocity v = [*(?) - x(t - 8t)]/dt, where 8? = 0.1 s. Adding a third time point for position will include acceleration, and so on, for higher derivatives. The ordinate in Fig. 6 is the number (??) of time points used in the linear fit (n = 1 uses the positions alone, n = 2 includes the velocities in each coordinate, n = 3 includes velocities and accelerations, etc.). Adding velocity terms improves the average correlations significantly in all coordinate systems. The additional increase is not significant for higher derivatives, despite the fact that adding more vari­ables can only increase the value of the correlation. Shoham (2001) found similar results when modeling the encoding properties of motor cortical neurons during voluntary arm movements. © 2 0 0 4 N RC Canada Stein et al. 763 Fig. 6. Mean correlation coefficients for all 64 cells, relating the number of spikes to kinematic variables in the 3 coordinate sys­tems. The coefficients were significantly lower in the Cartesian coordinates than in the angular coordinates (paired t test), but the values were not significantly different in the polar and angular coordinate systems. The correlation was calculated with the cur­rent time points (1), or after including previous time points (2­6), each separated by 0.1 s. This is equivalent to including the effects of velocity (2), acceleration (3), and higher differentials (4-6). Including velocity increased the correlation significantly in all coordinate frames, but increasing the number of time points did not significantly increase the correlation. 0.8r ■ Cartesian ■ Polar HI Angles I 0.6 (D O o c o CD L_ L_ o o 0.4 0.2 0 3 4 Order For all orders, the correlation was significantly less for Cartesian coordinates than for angular coordinates, but the difference was generally not significant between the polar and angular coordinates. Note that the angular coordinate system included 3 angles and their derivatives, whereas the polar coordinate system included only 2 coordinates (r and p) and their derivatives. However, as pointed out above, the 3 joint angles were not independent. The Discussion will consider the implications of these results of representing sensory information in different coordinate frames. Figure 7 shows the actual number of spikes, in intervals of 0.1 s, plotted against the predicted number, calculated from the lineai' model with position, velocity, and acceleration terms in polar coordinates for three cells (units 14.2, 22, and 75 in Fig. 2). Unit 22 (top) fired at the lowest rate. The main nonlinearity is the saturation near zero, because less-than- zero spikes are not possible. Unit 75 fired at the fastest rate and showed a curvature at the upper end, suggesting that it was approaching its maximum firing rate. Unit 14.2, the unit shown in previous figures, demonstrated intermediate behav­ior. A cubic polynomial fitted the curvature for these and other cells, as was also seen in the motor cortex (Shoham 2001) and in a simple neural model (described in Appendix A). Including the nonlinear term increased the correlation coefficients of unit 14.2 from 0.78 to 0.84, unit 22 from 0.64 to 0.71, and unit 75 from 0.88 to 0.91. Figure 8 shows the results for all the cells in all coordi­nate systems. The increase is significant in all coordinate systems. Again, the values are lower for the Cartesian than Fig. 7. The actual number of spikes in 0.1-s intervals increases nonlinearly with the predicted number of spikes (left) for three typical units. The data can be well fitted with a cubic (third or­der) polynomial (solid lines). The variance (right) increases with the actual number of spikes initially, but may reach a peak and decline as the number increases. Data arc fitted with a parabolic (second order) polynomial. Imp., impulse. the other coordinate systems, and there is no significant difference between the polar and angular representations. Although the fit is improved, a percentage of the discharge is still not accounted for and can be assumed to be "neural noise." Figure 8 also shows the variance of the deviation from the fitted Wiener cascade. The variance initially in­creases approximately linearly with the firing rate, but can reach a peak and even decline at high levels. The simple neural model described in Appendix A can also simulate this behavior. Simulated data will be compared with the biologi­cal results in the Discussion. The neural noise might be associated with the stochastic generation of nerve impulses, and filtering the time histo­grams can reduce this noise. If so, the filtering for the EPSP waveform would reduce the noise and improve the correla­tion. Therefore, only values in the white bars in Fig. 8 were filtered. The filtering did improve the correlations, and the value for the filtered rates in angular coordinates almost reached 0.7 on average over the population of 64 cells, com­pared to an initial value of less than 0.5 (Fig. 6). How well does this model explain the behavior of the sen­sory cells'? Figure 9 shows the time course of the firing rate over the whole random trial for the same three units that were considered in Fig. 7. Dotted lines show the predictions of the nonlinear model after filtering for each unit. Although there are still deviations between the actual and predicted firing rates, the major features of the bursting patterns are followed quite faithfully. Discussion We recorded about 100 sensory neurons with micro­electrode arrays, while the limb was moved randomly through most of its passive range of motion in the sagittal © 2 0 0 4 N RC Canada 764 Can. J. Physiol. Pharmacol. Vol. 82, 2004 Fig. 8. The correlation coefficient in all coordinate systems can he increased using the nonlinear fit (gray bars), shown in Fig. 7 (left), rather than the linear fit (black bars). Filtering the output can further increase the correlation (white bars) by reducing the neural noise, which is represented by the variance in Fig. 7 (right). 0.8r plane (Fig. 1). From these recordings, we established how populations of neurons code global variables, such as foot position in single animals, without introducing the inter­animal variability inherent when collecting neural popula­tions from many animals. We also tried to develop a general model for the variety of neurons in the DRG. Essentially, the model contains 3 parts: a linear part related to the posi­tion of the limb in space, and its derivatives; a nonlinear part related to impulse generation; and a residual component of neural noise. Each of these will be discussed separately in relation to outstanding issues of neural encoding. The linear portion and coordinate systems for sensory coding Muscle receptors have long been known to respond lin­early over a certain range of lengths (Matthews 1972; Terzuolo and Washizu 1962). This study is the first to test many sensory receptors in terms of different coordinate sys­tems. We assumed a priori that the system of joint angles would be the most natural, particularly for muscle receptors that span one or more joints. Cartesian and polar coordinates can only be obtained with calculations using nonlinear equa­tions from joint angles and segment lengths. The relation be­tween muscle length and joint angle can be quite nonlinear in certain ranges (Goslow et al. 1973). However, if the re­ceptors respond reasonably linearly to joint angles and their derivatives over a restricted range, then the correlations should be higher for this system than for Cartesian or polar coordinate systems. The system of joint angles did give sig­nificantly higher correlations than the Cartesian coordinates, but was generally not significantly different from the polar coordinates. The lack of a significant difference is interest­ing for several reasons: First, polar coordinates define the foot position in an ex­trinsic reference frame (with respect to the hip), whereas joint angles are measured in intrinsic coordinates that are Fig. 9. Time course of the actual (solid line) and predicted (dot­ted line) firing rates for the three neurons shown in Fig. 7. more closely related to receptor variables, such as muscle stretch (Bosco et al. 2000), as mentioned above. Second, there are 3 joint angles and only 2 polar coordi­nates; an extra fitted parameter might be expected to im­prove the fit and, strictly speaking, should have been penalized in the comparison. However, the correlation be­tween angles means that there are effectively only 2 inde­pendent angular variables. Third, research over the past 20 years has suggested that the motor cortex plans movements in terms of distances and directions of an endpoint (hand or foot), essentially a polar coordinate system (Chapin et al. 1999; Georgopoulos et al. 1982, 1992; Nicolelis 2001; Schwartz 1994;"). Others have argued that many cortical neurons code for muscle activity independent of the direction of the endpoint movement (Mussa-lvaldi 1988; Scott and Kalaska 1997; Kakei et al. 1999; Scott and Norman 2003). Without going into the de­tails of this controversy, the question then arises: Can joint- based receptors provide appropriate sensory information to guide movements of an endpoint (Todorov 2000)'? Recording from cells in the dorsal spinocerebellar tract (DSCT), Bosco and Poppele ( 2001) found that the discharge of most cells at this level was already well described in terms of movements of the endpoint, rather than of single joints. They suggested that convergence from a number of sensory receptors could provide this global information. Our results show that, even at the level of the sensory re­ceptors, responses can be well described in terms of whole- © 2 0 0 4 N RC Canada Stein et al. 765 limb movements. Individual cutaneous receptors can re­spond exquisitely to pressure to or deformation of a particu­lar part of the skin, but they also respond when the whole limb is moved, which stretches the skin of their receptive fields. Similarly, muscle receptors respond to movements of particular joints, but joint movements are linked by the bio­mechanics of muscles, ligaments, etc. Thus, at least in the passive state, muscle receptors are well correlated to move­ments of the whole limb (Loeb et al. 1985; Prochazka and Gorassini 1998a). The situation is more complex during ac­tive movements; the contraction of a and y motoneurons modify the discharge of muscle receptors, such as Golgi ten­don organs and muscle spindles (Prochazka 1989). Our ex­periments are being extended to behaving animals to see how active movements modify the results (Prochazka et al. 2003; Weber et al. 2002). Important extensions of these ex­periments include ones that will explore 3-dimensional movements, as well as ones where joint angles will be al­lowed to vary while the endpoint is fixed. The latter experi­ment may allow a direct dissociation of the endpoint/joint- angle ambiguity, but is likely to introduce new ambiguities because of the pressures applied by a highly unnatural con­straint. Finally, cutaneous, muscle, and joint receptors all contrib­ute to a proprioceptive sense of the limb in space (Gandevia et al. 1983). These experiments show how a generalized sense of limb position may arise already from the discharge of sensory receptors. This sensory information is processed through several systems. For example, the dorsal columns may preserve the specificity of the individual neural re­sponses for sensory perception of the nature of the stimulus and the area of skin that is being stimulated, whereas the DSCT combines responses in an appropriate way for the control of limb movements. This concept of multiple repre­sentations is analogous to that found in the visual system, where color, form, and distance (stereopsis) are processed separately (Zeki 1993), and is consistent with the multiple representations of the limb in different sensory and motor cortical areas (Rothwell 1994). The nonlinear portion and impulse generation The generation of nerve impulses inherently involves nonlinearities, such as the saturation nonlinearities for strong inputs (near-maximal firing rates determined by the refrac­tory period) and subthreshold inputs (rates near zero). A simple model is described in Appendix A, in which sub­threshold inputs sum linearly and decay exponentially with time. When a threshold is reached, the membrane potential is reset to the resting value and the process begins again. Longer time-constant processes, such as after hyperpolari­zations, which sum following successive action potentials, or slow processes of sensory adaptation, can be added, but the essential character of the nonlinear process does not change. In our experiments (Fig. 7) and in the simulations (Fig. 10), a cubic polynomial provided a good representation. In the motor cortex, Shoham (2001) tried polynomials up to ninth order, but also settled on a cubic polynomial. Neural noise and sensory discrimination The variability in neural responses will ultimately limit the ability to discriminate among different stimulus intensi­Fig. 10. Number of spikes produced by the neural model de­scribed in Appendix A for the random input. The actual number of spikes in 0.1 s and its variance arc shown on the right, and can be compared with the experimental data shown in Fig. 7. The fitted curves arc cubic and parabolic polynomials, as in Fig. 7. Neural model o o ties. The variance of the neural response initially increased linearly with the mean rate (Fig. 7). This has also been found to be the case in the motor system (Lee et al. 1998; Maynard et al. 1999; Shadlen and Newsome 1998). For some neurons in our sample, the variance reached a peak and declined (Fig. 6). These were typically neurons that showed saturation at high firing rates. As a result, large changes in input don't vary the firing rate substantially. The same result can be simulated in the neural model, with the variance initially increasing linearly with mean rate, but then reaching a peak and declining (see Fig. 10 and Appendix A). The idea of just-noticeable differences in sensory systems has a long history, going back to Weber and Fechner in the middle of the 19th century (Kandel et al. 1991). The Weber- Fechner law postulated that the ratio D/S was a constant, where D is the just-noticeable difference and S is the stimu­lus intensity. Experimentally, the ratio D/S is typically a de­clining function of stimulus intensity (Stevens 1961). In our cells, the variance was initially proportional to the mean rate, so the ratio of the standard deviation to the mean de­clined as the square root of the mean intensity. Studies over the past 40 years have suggested that, generally, the just- noticeable difference is a power function of the stimulus in­tensity (Stevens 1961). For our neurons, the power function would be 0.5 until the response begins to saturate. The expo­nent in the power function is also affected by nonlinearities in the sensory transduction processes, which have not been considered here. Comparisons of neural responses to the psychophysical estimates of stimulus intensity have been largely confined to the use of punctate stimuli that excite one or a few sensory receptors (e.g., Werner and Mount- castle 1965). Being able to record from a large population of sensory neurons should permit a direct comparison in mo­ © 2 0 0 4 N RC Canada 766 Can. J. Physiol. Pharmacol. Vol. 82, 2004 dalities, such as joint angle or pressure on an area of skin, that excite a large number of receptors. In conclusion, the ability to record simultaneously from many sensory neurons should permit many basic neurosci­ence questions to be addressed that have to date been difficult to study. We are currently extending these record­ings from the acute state to a chronic behaving animal (Prochazka et al. 2003; Weber et al. 2002), which should al­low study of other basic questions. In addition, the applica­tion of these methods to control stimulation of human muscles in paralyzed limbs after a spinal cord injury or stroke may be possible. This remains an exciting future ave­nue for exploration. Acknowledgements This work was supported by grants to R.B.S. from the Ca­nadian Institutes of Health Research, and to R.A.N. from NIH. Current addresses; Dr. Y. Aoyagi, Department of Reha­bilitation Medicine, Kawasaki Medical School, Okayama, Japan; Dr. D.J. Weber, Faculty of Physical Education, Uni­versity of Alberta, Edmonton; Dr. S. Shoham, Department of Molecular Biology, Princeton University, Princeton N.J. References Aoyagi. Y.. Stein. R.B.. Branner. A.. Pearson. K.G.. and Normann. R.A. 2003. Capabilities of a penetrating mieroeleetrode array for recording single units in dorsal root ganglia of the eat. J. Neurosci. Methods 128: 9-20. Bosco. G.. and Poppele. R.E. 2001. Proprioception from a spinocerebellar perspective. Physiol. Rev. 81: 539-568. Bosco. G.. Rankin. A., and Poppele. R.E. 1996. Representation of passive hindlimb postures in cat spinocerebellar activity. J. Neurophysiol. 83: 715-726. Bosco. G.. Poppele. R.E.. and Eian. J. 2000. Reference frames for spinal proprioception: limb endpoint based or joint-level based? J. Neurophysiol. 83: 2931-2945. Chander. D.. and Chichilnisky. E.J. 2001. Adaptation to temporal contrast in primate and salamandar retina. J. Neurosci. 15: 9 904-9 916. Chapin. J.K.. Moxon. K.A.. Markowitz. R.S.. and Nicolelis. M.A.L. 1999. Real-time control of a robot arm using simulta­neously recorded neurons in the motor cortex. Nat. Neurosci. 2: 664-670. Gandevia. S.C.. Hall. L.A.. McCloskey. D.I.. and Potter. E.K. 1983. Proprioceptive sensation at the terminal joint of the mid­dle finger. J. Physiol. 335: 507-517. 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. Georgopoulos. A.P.. Ashe. J.. Smyrnis. N.. and Taira. M. 1992. The motor cortex and the coding of force. Science (Washington. DC). 256: 1692-1695. Goslow. G.E.. Reinking. R.M.. and Stuart. D.G. 1973. The cat step cycle: hind limb joint angles and muscle lengths during unre­strained locomotion. J. Morphol. 141: 1-42. Guillory. K.S.. and Normann. R.A. 1999. A 100-channel system for real time detection and storage of extracellular waveforms. J. Neurosci. Methods. 91: 21-29. Hasan. Z.. and Houk. J.C. 1975. Analysis of response properties of deefferented mammalian spindle receptors based on frequency response. J. Neurophysiol. 38: 663-672. Kakei. S.. Hoffman. D.S.. and Strick. PL. 1999. Muscle and move­ment representations in the primary motor cortex. Science (Washington. D.C.). 285: 2 136-2 139. Kandel. E.R.. Schwartz. J.H.. and Jessell. T.M. 1991. Principles of neural science. Elsevier. Amsterdam, the Netherlands. Kearney. R.E.. and Hunter. I.W. 1990. System identification of hu­man joint dynamics. Crit. Rev. Biomed. Eng. 18: 55-87. Korenberg. M.J.. and Hunter. I.W. 1999. Two methods for identify­ing Wiener cascades having noninvertible static nonlinearities. Ann. Biomed. Eng. 27: 793-804. Lee. D.. Port. N.L.. Kruse. W.. and Georgopoulos. A.P. 1998. Vari­ability and correlated noise in the discharge of neurons in motor and parietal areas of the primate cortex. J. Neurosci. 18: 1161- 1170. Loeb. G.E.. Bak. M.J.. and Duysens. J. 1977. Long term unit re­cording from somatosensory neurons in the spinal ganglia of the freely walking cat. Science (Washington. D.C.). 197: 1192- 1194. Loeb. G.E.. Hoffer. J.A.. and Pratt. C.A. 1985. Activity of spindle afferents from cat anterior thigh muscles. I. Identification and patterns during normal locomotion. J Neurophysiol. 54: 549­564. Matthews. PB.C. 1972. Mammalian muscle receptors and their central actions. Arnold. London. Matthews. P.B.C.. and Stein. R.B. 1969. The sensitivity of muscle spindle afferents to small sinusoidal changes in length. J. Physiol. 200: 723-743. Maynard. E.M.. Hatsopoulos. N.G.. Ojakangas. C.L.. Acuna. B.D.. Sanes. J.N.. Normann. R.A.. and Donoghue. J.P. 1999. Neuronal interactions improve cortical population coding of movement di­rection. J. Neurosci. 19: 8 083-8 093. Mussa-Ivaldi. F.A. 1988. Do neurons in the motor cortex encode movement direction? An alternative hypothesis. Neurosci. Lett. 91: 196-111. Nicolelis. M.A.L. 2001. Actions from thoughts. Nature (Lond.h 409: 403-407. Otten. E.. Hulliger. M.. and Scheepstra. K.A. 1995. A model study on the influence of a slowly activating potassium conductance on repetitive firing patterns of muscle spindle primary endings. J. Theor. Biol. 173(1): 67-78. Perkel. D.H.. and Bullock. T.H. 1968. Neural coding. Neurosci. Res. Prog. Bull. 6: 221-348. Prochazka. A. 1989. Sensorimotor gain control: a basic strategy of motor systems? Prog. Neurobiol. 33: 281-307. Prochazka. A., and Gorassini. M. 1998a. Ensemble firing of mus­cle spindle afferents recorded during normal locomotion in cats. J. Physiol. 507: 293-304. Prochazka. A., and Gorassini. M. 1998fc. Models of ensemble fir­ing of muscle spindle afferents recorded during normal locomo­tion in cats. J. Physiol. 507: 277-291. Prochazka. A.. Westerman. R.A.. and Ziccone. S.P. 1976. Dis­charges of single hindlimb afferents in the freely moving cat. J. Neurophysiol. 39: 1090-1104. Prochazka. A.. Weber. D.J.. Stein. R.B.. and Normann. R.A. 2003. Global kinematics can be deduced from the firing of small num­bers of sensory afferents during normal locomotion. Soc. Neurosci. Abst.: 267.7. Richmond. B.J.. Optican. L.M.. Podell. M.. and Spitzer. H. 1987. Temporal encoding of two-dimensional patterns by single units © 2 0 0 4 N RC Canada in primate inferior temporal cortex. I. Response characteristics. J. Neurophysiol. 57: 132-146. Rieke, F„ Warland, D., de Ruyter van Steveninck, R.R., and Bialek, W. 1997. Spikes: exploring the neural code. MIT Press, Cambridge, Mass. Rothwell, J.C. 1994. Control of human voluntary movement. Chap­man & Hall, London. Rousche, P.J., and Normann, R.A. 1992. A method for pneumati­cally inserting an array of penetrating electrodes into cortical tissue. Ann. Biomed. Eng. 20: 413-422. Schwartz, A.B. 1992. Motor cortical activity during drawing move­ments: single-unit activity during sinusoid tracing. J. Neurophysiol. 68: 528-541. Schwartz, A.B. 1994. Direct cortical representation of drawing. Science (Washington, D.C.), 265: 540-542. Scott, S.H., and Kalaska, J.F. 1997. Reaching movements with similar hand paths but different arm orientations. I. Activity of individual cells in motor cortex. J. Neurophysiol. 77: 826-852. Scott, S.H., and Norman, K.E. 2003. Computational approaches to motor control and their potential role for interpreting motor dys­function. Curr. Opin. Neurobiol. 16: 693-698. Serruya, M.D., Hatsopoulos, N.G., Paninski, L., Fellows, M.R., and Donoghue, J.P 2002. Instant neural control of a movement signal. Nature (Lond.), 416: 141-142. Shadlen, M.N., and Newsome, W.T. 1998. The variable discharge of cortical neurons: implications for connectivity, computation and information coding. J. Neurosci. 18: 3 870-3 896. Shoham, S. 2001. Advances towards an implantable motor cortical interface. University of Utah, Salt Lake City, Utah. Shoham, S., Fellows, M.R., and Normann, R.A. 2003. Robust auto­matic spike sorting using mixtures of multivariate t- distributions. J. Neurosci. Methods, 127: 111-122. Stein, R.B. 1965. A theoretical analysis of neuronal variability. Biophys. J. 5: 173-194. Stein, R.B. 1967. The information capacity of nerve cells using a frequency code. Biophys. J. 7: 797-826. Stein, R.B. 1970. The role of spike trains in transmitting and dis­torting sensory signals. In The Neurosciences, Second Study Program. Rockefeller University Press, New York. pp. 597-604. Stein, R.B., and Weber, D.J. 2004. Editing trains of action poten­tials from multi-electrode arrays. J. Neurosci. Methods. 134: 91­100. Stevens, S.S. 1961. The psychophysics of sensory function. In Sen­sory communication. MIT Press, Cambridge Mass. pp. 1-33. Taylor, D.M., Helms Tillery, S.I., and Schwartz, A.B. 2002. Direct cortical control of 3D neuroprosthetic devices. Science (Wash­ington, DC), 296: 1829-1832. Terzuolo, C.A., and Washizu, Y. 1962. Relation between stimulus strength, generator potential and impulse frequency in stretch re­ceptor of crustacea. J. Neurophysiol. 25: 56-66. Todorov, E. 2000. Direct cortical control of muscle activation in voluntary arm movements: a model. Nat. Neurosci. 3: 391-398. Weber, D.J., Stein, R.B., Aoyagi, Y„ and Prochazka, A. 2002. Chronic multi-unit recording of sensory neuronal activity in the cat dorsal root ganglion. Soc. Neurosci. Abst. 348.7. Werner, G., and Mountcastle, V.B. 1965. Neural activity in mechanoreceptive cutaneous afferents: stimulus-response rela­tions, Weber functions, and information transmission. J. Neurophysiol. 28: 359-397. Wessberg, J., Stambaugh, C.R., Kralik, J.D., Beck, P.D., Laubach, M., Chapin, J.K., Kim, J., Biggs, S.J., Srinivasan, M.A., and Nicolelis, M.A. 2000. Real-time prediction of hand trajectory by Stein et al. ensembles of cortical neurons in primates. Nature (Lond.), 408: 361-365. Zeki, S. 1993. A vision of the brain. Blackwell Scientific Pub­lishers, Oxford. Fig. 11. Predictions from a stochastic process (see Appendix A), in which it takes one (Poisson process), 2, 4, or 8 pulses (EPSPs) to produce an output (spike). The relationship between the variance and mean number of spikes is exact for the Poisson process, but was simulated for 2 (x), 4 (O), or 8 (□) pulses. The solid lines on the right are the known asymptotic behaviors, as the mean number of spikes becomes very large. The solid curves on the bottom are the calculated values if the process has no variation. 767 Appendix A This section will briefly present analytic and computed re­sults from a neural model that can be compared with the ex­perimental data. In particular, we will show that even a simple neural model predicts the cubic input-output func­tions and the dependence of the variance on the mean that was found in Fig. 7. Because the generation of nerve im­pulses can be considered to be a stochastic point process, certain predictions can be made directly from the theory of these processes (Cox and Miller 1965). The simplest exam­ple is the Poisson process. In this process, variance increases proportionally to the mean number of pulses (spikes) that occur in time t (see Fig. 11). At the other extreme is a neu­ron that fires completely regularly with a mean of / + d im­pulses in time r, where / is an integer and d is a fraction between zero and one. It immediately follows that the vari­ance increases to a maximum, when d = 0.5, and goes to zero when d = 0. Various intermediate examples can be con­sidered. For example, if the neuron integrates inputs per­fectly and it takes n EPSPs to produce a spike, then the interval distribution is given by a gamma distribution of the order of n. A gamma distribution of the order of 1 is the Poisson process considered above. Results for n = 2, 4, and 8 were calculated by selecting intervals at random from the probability density function for a gamma distribution (Cox © 2 0 0 4 N RC Canada 768 Can. J. Physiol. Pharmacol. Vol. 82, 2004 and Miller 1965). The point process was converted to a his­togram of the number of spikes in time t, where t = 0.1 s, as it was in our experiments. Finally, it can be shown that the variance will increase proportionally with the mean for large mean numbers, and that the slope of the line is 1/n (Cox and Miller 1965). Figure 11 also shows these asymptotic values. Although easily calculated, these values are for the steady-state situation, and may not apply to time-varying in­puts. Nor do they include the effects of nonlinearities, which are introduced by the refractory period of a neuron and the decay of EPSPs over time. The model was implemented with a Matlab program and includes 5 components. First, a band-limited input time series was generated by low-pass fil­tering uniformly distributed random numbers (third order Butterworth filter, 2 Hz cutoff). Second, the filtered signal was scaled and offset to provide different levels of stimuli. Third, the neural model consisted of a "leaky integrator," with a time constant of 0.1 s, when the output was less than 1.0. Fourth, when the output exceeded 1.0, a pulse (spike) was generated and the model was reset to zero for a refrac­tory period of 5 ms. Fifth, an independent source of noise (random numbers) was added to the neural model to simu­late the intrinsic variability observed in actual neurons. Figure 10 shows 20 s of the random input and a histogram of spikes that were produced. The best linear model was cal­culated using the least mean squares method, as it was for the experimental data, and the actual number of spikes was compared with the predictions. A cubic polynomial fitted the relationship well. The variance increased with the mean, but then began to decline as the actual number of spikes began to saturate. The data compare very well with those observed experimentally in Fig. 7. Clearly, there are many other non­linear features in neurons. For example, an after hyper­polarization was added to the model that could sum with repetitive activity, but this did not change the character of the responses. In conclusion, we suggest that the model pre­sented here may be the simplest one to describe the encod­ing mechanisms of DRG neurons. Reference Cox, D.R., and Miller, H.D. 1965. The theory of stochastic pro­cesses. Methuen & Co., London. © 2 0 0 4 N RC Canada