Abstract
The premotor cortex (PM) is known to be a site of visuo-somatosensory integration for the production of movement. We sought to better understand the ventral PM (PMv) by modeling its signal encoding in greater detail. Neuronal firing data was obtained from 110 PMv neurons in two male rhesus macaques executing four reach-grasp-manipulate tasks. We found that in the large majority of neurons (∼90%) the firing patterns across the four tasks could be explained by assuming that a high-dimensional position/configuration trajectory-like signal evolving ∼250 ms before movement was encoded within a multidimensional Gaussian field (MGF). Our findings are consistent with the possibility that PMv neurons process a visually specified reference command for the intended arm/hand position trajectory with respect to a proprioceptively or visually sensed initial configuration. The estimated MGF were (hyper) disc-like, such that each neuron's firing modulated strongly only with commands that evolved along a single direction within position/configuration space. Thus, many neurons appeared to be tuned to slices of this input signal space that as a collection appeared to well cover the space. The MGF encoding models appear to be consistent with the arm-referent, bell-shaped, visual target tuning curves and target selectivity patterns observed in PMV visual-motor neurons. These findings suggest that PMv may implement a lookup table-like mechanism that helps translate intended movement trajectory into time-varying patterns of activation in motor cortex and spinal cord. MGFs provide an improved nonlinear framework for potentially decoding visually specified, intended multijoint arm/hand trajectories well in advance of movement.
Introduction
The ventral premotor cortex (PMv) receives, among other signals, proprioceptive and visual information about arm/hand position/configuration as well as visual information about potential targets for reaching and grasping. It is active during movement preparation and execution of movement, as well as during the application of contact force (Table 1). Many PMv neurons exhibit movement-specific activity considerably before movement onset, often at the time of target identification or presentation. Therefore, PMv potentially provides signals useful for very early determination of motor intent. Understanding how information is encoded in PM has gained impetus from interest in extracting cortical information for controlling prosthetic limbs (Donoghue, 2002; Schwartz, 2004; Santhanam et al., 2006; Aggarwal, 2011). As decoding accuracy can be enhanced by quantitative understanding of both the signals encoded and the encoding process, faithful models of PMv neuronal firing should be of value.
Principal investigations characterizing PMv function
PMv is structurally heterogeneous and resides within a potentially complex regional sensorimotor control network (Fig. 1). Still, it is observed that during similar movements many neurons exhibit distinct, reproducible firing patterns. This consistency raises the possibility that PMv neuronal activity might be describable by a fixed function of movement-related signals. Although linear models that relate kinematic and force signals to neural firing patterns have enjoyed some success in motor cortex (see, e.g., Hepp-Reymond et al., 1999; Shoham et al., 2005), such models do not account well for the behavior of a significant fraction of the neurons. As an initial approach, we investigated whether much of the firing pattern variation might be explained by nonlinear encoding of simple signal content. Taking note of the previous finding of bell-shaped firing intensity tuning curves in relation to visual target location (Graziano et al., 1997), we examined Gaussian encoding in particular. And in light of the experimentally observed salience of target and limb location/configuration information in affecting PMv neuronal firing behavior (Raos et al., 2006), we focused first on kinematic rather than force-related signal content. Specifically, we tested the hypothesis that PMv neurons encode multidimensional kinematic signals within multidimensional Gaussian fields (MGFs). For comparison, we evaluated the performance of linear and two simpler nonlinear multidimensional models. MGFs are generalizations of one-dimensional Gaussian tuning curves (Wei and Stocker, 2012) and are analogous to receptive fields in primary sensory areas (Hubel and Wiesel, 1962; Levick, 1972). We examined 110 PMv neurons recorded in two primates that executed four reach-grasp-manipulate tasks and tested the dependence of firing pattern on four kinematic signals. We found that models that considered PMv neurons to encode time-advanced, limb position/configuration trajectory-like signals within rank-1 or full-rank MGF with at least five dimensions were significantly more efficient than the other models considered. This finding is consistent with the possibility that PMv function includes prominently the processing of a visually and proprioceptively specified, feedforward reference command for intended limb position/configuration.
Principal structure and regional connectivity of PMv. In the macaque, PMv consists of at least two subareas with somewhat different input sources: F4 lies caudally, adjacent to motor cortex; the F5 complex, consisting of subparts F5a, F5c, and F5p, lies rostrally. The human analogs are designated PMvc and PMvr, respectively. F5a receives preferentially visual information from the AIP and projects to F5c, F5p, and F4 (Gerbella et al., 2011). F4 also receives multimodal sensory input directly from the VIP (Fogassi et al., 1996; Luppino et al., 1999). Both F4 and F5 are also bidirectionally interconnected with the secondary somatosensory area (S2), the dorsal premotor cortex (PMd), the supplementary area (SMA), areas in prefrontal (PF), cingulate (CG) cortices, cerebellum (CBLLM), and basal ganglia (BG). Principal outputs are to brain stem (BS), spinal cord (SC), primary motor cortex (M1), and SMA (Dum and Strick, 2002; Dancause et al., 2006).
Materials and Methods
Experimental tasks and data collection.
The details of the experimental and recording methods were published previously (Mollazadeh et al., 2011). Briefly, two male rhesus monkeys (Macaca mulatta; M1 and M2) were trained to perform four reach-grasp-manipulate tasks (Fig. 2): (1) reach and pull a mallet; (2) reach and pull a rod; (3) reach and push a button; and (4) reach and rotate a ball (Fig. 2A,B). At the start of the trial, the primate held the centrally located home object. At a random interval between 330 and 1100 ms, a second visual cue identified one of four target objects located on a 13 cm radius circle around home. The monkey was to then reach, grasp, and manipulate the designated target object and to follow with a steady hold. Although the button was not “grasped” in the usual sense, an analogous particular hand configuration with curled fingers and extended thumb was consistently adopted before the push. The target object was held for 1000 ms to receive a food pellet reward. The trials were aborted if the monkey prematurely released the home object or failed to release the home object within 1000 ms. The monkeys could see both the target and their own arms throughout the movement. For M1, an average of 112 and 60 trials across tasks were conducted on days 1 and 2, respectively. Because we could not be sure that electrode positions remained constant between the two sessions, the neurons recorded during M1's two sessions were treated as being distinct from each other. For M2, 160 trials/task were completed in a single recording session. Both kinematic and neuronal data were recorded while the monkey performed these tasks. In total, 123 PMv neurons were recorded from M1 and M2 using floating microelectrode arrays (MicroProbes for Life Sciences). Each floating microelectrode array consisted of 16 parylene-C-insulated platinum/iridium recording microelectrodes arranged in a 4 × 4 triangular matrix on a 1.95 × 2.45 mm ceramic chip. The lengths of these electrodes vary from 1.0 to 6.0 mm in monkey M1 and from 1.5 to 8.0 mm in monkey M2. The recording arrays were located at the boundary of the F4 and F5 subregions of PMv. Therefore, some of the electrodes were located in F4, whereas others were in F5. Thirty sensors were placed on the monkey's hand that recorded their x, y, and z positions. The position of the sensors on the hand and their movement during 4 tasks are shown in Figure 2.
Experimental setting. In our experiment, the monkeys performed four reach-grasp-manipulate tasks: (1) reach and pull a mallet (perpendicular cylinder); (2) reach and pull a rod (coaxial cylinder); (3) reach and push a button; and (4) reach and rotate a ball. A, Position of objects. B, Position of objects with monkey's hand C, Position of sensors on the hand. The figure is modified and reproduced with permission from the following: Mollazadeh et al. (2011; their Fig. 1) and Aggarwal (2011; their Fig. 3.10).
Processing of neuronal recordings.
The neuronal spiking data consist of extracellularly recorded continuous voltage profiles. These were analyzed and spike sorted using the Offline Plexon Sorter (2003, Plexon; http://www.plexon.com) to yield spike times tk for each neuron. Neuronal spike trains with interspike intervals <1 ms were considered to be generated from more than one neuron (n = 13) and were removed from the dataset. This process enabled creation of idealized spike trains ψ(t) = Σkδ(t − tk), where δ(·) is the Dirac δ function that represents a single spike. We denote the spike train from neuron i for trial j on task l as ψijl(t). All trials were aligned in time so that time t = 0 represented the onset of movement. We noted that, even while motor performance remained quite similar from trial to trial within the same task, the neuronal firing intensity varied in a seemingly random manner. By computing running estimates of the mean and variance of the number of spikes using windows of length 5 ms (the bin size we chose for the study), verified the proportionality of these two parameters, indicating that we could reasonably view the spikes as being generated by time varying Poisson process. Therefore, we considered our neural spike trains to be the sum of two component inhomogeneous processes. The first is an inhomogeneous (time-varying rate) Poisson process with rate ζijl(t) that may vary from trial to trial. The second is an inhomogeneous Poisson process with rate ζi0l(t) that remains the same for all trials at a given task. Thus, the rate of the total process λijl(t) is given by Equation 1 (Cox and Isham, 1980; Leemis, 2003):
The first component is considered to correspond to background firing activity that may change according to moment to moment variations in the state of the neuronal network involving PMv, or to other variables not systematically related to the experimental tasks or trials. The second component is attributed to the neural processing that is systematically related to the experimental tasks. The simple additive partitioning was done for parsimony in the absence of more specific information about trial-by-trial firing variation.
Under these assumptions, we could construct virtual task-related neuronal spike trains ψil(t) for each neuron and task by summing ψijl(t) across trials as follows:
where the symbol ≜ specifies a definition. As all four tasks had slightly different number of trials, and data were divided into training and testing sets, n is half of the minimum (over tasks) number of trials done. Further, because of different time duration of trials, a common window (−0.75, 2 s) relative to movement onset was identified, and the data were averaged only inside this window. Of this 2.75 s window, a subwindow of 1.5 s was used for fitting the models, depending on an assumed delay value. ψil(t) should therefore be generated by a Poisson process with rate
In Equation 3b, we approximated ∑j=1nζijl(t) ≃ ci after verifying that there was no systematic relationship of background firing rate to time, task, or trial. As defined, virtual task-related neural spike trains evidently have a much higher rate than those of individual trials. This facilitates more accurate estimation of λil(t). The ratio of the rate λil(t) to its standard deviation, a form of signal-to-noise ratio, is
We then derived a task-related neuronal response signal ril(t) as follows:
where “*” represents the convolution operation and h(t) is the impulse response of a smoothing filter. ril(t) is our estimate of λil(t) for the i-th neuron in task l. It is obtained by low-pass filtering ψil(t) with a two pole Butterworth filter with cutoff frequency 30 Hz and then scaling the result by
Finally, we constructed a total neuronal response signal ri(t) ≜ [ri1(t), ri2(t), ri3(t), ri4(t)] by concatenating the signals ril(t) from all four tasks. This signal represented the complete behavior of neuron i's firing rate in the experiment. The distribution of average neuronal firing rates is shown in Figure 3.
Firing rate. Empirical distribution of mean firing rate across 110 neurons recorded from PMv of two monkeys during reach-grasp-manipulate task.
Processing of kinematic data and specification of candidate behavioral signals.
We noted that the monkey's arm trajectories were quite similar across trials of the same task (r2 > 0.92, for all tasks). Therefore, we averaged the recorded kinematic data over n trials to derive a task-related prototypical trajectory xl(t) for each task l:
where we denote the kinematic trajectory data for trial j and task l by xjl(t). To match the spike train data, xjl(t) from different trials was aligned such that xjl(0) represented the movement onset in trial j. Finally, we generated a total prototypical position/configuration trajectory vector as x(t) ≜ [x1(t), x2(t), x3(t), x4(t)] by concatenating the signals xl(t) from all four tasks.
x(t) was used to derive five types of prototypical behavioral signals as follows:
The 90-dimensional Cartesian position/configuration vector signal x(t): This vector contains the 3 d Cartesian coordinates of 30 sensors on the hand in workspace coordinates.
The 21-dimensional joint angle vector signal, q(t): In our experiment, the angles of 21 joints of the arm and hand were computed from the positions of the 30 hand sensors. The shoulder and elbow joint angles were estimated trigonometrically.
The 21-dimensional joint angular velocity vector signal, q̇(t): Numerically differentiated from joint angle data.
The 21-dimensional joint acceleration vector signal, q̈(t): Numerically differentiated from angular velocity data.
The 42-dimensional joint angle and joint angular velocity signal [q(t), q̇(t)]T.
The prototypical signals were then used as waveforms to specify a set S of five types of alternative behavioral signals s(t) : s(t) ≜ x(t − d) − xinitial, s(t) ≜ q(t − d)− qinitial, s(t) ≜ q̇(t − d), s(t) ≜ q̈(t − d), or s(t) ≜ [q(t − d) − qinitial, q̇(t − d)]T, where xinitial and qinitial were the arm/hand position/configuration at the home position and d is an arbitrary delay relative to the observed kinematics. We considered the possibilities that d could be positive or negative. Specifically, we examined delay/time advance values chosen from among 15 possibilities: d ∈ D = {−500, −450, …, 200 ms}. We selected these five types of signals as possible candidates for encoding within PMv neurons because position-like, velocity-like, and acceleration-like signals have been observed in motor cortex during movement (Georgopoulos et al., 1982, 1988; Paninski et al., 2004). The last signal allows the possibility that firing intensity could be driven by arbitrary linear combinations of joint angle and angular velocity. These five candidates specify only some of the possible waveforms for signals encoded in PMv. As we considered possible latencies d that would preclude interpretation of the behavioral signal as representing simple afferent feedback, any interpretation of the motor control role of these signals (e.g., as intended, predicted, or sensed kinematics) depends on determination of d and further analysis.
It is also important to note that the monkeys in this experiment executed four different tasks that collectively did not exploit the full range of motions possible at 21 joints. For example, none of the tasks involved spreading of the fingers or moving only the fourth digit. Therefore, the high-dimensional kinematic data are partially redundant. To reduce this redundancy, we performed principal components analysis (Jolliffe, 2002) separately on each of the five prototypical kinematic signals and found that the first five principal components (PC5) accounted for at least 85% of the total variance of each signal type. Therefore, we used the first five PCs as a five-dimensional vector signal s(t) ≜ PC5(s(t)) in place of s(t) to estimate the MGF more efficiently. In Figure 4, these are plotted for the position signal in Cartesian coordinates (s(t) = x(t) = PC5(x(t))). x(t) accounts for >95% of variance in x(t) across all tasks. The first three PCs predominantly capture translation of the hand in the x, z, and y directions, respectively. PC4 approximately captures wrist rotation and thumb movement, whereas PC5 approximately captures wrist rotation and finger movement. Because sensor movements during hand translation are larger than with wrist rotation or thumb and fingers movement, PC1 and PC2 account for most of the variance in x(t). However, some neuronal firing clearly depends most strongly on rotation and/or thumb and fingers movement. Therefore, after choosing first five PCs, we further normalized the variances so that the contributions of all PC were equal.
Principal components of change in position signal in Cartesian coordinates x(t). A, Cumulative variance percentage versus sorted principal components. The first five PCs account for nearly 95% variance in the data. B, Change in arm/hand position/configuration along the first five principal directions (moving from green to yellow).
Modeling neuronal firing intensity using sensorimotor fields.
As discussed before, we assumed that the task-related neuronal spiking activity can be described by an inhomogeneous Poisson process with rate, λil(t) (Cox and Isham, 1980). We also assumed that the task dependence of λil(t) obtains entirely through the encoding of task-dependent prototypical kinematic signals within time-, task- and trial-invariant generic neuronal sensorimotor fields λ(s(t), Θ). For each neuron, the field specifies firing intensity (rate) as a function of a behavioral input signal:
where Θi is the vector of fixed parameters that tailors the generic model sensorimotor field to neuron i. Accordingly, we represent the task-specific neuronal firing rate λil(t) as determined by the behavioral signal sl(t) as:
Nonlinear sensorimotor field classes.
Although linear encoding models have been used in primary motor cortex, the limitations of these models have been noted (Shoham et al., 2005; Aggarwal, 2011). Therefore, we evaluated the possibility that, for all task-related behavioral signals sl(t) (or sl(t), if a principal components representation is used), λil(t) belongs to one member of the following family (M) of six sensorimotor field classes as follows:
Linear fields
Square-root linear fields
Log-linear fields
Rank-1, MGF: parallel (hyper) plane
Full-rank, symmetric MGF: (hyper) spherical
Full-rank, general MGF: (hyper) ellipsoidal
Model classes 1–6 contain an increasing number of free parameters to be estimated. ci is the background neuronal firing rate from Equation 3b, that is taken to be independent of s(t). The second, nonlinear term in each of the models corresponds to ζi0l(t). In the second term, αi is an arbitrary scalar related to the offset of the field from the zero point of s(t), ki is the maximum additional firing rate due to s(t), σi is a scalar representing radial spread of the (hyper)spherical field and corresponds to the standard deviation of the Gaussian, bi is an arbitrary vector that specifies the directional orientation of field, μi represents the center of the Gaussian field with respect to the s(t) zero point, and Bi2 is a positive definite (all eigenvalues > 0), symmetric matrix.
For convenience, we took the initial, stationary arm/hand position/configuration s(0) (at home position) to be the zero vector. The spread of a rank-1 MGF can be described by σ =
The L, sqrt-L, and log-L fields are simple models used by several investigators for modeling neural firing in the primary motor cortex and the other brain areas (Sarma, et al., 2008, 2010, 2012; Kang et al., 2015). For PMv, however, these models are slightly unattractive a priori because they would not easily account for experimentally observed bell-shaped tuning curves. Still, they are reasonable candidates to evaluate before considering more complex models. By contrast, the MGFs have several attractive features for realistic modeling. First, they do not allow unlimited firing intensity as do linear, square root-linear, and log-linear fields. Second, they predict bell-shaped tuning curves for many inputs while they still allow for a range of possible realistic tuning properties discussed below. Finally, it is worth noting that, under certain fairly broad conditions, approximate rank-1 MGF encoding can be afforded by an elementary neural network (Fig. 5A). If neuron A receives convergent multichannel input represented by the vector X(t) and has a steady-state input–output firing probability relationship S̃(x), then when X(t) varies slowly with respect to internal neuronal and local network dynamics, the output of neuron A can be represented as S̃(bTX(t)), where b is the vector of input connection strengths that relates X(t) to neuron A. We will assume that S̃(x) is monotonically increasing as are most neuronal input–output relations (Fig. 5Bi). If each of the B neurons in isolation has approximately sigmoidal steady-state input–output firing probability relationships Si(βiu − θi) and if these neurons are also mutually inhibitory, then it is plausible that the output of the ith B neuron may be represented crudely approximately as follows:
where
Here, βi and θi are parameters that determine the peak slope and minimum effective threshold of Yi(u). γij is the inhibitory strength of the jth B neuron on the ith and Ii(u) is the total inhibitory input into neuron i. It has been noted in area M1 (Fromm and Evarts, 1981) that larger layer 5 pyramidal neurons tend to have somewhat higher input thresholds and outputs that saturate more gradually with increasing input. Essentially, in primary motor cortex there appears to be at least a weak “size principle.” We consider the possibility that the same is true in PMv. In Figure 5A, we let B1, B2, and B3 represent successively larger neurons. Under the additional asymmetry condition that larger neurons inhibit smaller neighbors more strongly while being inhibited less strongly by them (i.e., here γij > γji when j > i), then it may be verified by numerical solution of Equation 14 that for many values of βi, θi, and γij, Y1(S̃(bTX(t))) and Y2(S̃(bTX(t))) are bell-shaped functions of bTX(t), whereas Y3(S̃(bTX(t))) is sigmoidal, similar to one side of a Gaussian (e.g., Fig. 5Bii). Thus, Y1(S̃(x)) and Y2(S̃(x)) approximate rank-1 MGF encoding of X(t), and Y3(S̃(x)) approximates rank-1 MGF encoding before it saturates. The foregoing analysis is very limited, especially as it neglects fast neuronal and network dynamics. However, it supports the general feasibility of rank-1 MGF encoding to the extent that a neuronal network meets the relatively modest conditions given. In this regard, it may be relevant that, in areas F5c and F5p, layer 5 pyramidal cells are present in two or three distinct sizes, respectively (Belmalih et al., 2009) and that layer 5 pyramidal neurons are seen to participate in mutually inhibitory (among other types of) networks in somatosensory, frontal and visual cortices (Deuchars and Thompson, 1995; Kawaguchi and Kubota, 1997; Xiang et al., 1998).
Possible neural implementation of approximate rank-1 MGF encoding. A, Hypothetical neural network consisting of an input neuron A that receives converging signals, collectively X(t). Neuron A is assumed to have a sigmoidal steady-state response function S̃(x) that is substantially linear over a sizeable range and to excite neurons B1, B2, and B3. The latter are assumed to be mutually inhibitory with a size principle. Bi, S̃(x) is taken to be 10 × tanh(0.15x), which gradually saturates as x approaches 10. Si(βiu − θi) are of form (1 + tanh(βiu − θ i))/2. This function is shown for βi = 1, θi = 5 after scaling by 10 for clarity. Bii, Bell-shaped and sigmoidal net output responses from the B neurons are Yi(S̃(x)) for parameter values βi = 1.0, 0.8, 0.6; θi = 2.5, 4.0, 4.2; and γ12 = 1.0, γ13 = 1.5, γ21 = 0.002, γ23 = 0.85, γ31 = 0.001, γ32 = 0.08. The peak values, locations, and widths of the Yi(S̃(x)) vary with different values of connection strength parameters βi, θi, and γij.
The characterizations, “parallel (hyper)plane,” “(hyper)spherical,” and “(hyper)ellipsoidal” refer to the surfaces of constant response. Gaussian response profiles are predicted whenever input signals cut across these regions. The difference in the geometries of the three MGF classes can be appreciated as follows. First, Equation 12 is merely a uniformly symmetric special case of Equation 13 that has significantly fewer defining parameters. In Equations 12 and 13, the peak firing rate parameter value ki + ci occurs whenever s(t) = μi and falls off as a Gaussian curve with distance as s(t) diverges from the point μi along any straight line in input signal space. The regions where neuronal firing rate remains greater than any arbitrary value λ0 (lying between ci and ki + ci) are those where s(t) remains within, respectively, some (hyper)sphere or (hyper) ellipsoid around μi. The neuron could be considered to be “tuned” or sensitive to signals in these localized regions. In contrast, Equation 11 is a nonuniformly symmetric, limiting case of Equation 13 that also has significantly fewer defining parameters. In Equation 11, the peak value of firing rate occurs wherever s(t) = μib̃i+p, with
and p is any vector perpendicular to bi. It may be noted that μib̃i is a point in the direction of bi, with distance μi from the origin. This means that peak firing occurs along one or more lines in signal space rather than at a single point. In multidimensional space, there are many possible directions and magnitudes for the vector p. Therefore, s(t) may deviate perpendicularly from bi in many different directions while yielding a maximal (or some other constant level of) neuronal firing. For Equation 11, the regions where neuronal firing remains greater than some value λ0 (between ci and ki + ci) are infinitely wide, slab-like regions that lie between two parallel (hyper)planes perpendicular to bi. Conceptually, these regions can be viewed as the limit of large, flat, disc-like ellipsoids that are much narrower along one direction in input space than along all other directions. When s(t) diverges from μib̃i in any direction that is not perpendicular to bi, firing rate again falls off as a Gaussian curve with distance as with Equations 13 and 12. Firing is predicted to modulate most strongly when s(t) varies along the direction of b.
The qualitative features of rank-1 MGF structure can be appreciated most easily in 2 dimensions where the field λ(s, Θ) appears as an infinitely long ridge on the input signal plane (Fig. 6, blue surface). The regions of high firing rate are strips in the s plane that lie between two lines perpendicular to b on either side of μb̃. Figure 6A shows the relationship between several possible simple trajectories for s(t) and a rank-1 MGF defined by b and μ. Our behavioral input signals all began at the origin and radiated outward in PC space, although not generally along straight lines. Those signals related to position traveled to endpoints different from the origin, whereas those related to velocity or acceleration returned to the origin at the end of movement. For simplicity, in this figure we have used straight, radial trajectories from the origin to represent possible s(t). These are most similar to observed position/configuration related signals x(t) and q(t). For the two rightmost trajectories, maxt(b̃Ts(t)) > μ. In this case, the MGF will result in significant nonlinear modulation of firing rate over the course of trajectory. For the middle two trajectories, firing will be nontrivial but constant as they travel perpendicular to b. With respect to these first four trajectories, the MGF can be considered “local” as it gives rise to significant firing. On the other hand, the two leftmost trajectories will be associated with much less firing because they lead away from the MGF. Here, approximately maxt(b̃Ts(t)) < μ − σ and the MGF can be considered relatively “remote.” As to be discussed further below, Figure 6B shows a potentially important special “near linear” case of a local MGF where the full set of radial trajectories lies close to one side of a rank-1 MGF. This occurs approximately when μ − 2σ < mint(b̃Ts(t)) < maxt(b̃Ts(t)) < μ. In this case, the firing rate of the neuron will depend largely on the projection of the trajectory on the direction of b̃. For radial trajectories, this projection will be proportional to the cosine of the angle between the trajectory and b. Cosine tuning of neuronal firing with respect to movement direction has been observed experimentally in motor cortex by several investigators (Schwartz et al., 1988; Moran and Schwartz, 1999). Figure 6C shows an MGF that can be considered “remote” with respect to the input because all s(t) trajectories lie far away from the MGF center. As a result, the neuron fires infrequently and modulates little for all of the modest-sized inputs shown. It is evidently tuned most specifically to large signals in the b direction.
Types of rank-1 MGF relative to input in two dimensions. A, Local field, the maximum extent of some input trajectories lies significantly beyond the center of the MGF. B, Near linear field, the extent of the trajectories is small compared with SD of MGF, and the relationship is such that the field is approximately planar. C, Remote field, the entire set of input trajectories lies far away from MGF center.
Finally, it is useful to note that certain full-rank, positive definite B2 matrices are nearly rank-1 because they have a single eigenvalue whose magnitude is much larger than those of all the others. Such matrices define fields with hyper-ellipsoid tuning regions that are flat, very wide, yet finite (pancake-like) disks that approximate the infinite slab-like tuning regions of Equation 11. In 2 dimensions, such regions are long, narrow ellipses that near the center approximate two parallel lines. Accordingly, these full-rank fields may be indistinguishable from rank-1 fields in practice.
Model parameter estimation.
We estimated the parameters Θi of different sensorimotor fields using maximum likelihood estimation. One can show that, for multitrial data with a large number of trials n, the neuronal the firing rate ri(t) (Eq. 4) is approximately normally distributed with mean
The above maximization problem has several local maxima. Therefore, we used stochastic maximization techniques (Spall, 2003) to obtain a globally optimal estimate Θ̂ for each neuron. Also, during our search, we constrained the possible parameter vector estimates as Θa ≤ Θ̂ ≤ Θb for some fixed Θa, Θb to ensure that our search results fell within a biologically feasible range (Table 2). For instance, we included the constraint 0 < ki < 100 Hz as neurons do not typically fire above 100 Hz.
Parameter bounds
The remaining model parameter d was estimated in a slightly different manner because of the possibility that its value differed between subpopulations in PMv rather than between individual neurons. In principle, delays might be due to mechanisms within PMv and therefore might vary significantly between each PMv neuron. However, in the absence of a specific putative mechanism, we considered this possibility to be unlikely. Rather, given that PMv receives inputs principally from anterior intraparietal area (AIP), SII, and perhaps cerebellum, we evaluated the possibility that the input signals might be present in up to three different PMv neuronal subpopulations with significantly different delay values. To do this, we considered each of the 15 values of d ∈ D to specify a different model. We then sought to determine whether the assumption of a mixture of models was superior to assuming just one (below). We acknowledge that the signals from different input sources might also differ in more ways than just delay and could be present in different combinations. However, in light of the apparently simple tuning of PMv neurons in relation to potential prehension target location (Graziano et al., 1997, 2002), we did not examine these more complicated possibilities here.
Model selection.
Because our model classes did not have the same number of free parameters, we selected models using the Akaike Information Criterion (AIC) (Akaike, 1974) on twofold cross-validation spiking data (Eq. 16). The model with the lowest AIC value is considered the better model. The AIC penalizes a larger number of parameters in a model unless it results in a sufficiently greater model likelihood. In this sense, low AIC values reflect greater model efficiency in representing the data. To compute AIC values, we first formed the total neuronal spike train ψi(t) = [ψi1 (t), ψi2 (t), ψi3 (t), ψi4 (t)] for each neuron using the testing dataset, and then found the total predicted firing intensity across trials λ̂i(t) = λ(s(t), Θ̂i) for each neuron using 450 (6 × 5 × 15) candidate sensorimotor field models encompassing from all model classes (M), kinematic signal types (S), and delay values (D). We then computed the likelihood Li that ψi(t) was generated by an inhomogeneous Poisson process with rate λi(t). This yielded the AIC values for each neuron and model using both training (a) and testing (e) data.
A preliminary analysis of the training dataset demonstrated that for each model m and neuron i certain signals sm,ia and delay values dm,ia produced the lowest AICm,s,d,ia relative to all other signals and delays. Using this identified optimal signal and delay for each neuron, it was determined that among all model classes, FR-MGF produced the lowest mean AIC value across all neurons i in the testing set. Therefore, the FR-MGF model class was considered to be the nominally most efficient. Further analysis was directed toward determining the significance of any difference between FR-MGF and the other model classes. For this purpose, ΔAIC values were computed for each class and each neuron as the following difference:
ΔAICm,ie is positive for neuron i whenever the FR-MGF model class is more efficient than model class m ≠ FR-MGF on the test data. It is negative when m is less efficient. The empirical population mean Δ̅A̅I̅C̅me =
After choosing model class, the prototypical signal that works best for PMv neurons was selected. Briefly, first ds,ia ≜ argmind∈D AICm=m̂,s,d,ia was estimated and ΔAICs,ie ≜ ΔAICm=m̂,s,d=ds,ie,ia − AICm = m̂,s=pos-x,d=ds=pos-x,ia,ie values were calculated for each signal and neuron. Then the mean Δ̅A̅I̅C̅se =
After choosing both model and signal class, the delay or delays that works best for PMv neurons were determined. As discussed above, we entertained the possibility that s(t) might be encoded with up to three different delay values. Therefore, PMv neurons were evaluated as consisting of one of the following: (1) one group of homogeneous neurons; (2) two subgroups; or (3) three subgroups of neurons based upon model performance across different delays. Presumably, if PMv consists of two or three subpopulations of neurons having two or three different signal delays, then when the whole population is modeled using any single delay value, two or three groups of neurons should be observed based on quality of fit. Those neurons having signal delay values closest to the value assumed for the model should be fit best; those with delay values farthest away from the model should be fit least well. Accordingly, two or three corresponding subpopulations of AIC values should be observed having different means. More generally, we can ask whether or not AIC analysis supports the assumption of two or three subpopulations, each having a different pattern of means across a range of delays. To this end, individual neuronal AICm=m̂,s=ŝ,d,ia,e values were computed for each neuron using the 15 different delay values. This yielded two 15 element vectors of AIC values for each neuron for training and test set, respectively. This vector summarizes the performance of each neuron for any delay. Using this vector, first, delay-insensitive neurons were identified using k-means clustering (Hartigan and Wong, 1979), as the group of neurons that show almost equal AIC values across the 15 delays on both training and testing set, and were discarded. The remaining neurons were considered delay sensitive. Among these neurons, the neurons that have similar signal delay should show a similarly oriented vector of AIC values. Therefore, PMv population is divided into 1, 2, or 3 groups based upon the orientation of these vectors again using k-means clustering (Hartigan and Wong, 1979) using training data. After finding 1, 2, or 3 different empirical groups of neurons, we computed the means of AIC values within each subgroup at each of the 15 delay values.
Here, kj is the number of neurons in jth subgroup. The delay that yielded the minimum A̅I̅C̅sga(d) was considered optimal for the subgroup j. This process yielded 6 different optimal delays (1 when considering only one group, 2 when considering two groups, and 3 when considering three groups). Using the 1, 2, or 3 groupings with 1, 2, or 3 group optimal delays, we recomputed the single “delay-tailored” mean ΔAIC for assuming 1, 2, or 3 groupings with respect to assuming 1 grouping across the entire population using test data. For instance, let's assume in case when we assume 3 groupings, k-means yielded groupings with k1, k2, and k3 neurons with d1, d2, and d3 optimal delays, the “delay-tailored” mean ΔAIC was computed as follows:
The subgrouping and the delays that best work for PMv neurons were selected by analyzing “delay tailored” mean ΔAIC calculated separately for each subgrouping.
Model performance evaluation.
After selecting the sensorimotor field type, signal type, and delay, the absolute performance of our models was evaluated using the time rescaling theorem and the KS-statistic for point process models (Brown et al., 2002). In our study, this statistic determines the confidence with which an observed spike train can be said to be generated from an inhomogeneous Poisson process with rate λ(s(t), Θ̂i). A KS-statistic normalized by 1.63/sqrt(n) (NKS) that is >1 indicates that the chance that the spike train was generated from such a process is <1%. Similarly, a NKS >0.83 indicates this chance is <5%. Therefore, for a model of PMv to be considered statistically indistinguishable from a correct model, 99% of neuronal spike trains must have NKS < 1. A model of PMv firing for which <95% of neuronal spike trains have NKS < 0.83, or for which <99% have NKS < 1, can be therefore considered significantly different from a correct model. We computed the percentage of neurons with NKS < 1 to give us a sense of how close our selected field models are to one that is theoretically correct. Still, it is possible for two models to be equally correct across tasks in a mean sense, whereas one of the two may be more precise in terms of less variance from the correct average predicted neuronal activity. Therefore, we also assessed the performance of the MGF models by examining visually the fidelity with which neuronal firing patterns are predicted by λ(sl(t), Θ̂i) and by attempting to reconstruct sl(t) from ril(t) as described next.
Kinematic reconstruction using nonlinear sensorimotor field models.
We also sought to determine whether the selected nonlinear field model(s) could be used to reconstruct behavioral signals from neuronal firing datasets. Essentially, this amounts to seeking to invert a model from Equations 8–12 to solve for sl(t) when using ril(t) to estimate λil(t). However, any observed ril(t) is only probabilistically related to λil(t) as described before. Moreover, Equations 8–12 are nonlinear and high-dimensional such that there are in general more than one sl(t) that yield a given λil(t). Therefore, we approached the inversion indirectly by first computing the likelihoods of observing ril(t) at time t for a fixed Θ̂i assuming different values for sl(t). As the likelihood function may have several local maxima, the likelihood was calculated over the set of all possible sl(t) values in the training set across all tasks. We then chose as our estimate of sl(t), the value ŝl(t) that provided the maximum likelihood over this set. This entire process was repeated at every time point t at intervals of 5 ms. Once the maximum likelihood estimate of s(t) was obtained, it was smoothed bidirectionally using a low-pass Butterworth filter with 10 Hz cutoff frequency.
Results
We first showed that rank-1 MGF and full-rank MGF are the most efficient model classes for representing PMv neuronal firing patterns. Second, behavioral signals x(t) and q(t) outperformed q̇(t), q̈(t) and [q(t), q̇(t)]T in rank-1 MGF model fitting. Third, kinematic data appeared to be represented in the PMv neuronal population primarily at one delay value of approximately −250 ms. Fourth, the fitted MGF appeared to cover the multidimensional behavioral signal space with sufficient overlap between individual MGF such that a number of PMv neurons would be activated simultaneously for most values of sl(t). Finally, we were able to reconstruct kinematic trajectories from cross-validation firing data with reasonable fidelity using both rank-1 and full-rank MGF.
Performance of different sensorimotor field classes
Figure 7A–C plots the empirical cumulative distribution of Δ̅A̅I̅C̅m,ie for all model classes m ∈ M over all neurons i = 1, …, k, and Δ̅A̅I̅C̅me for all model classes m ∈ M. It is seen, given the test data, that both rank-1 and full-rank MGF model classes (R1-MGF and FR-MGF, respectively) are significantly more likely to be correct model classes than are the L, sqrt-L, log-L, and Sph-MGF model classes. There is no significant difference in ΔAIC between the rank-1 and full-rank MGF model classes. These findings indicate that, at least for this dataset, R1-MGF models (8 parameters) are comparable with FR-MGF models (22 parameters) in data modeling efficiency. This suggests that any improvement in modeling seen for FR-MGF over that for R1-MGF is appropriate for the increased number of parameters of FR-MGF. On the other hand, Sph-MGF models appear to be inherently much less efficient. Together, these observations are strongly consistent with the likelihood that PMv fields are indeed (hyper)disc-like (or narrow ellipse-like for 2 dimensions) in shape. That is, they appear to be finite in size and much wider in all directions than in one. As discussed above and verified here, R1-MGF models therefore provide effective lower dimensional approximations to (hyper)disc-like (hyper)ellipsoids. For ease of evaluation and presentation, we selected R1-MGF for most further analyses.
Relative performance of different model classes. Ai–Ci, Distribution of ΔAIC over all PMv neurons calculated using cross-validation data for both rank-1 and full-rank MGF. The delay and signal type used to calculate the ΔAIC are the ones that give minimum ΔAIC on training data for each neuron. Aii–Cii, Mean ΔAIC for the six sensorimotor fields calculated overall all PMv neurons. A, Results for both Monkey M1 and M2 combined. B, Monkey M1. C, Monkey M2. Error bars show 95% confidence bounds. *p < 0.05.
Performance of different kinematic signals
Figure 8A–C plots the empirical cumulative distribution of Δ̅A̅I̅C̅s,ie for all signal classes s ∈ S over all neurons i = 1, …, k, and Δ̅A̅I̅C̅se for all signal classes s ∈ S. It can be appreciated that both position/configuration x(t) and joint position/configuration q(t) signals are significantly more likely to be present than are q̇(t), q̈(t), and probably [q(t), q̇(t)]. With the latter signal, we found a significant difference in one monkey, but not in the other. And in general, models using [q(t), q̇(t)] were nearly as efficient as those including x(t) or q(t) alone. We view this as being consistent with the possible presence of a small velocity-related signal, or with a larger velocity-related signal in some subset of PMv neurons, in addition to x(t) or q(t). We found no statistically significant difference between the modeling provided by x(t) and q(t). This is consistent with the possibility that positions in Cartesian and joint space in these tasks could be related by substantially linear transformations. To the extent that our MGF indicated neuronal sensitivity to relatively limited regions in both Cartesian and joint space, linear transformations are likely to be often adequate. In such a situation, the change in coordinate frame merely yields different optimal values for parameters αi, μi, bi, and Bi. Because we had direct measurement of x(t), we used it in preference to q(t) for all further evaluation. Still, we expect that using q(t) would have provided comparable results.
Relative performance of different signal types. Ai–Ci, Distribution of ΔAIC over all PMv neurons for best delay and for each position/configuration-x, position/configuration-q, velocity, and acceleration behavioral signal. Aii–Cii, Mean ΔAIC for the four behavioral signals calculated overall all PMv neurons. A, Results for Monkeys M1 and M2 combined. B, Results for Monkey M1. C, Results for Monkey M2. Error bars show 95% confidence intervals. *p < 0.05.
Performance using different delays
In Figure 9A–Ci, subgroup mean A̅I̅C̅sge values are plotted under each subgrouping for both the monkeys combined, M1 and M2, respectively. All three datasets demonstrate that, when assuming one homogeneous PMv population, the best (minimum) mean A̅I̅C̅sge values were obtained by assuming a delay value of −250 ms. When assuming two or three subpopulations within PMv, the trends in the mean A̅I̅C̅sge values generally showed broader delay ranges for lowest AIC values. Still, for at least one subpopulation, −250 ms is always within the range of delay values associated with near minimal A̅I̅C̅sge values. Figure 9A–Cii summarizes the distribution of the delay-tailored population A̅I̅C̅dt values for both monkeys combined (M1 and M2, respectively). We found that supposing that PMv consists of one homogeneous population with respect to delay gave a A̅I̅C̅dte value statistically equivalent to assuming two or three subpopulation within PMv. Hence, there was no statistical evidence to support or reject the assertion that there is more than one delay at which x(t) is encoded within PMv. On the other hand, the findings were consistent with the prominent presence of at least one signal of form x(t − d) that occurs at ∼−250 ms relative to movement.
Relative performance using different signal delays. Ai–Ci, Means of single neuron AIC values over PMv neurons, assuming 1, 2, and 3 subpopulations within PMv for both monkeys combined (Monkeys M1 and M2, respectively). Aii–Cii, Summary of distributions of delay-tailored population AIC means calculated assuming 1, 2, and 3 different subpopulations within PMv and choosing the best delays for the each subpopulation. Error bars show 95% confidence intervals.
Model fits of neuronal activity
Figure 10A, B illustrates differences in L, log-L fields, R1-MGF, and FR-MGF estimated for two example PMv neurons from Monkey M1 and Monkey M2, respectively, while the monkeys performed the four different reach-grasp-manipulate tasks. We found that, in 75% of neurons, firing was driven mainly by the field as ki was at least as twice large as ci. The two neurons whose activities are shown next are in this category. For 88% of neurons, the contributions of field-related and background firing were at least comparable (ki at least as large as ci). For the remaining 12%, the contribution of background activity was greater than that of the field. Subpanels i, ii, and iii plot vertically the value of the neuron's estimated linear fields, log-linear fields, and R1-MGF λ(s, Θ̂i), respectively, at all points s within an arbitrarily selected slicing plane (shown as the x, y plane) that contains each respective field's b vector (data not shown) directed along the x-axis as in Figure 6. Panel iv plots vertically the value of the neuron's estimated FR-MGF, λ(s, Θ̂i), at all points s within a slicing plane that contains first two eigenvectors of Bi directed along the x and y axes, respectively. In each case, the slicing plane provides a two-dimensional cross-section of the five-dimensional sensorimotor field. Because the behaviors of L, log-L fields, and R1-MGF in the three dimensions not depicted are identical to that shown along the y-axis, panels i and ii can be used to completely understand the characters of these fields. However, as the widths of the FR-MGF in the dimensions defined by the third, fourth, and fifth eigenvectors are not in general identical to that shown along the second eigenvector (along the y-axis), panel iv cannot depict the character of FR-MGF fully. Still, the panel confirms that, like the R1-MGF, the FR-MGF is narrow in one direction and much wider along a second dimension. Although it is not shown, it may be noted that the estimated FR-MGF becomes progressively much wider in the remaining 3 dimensions associated with progressively smaller eigenvalues.
Comparison of L fields, log-L fields, R1-MGF, and FR-MGF and their predicted firing responses for two example neurons (A, B) from Monkeys M1 and M2, respectively. Evaluation of (i) L fields, (ii) log-L fields, and (iii) R1-MGF above a slicing plane that contains each field's bi vector. iv, Evaluation of FR-MGF above the slicing plane that contains the first two eigenvectors of Bi. i, ii, iii, iv, The firing rates predicted by respective models, λ(sl(t), Θ̂i), are indicated on the fields for the four tasks (brown, green, magenta, and black lines) above the projections (data not shown) of sl(t) on the slicing plane. v, vi, vii, viii, λ(sl(t), Θ̂i) is overlaid on ril(t) measured during the four tasks.
Superimposed on each field are four curves showing the paths of the predicted task-related firing rates, λil(t). Each predicted firing rate curve in principle lies above the projection (data not shown) onto the slicing plane of its corresponding task signal sl(t). These task-related signal projections radiate from the origin and are analogous to the red radial trajectories in Figure 6. As can be seen from Figure 6A, B, both R1-MGF and FR-MGF have a maximum predicted rate and are compact. In particular, they are “local” with respect to the input signals in the sense of Figure 5. The neuron becomes much more active only when the input signal projection (data not shown) crosses the narrow strip-like/narrow-ellipse-like region under the ridge/mountain of either R1-MGF or FR-MGF. In contrast, L and log-L fields (Fig. 10 Ai,Aii,Bi,Bii) predicts unrealistically intense neural firing over a vast region of kinematic space. The evident local similarity of shape between the estimated R1-MGF and FR-MGF explains why these fields could be expected to often predict neural firing similarly. Panels v, vi, vii, and viii in Figure 10 plot the rates predicted by L, log-L fields, R1-MGF, and FR-MGF together with the empirical firing rate ril(t) for each of the four tasks. Time t = 0 marks movement onset. In the two neurons, respectively, the two R1-MGF (NKS = 0.48, 0.42, AIC = 4.4e3, 5.2e3) or FR-MGF (NKS = 0.41, 0.55, AIC = 4.2e3, 5.1e3) predict the observed firing of the neurons better than linear fields (NKS = 2.5, 0.68, AIC = 5.3e3, 5.5e3) and log-linear fields (NKS = 0.72, 0.58, AIC = 4.7e3, 5.5e3). The p value indicating the probability that R1-MGF and log-linear models do not differ in their ability to fit the data can be calculated as exp(ΔAIC/2), = exp(−108), exp(−128) ≪ 0.05. So, there is a highly significant difference in model class quality for both the neurons. In absolute terms, neurons whose firing patterns when modeled across all tasks generated a NKS < 0.55, as did the neurons shown in Figure 10, had at least a 40% chance of having an intensity function exactly as given by the model. For model classes FR-MGF, R1-MGF, log-L, and L, the percentage of neurons with NKS < 0.55 were 56%, 51%, 38%, and 32%, respectively. These percentages therefore approximate the fraction of neurons whose firing was modeled as well or better than those in Figure 10. A complementary analysis reveals that the r2 values for the fits in Figure 10 are as follows: 0.89,0.71; 0.8,0.52; 0.74,0.23; and 0.51,0.23 for FR-MGF, R1-MGF, log-L, and L models, respectively. For these four model classes, the percentages of neurons with r2 values at least great as 0.5 are 36%, 25%, 16%, and 3%, respectively. According to both analyses, the FR-MGF and R1-MGF were clearly superior in fitting the neuronal firing data. The fit provided by FR-MGF typically appears slightly better than that provided by R1-MGF. This is not surprising given its greater number of parameters. Finally, for the four model classes, the percentage of neurons with NKS < 1 were 93%, 89%, 79%, and 73%, respectively. Because each is <99%, each of the models is seen to be statistically significantly different from an ideal model. However, the FR-MGF and R1-MGF were closer to ideal.
Parameter sensitivity and specificity analysis
Once we decided which MGF model class was preferable, we examined the sensitivity of model performance in the neighborhood of the optimal parameter estimates. For this, we examined the effects of both 10% and 20% changes in each parameter value on the range and median value of the ΔAIC values. The latter corresponds to any horizontal shift in the 50% point of the CDF. Figure 11A, B plots the mean and median of ΔAIC values on changing the parameters by 10% and 20% across the neuronal population. We see that the fits are most sensitive to b1, b2, b3, and α and least sensitive to b4, b5, k, and c. This is consistent with a particular importance of field directional orientation (particularly with respect to the first 3 PC) and location with respect to the input signal.
Sensitivity of rank-1 models to change in parameters. Mean (A) and median (B) of empirical distribution.
Further, we found that while MGF representations that produced firing fits of at least the quality shown in Figure 10 could be found for at least half of PMv neurons, these representations as a whole were not strongly unique in up to 30% of neurons. Some parameters, such as the background firing level ci and maximum signal responsiveness ki, could be varied with relatively little loss of model likelihood. Still, others were found to be at least strongly optimal locally in the sense that as little as 10% and 20% change in especially the field location and orientation values yielded significant and progressive degradation of model likelihood (Fig. 11).
Location, size, and orientation of MGF
The location, size, and orientation of its MGF determine a neuron's sensitivity region within the input signal space. Figure 12A–C shows the distributions of μi and σi from Equation 15 for the 110 neurons from both monkeys: 61 neurons from Monkey M1 and 49 neurons from Monkey M2. Probability densities for the distribution for μi and σi were estimated using a kernel density estimation approach (Rosenblatt, 1958) with Gaussian kernels of width h = 0.5. We found that the range of μi values was severalfold larger than most values of σi. This indicates that, along their b vector directions, many R1-MGFs are narrow with respect to their distance from the origin. This suggests that Figure 6 provides realistic depictions of the relationships between R1-MGF width and their locations with respect to the origin of the input signal space. The MGF in our dataset were scattered along their b vector directions with little overlap with most of the other similarly oriented fields.
Distribution of μ and σ of 110 R1-MGFs. Ai–Ci, Empirically estimated probability density function of the distribution of R1-MGF field centers μ, along with μi for each neuron i (black dots). Aii–Cii, Probability density function for distribution of R1-MGF SDs σ, along with σi for each neuron i (black dots). A, Results for Monkeys M1 and M2 combined (110 MGFs). B, Results for Monkey M1 (61 MGFs). C, Results for Monkey M2 (49 MGFs).
As the long axes of a rank-1 MGF lie perpendicular to its b vector, the direction of the b vector determines the orientation of the MGF. In examining the distribution of bi vector directions (representable by the unit vectors b̃i), we found that in Monkey M1 the percentages of vectors having positive values in components 1–5 were 74%, 56%, 52%, 64%, and 44%, respectively. In Monkey M2, the percentages were 67%, 47%, 65%, 53%, and 55%, respectively. As most of these values are near 50%, the rank-1 MGF orientations were distributed within the five-dimensional input signal PC space with substantial uniformity.
To assess the possible functional significance of the observed distributions of μi, σi, and b̃i directions, we evaluated the R1-MGF in each monkey along the four radial vector input signal paths rj, j = 1, 2, 3, 4 that lie along the directions in five-dimensional PC space of the target position/configurations used in the four tasks. These paths are analogous to the red vectors in Figure 6 and correspond to the paths followed by trajectories modeled in Equation 21 below. The Gaussian functions shown in Figure 13 are cross sections of the MGF that can be described by parameters μsensij and σsensij. μsensij represents the location (i.e., possible input vector signal value) along path rj where neuron i will be maximally activated. σsensij indicates the spread of neuron i′ s region of sensitivity (or region of tuning) to signals that take values along path rj. These parameters are related to μi and σi, respectively, by the following:
where SEC(bi, rj) is the secant of the angle between the b vector of neuron i and rj. Evidently, the possible ranges for μsensij and σsensij are (μi, ∞) and (σi, ∞), respectively. In Figure 13, we see that the collection of R1-MGF provides sensitivity regions that cover behaviorally significant ranges of all four rj. Moreover, the overlapping sensitivity regions have multiple variances. Similar results were found in Monkey M1. We then chose four additional vector paths in random directions and repeated the analysis above. We found the same qualitative field distribution and coverage characteristics, which confirms the generality of these findings. If FR-MGF models had been used instead, we would have observed that many of the sensitivity regions would have lower peak intensities because of the fall-off in FR-MGF strength with distance from the field center in all directions. A number of the sensitivity cross-sections shown would become trivial in height. On the other hand, we expect that a larger sample of PMv neurons would provide denser and more highly overlapped coverage that would increase the number of nontrivial FR-MGF sensitivity regions. Together, these findings and considerations suggest that multiresolution coverage of signal paths by MGF is typical throughout the PMv input signal space. As a result, it appears likely that, whenever the multidimensional input signal changes value, a different unique collection of PMv neurons will be preferentially activated.
Signal path coverage of the R1-MGF in Monkey M2. A–D, Distributions of 49 sensitivity regions derived by evaluating 49 R1-MGF, respectively, along the target position/configuration radial vectors, rj, in the four tasks. When s(t) is taken to be x(t) as here, or q(t), the signal paths correspond to movements between arm/hand position/configuration in 5 d PC space. The starting and ending arm/hand position/configuration of these movements are depicted in Cartesian space as green and yellow hands, respectively. In each panel, the two small hand coordinate frames are located at the point along the x-axis that correspond to their positions in the rj vector direction. The green hand always represents the home arm/hand position/configuration that lies at the origin in the signal space. An 11 on the x-axis marks the target arm/hand position/configuration. For clarity, narrower (σ ≤ 4) and wider (σ > 4) sensitivity regions are plotted separately.
Kinematic reconstruction
In Figure 14, the time courses of the components of xl(t − d), l = 1, …, 4 are plotted for Monkey M2 along with the components of the maximum likelihood reconstruction of the behavioral signal x̂l(t − d) using fitted FR-MGF models and firing from all 49 neurons in the 80 cross-validation trials. Our reconstruction algorithm was blind to task type. We see that the first 5 PC of the kinematics were reconstructed with respectable fidelity (r2 = 0.99 for M2 and r2 = 0.94 for M1). We chose session 1 for reconstruction using M1 as this session had more trials (n = 56) and more neurons (n = 39) than session 2. The reconstruction accuracy using R1-MGF was very similar (r2 = 0.98 for M2 and r2 = 0.95 for M1), whereas the result using log-L fields was poorer (r2 = 0.88 for M2 and r2 = 0.84 for M1). Moreover, accuracy was more consistent across tasks when based on R1-MGF and FR-MGF than when based on log-L fields (Table 3). As may be deduced from Figure 6B, C, when trajectories are located strictly on one side of an R1-MGF, the field approximates a log-L field. Therefore, it may be expected that, for some subset of neurons, log-L fields may approach the performance of R1-MGF and hence of a similarly shaped FR-MGF on some tasks. The inability to generalize comparably across the space reflects the deficiency in nonlinearity of the log-L fields. L, sqrt-L, and Sph-MGF are geometrically even less similar to R1-MGF and disc-like FR-MGF than are log-L fields. As such, their reconstruction performance is uniformly poorer.
Reconstruction of kinematics. A–D, Reconstruction of five-dimensional kinematics using cross-validation firing rates for Monkey M2 for the four tasks. Dashed lines indicate the first five principal components of the actual kinematics time advanced by 250 ms, x(t + 250) − xinitial. Solid lines indicate the components of the five-dimensional reconstructed (estimated) behavioral signal ŝ(t) = x̂(t + 250) − x̂initial. Dotted black vertical lines approximately delimit the reach-and-grasp phases of the reach-grasp-manipulate signal separately. 0 indicates the onset of movement.
r2 values for reconstruction of x(t) for M2 and M1 using different sensorimotor field classes
Because x(t) carries >95% of the power (variance) of x(t), x(t), and hence x̂(t) closely approximates the actual arm/hand position/configuration trajectory. Of note, because d ≃ −250 ms, reconstruction of xl(t − d) based on neuronal firing ril(t) i = 1, …49 amounts to a prediction of the actual trajectory 250 ms in advance. In many tasks, poor reconstruction tended to occur in the period −250 < t ≤ 250 ms. In many tasks, poor reconstruction occurred, especially with the fourth and especially the fifth PC. This was not surprising as our models were least sensitive to these because of their lower signal-to-noise ratios.
It may be appreciated in Figure 14 that, especially for PC1-PC4, the changes in component values for approximately the period −250 ms < t ≤ −100 ms are reasonably monotonic and synchronous. This is followed by a second, largely monotonic and synchronous epoch approximately for −100 ms < t ≤ 0 ms during which most of the change occurs in PC3, PC4, and PC5. This pattern corresponds to a reasonably straight, radially directed reach associated segment, followed by a differently directed grasp-associated segment. These signal segments occur ∼250 ms before the actual movements. The remainder of the signal trajectory corresponds to the postgrasp manipulation. Because the magnitude of PC change in the associated with the grasp signal phase is generally substantially less than that in the first phase, the whole two phase signal can be roughly approximated by a single straight radial trajectory as follows:
where xtarget is the (fixed) visually specified, intended limb position/configuration after grasp based on object size, shape, orientation, and location, and ρ(t) is a scaling fraction that transitions from 0 to 1 over the time course of the command generation. Thus, the neuronal firing pattern is predicted to depend on the distance and direction between initial and target limb position/configurations. When all trajectories in a task begin from a common xinitial, then differences in neuronal firing patterns are influenced specifically by the target position/configuration. As xinitial can be supplied by proprioception, Equation 21 together with Equation 11 or Equation 13, represent simple, explicit models of visuo-proprioceptive integration that appear to be consistent with physiological data. Finally, it may be appreciated from Figure 14 that the net direction to the final target in signal space is revealed very early in the movement command. Only in the rod and ball tasks, are two components (PC4 and PC5) of the net direction significantly mis-specified within the first few milliseconds after comman generation onset.
Discussion
Our investigation demonstrates that MGF models are sufficient to account for much of the firing behavior of PMv neurons during selected reach-grasp-manipulate tasks. Importantly, a single MGF representation is generally predictive of a PMv neuron's firing across different tasks. If we consider a tuning curve to be the neuronal firing response as a function of a monotonic change in location along a path within the input signal space, then the MGF formulation determines inherently that many neuronal tuning curves will be bell-shaped. Combinations of two-dimensional Gaussian functions have been found to describe visual receptive fields (Hubel and Wiesel, 1962; Jones and Palmer, 1987; Duhamel et al., 1997), parietal receptive fields (Brotchie et al., 1995) and hippocampal place cell fields (Brown et al., 1998). On the other hand, heretofore Gaussian-based receptive fields identified experimentally have involved signals with at most 2 degrees of freedom. To our knowledge, this is the first description of high-dimensional Gaussian neuronal input fields.
We found rank-1 and disc-like (nearly rank-1), full-rank MGF to be comparably efficient in relating each neuron's input signal to its firing. Rank-1 MGFs are approximations that are especially adequate for modeling responses to signals that lie close to the neuron's estimated b vector. Because each side of a rank-1 MGF can be approximated by a log-linear field (Fig. 10), it is understandable that log-linear or even linear models can also provide good firing pattern fits for some tasks in some neurons. The apparent value of MGF is comparable or improved fit across tasks in all PMv neurons. It will be important to test MGF models' abilities to explain firing in tasks different from those used to estimate its parameters.
Assuming MGF, we found limb position/configuration to be more explanatory of neuronal firing than velocity or acceleration. This is consistent with the observed specificity of PMv neuronal firing during both movement and maintenance of grasps (Raos et al., 2006). The known input–output connectivity of PMv supports the likelihood of converging somatosensory information from SII (directly to PMv and indirectly via AIP and ventral intraparietal area (VIP)) and visual target information from AIP (Dancause et al., 2006) and VIP (Duhamel et al., 1997). However, our results indicate that trajectory-like information appears to unfold ∼250 ms before movement. Moreover, although SII receives ample sensory input from the limb, it does not appear to be a simple relay station. Rather, it is connected with working memory and has been implicated in transiently holding signals for subsequent comparison and decision making (Romo et al., 2002; Pleger et al., 2003). Premovement activity in PMv is also consistent with observations in many so-called “object-type” visual-motor neurons in AIP of firing increase on the order of 200–250 ms before movement onset (Murata et al., 2000). Together, our findings are consistent with the possibility that AIP and PMv are involved in generating a multidimensional feedforward reference trajectory for the limb position/configuration. Because just before movement the level of activity in PMv “motor” neurons is routinely quite low compared with the peak levels achieved during movement and holding (Raos et al., 2006), the command signal is likely to be incrementally feedforward with respect to some recently held limb position/configuration, presumably supplied from SII. This interpretation is consistent with observed shifts in PMv tuning curves with changes in hand position (Graziano et al., 1997) and our failure to find a continuous position-related signal that lagged body movement. It also supports and extends concepts presented previously (Fagg and Arbib, 1998).
Although the signal encoded in PMv may be much more like an intended position trajectory than either a pure velocity or acceleration waveform, its precise morphology is not established. In particular, a contribution of a velocity signal at comparatively low amplitude or within a subset of PMv neurons is not excluded. Filtered versions of the input signal that would account for dynamic processes, including, for example, rhythms and transients in the surrounding cortical network, might provide an even better explanation of the neural firing data, given that many of the higher-frequency features are not captured by our models (Fig. 10). Similarly, the incorporation of additional signal components related to cerebellar and/or force control circuits could be beneficial. As dynamic forces would be expected to correlate fairly strongly with angular velocity and/or acceleration, it appears that either PMv does not strongly encode such specifications or dynamic force command signals are comparatively small in this task. These considerations underscore the experimental design limitation that position, velocity, acceleration, and force were not systematically dissociated in the tasks. Figure 10 may then be best viewed as modeling the component of the encoded signal that was linearly related to position. Although such signals accounted for much of the firing variance here, more studies will be needed to understand the role of PMv conclusively.
Our analysis showed that strong modulation of any individual neuron's firing occurs primarily along a single direction within the signal space. In AIP and PMv, visuomotor neurons have been found to “prefer” certain final hand configurations (i.e., grip type) (Rizzolatti, 1998) in the sense that they fire most vigorously when certain target objects are viewed and/or engaged (Murata et al., 2000; Raos et al., 2006). Appropriately, Figure 6 illustrates that a neuron may be expected to fire more intensely for trajectories toward some intended final limb position/configurations. Figure 6A predicts that the strength of a neuron's preference in terms of differential firing intensity with respect to direction should increase as the movement command unfolds (toward a radius of μ). This is consistent with experimental findings (Raos et al., 2006). Moreover, in light of Figure 14, it may be appreciated that the general direction of the initial limb position/configuration may become evident very early in the movement command (i.e., when ρ(t) is small in Eq. 21). We consider that. in studies incorporating a passive observation period before movement (e.g., Murata et al., 2000; Raos et al., 2006), the AIP-PMv system may first become primed for movement by a moderate level of activity, corresponding to a substantially constant nonzero ρ(t), being maintained in visual-motor neurons. If so, then it would be expected that some neuronal selectivity would be revealed as early as target presentation (Raos et al., 2006). Also, if during observation the primate were presented with the same target at different locations in the workspace, a bell-shaped tuning curve should be elicited because of changing values for xtarget with fixed xinitial (Eqs. 11, 13, and 21). And if repeated with a changed initial position of the limb, then a corresponding shift would be predicted in this tuning curve (Graziano et al., 1997). Figure 6A also implies a somewhat log-linear responsiveness whenever signal trajectories are much shorter than μ. Figures 6B and 13 verify moreover, that for certain trajectory-MGF relationships, a substantially linear (affine) encoding results. Because MGF can approximate both linear and log-linear models over subregions of input signal space, it is likely that MGF can provide comparable, if not improved, modeling of firing in motor cortex and perhaps other cortical areas.
We noted that the b vector directions appear to be distributed widely. Also, the vast majority of the MGFs were “local” and many with μ and σ small or modest in size respect to the input signal, which implies strongly nonlinear modulation over the course of movement (Fig. 10). Moreover, the range of μ values is large with respect to most values of σ and the range of μ enveloped all observed s(t). Together, these indicate both that Figure 6 is a qualitatively reasonable and that any particular value of s(t) will cause near-maximal activation of only a distinct and limited subset of neurons. Given that PMv projects directly to motor cortex, it appears that changing small sets of strongly active neurons in PMv might generate corresponding time-varying shifts in muscle activations. That is, MGF may facilitate implementation of a high-dimensional lookup table. Neurons with fields remote to the input subspace of PC4 and PC5 would be expected to be concerned preferentially with proximal joints, as many in area F4, whereas those with fields remote from PC1, PC2, and PC3 preferentially with distal movements, as in area F5. However, it also appears that a nontrivial number of our recorded neurons is best suited to modulate strongly with coordinated proximal and distal actions.
We were able to reconstruct high-dimensional trajectories with reasonable fidelity 250 ms in advance using filtered neuronal firing data from a modest-sized dataset. The quality was generally slightly better for M2 than M1 as more data were available for M2. We used multitrial data and bidirectional smoothing. Therefore, this exercise does not indicate directly the performance achievable in real-time, causal decoding. However, it demonstrates the invertibility of a collection of MGF to yield an encoded signal. The potential added value of a more accurate nonlinear encoding model is reflected in the uniformly superior reconstruction performance of rank-1 MGF and near rank-1 MGF models in both animals.
Footnotes
This work was supported by National Science Foundation Grants EFRI-1137237 and ECC-1055560 to S.V.S. We thank the staff of Dr. Marc H. Schieber's laboratory at University of Rochester, Rochester, NY, for doing valuable discussions and providing us data from both monkeys.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. Rahul Agarwal, Department of Biomedical Engineering, Johns Hopkins University, 3400 N Charles Street, Baltimore, MD 21218. rahul.jhu{at}gmail.com
