## Abstract

Microstimulation is widely used in neurophysiology to characterize brain areas with behavior and in clinical therapeutics to treat neurological disorder. Current intensity and frequency, which respectively influence activation patterns in spatial and temporal domains, are typically selected to elicit a desired response, but their effective influence on behavior has not been thoroughly examined. We delivered microstimulation to the primate superior colliculus while systematically varying each parameter to capture effects of a large range of parameter space. We found that frequency was more effective in driving output properties, whereas properties changed gradually with intensity. Interestingly, when different parameter combinations were matched for total charge, effects on behavioral properties became seemingly equivalent. This study provides a first level resource for choosing desired parameter ranges to effectively manipulate behavior. It also provides insights into interchangeability of parameters, which can assist clinical microstimulation that looks to appropriately control behavior within designated constraints, such as power consumption.

- microstimulation
- deep brain stimulation
- superior colliculus
- motor control
- spatiotemporal

microstimulation in different brain structures has become a widely used technique to report causality between stimulation-induced modulation and resulting behavior. Microstimulation experimentation encompasses studies on motor control, cognition, neuroprostheses, and more (for reviews, see Cohen and Newsome 2004; Gandhi and Katnani 2011; Histed et al. 2012; Rebesco and Miller 2011). Such causal links have facilitated the transition of stimulation into the clinical setting as a therapeutic treatment to restore or suppress motor and/or psychiatric function in patients with different types of neurodegenerative diseases (Holtzheimer and Mayberg 2011; Perlmutter and Mink 2006). As a result of microstimulation's ability to generate responses similar to natural physiologic behavior, the methodology has become an accepted means to gauge brain function. The two stimulation parameters most often varied are current intensity and frequency. For example, skeletomotor studies have used ranges of current intensity and frequency to understand motor cortex signals (Graziano et al. 2002). On the other hand, deep brain stimulation varies both parameters to understand parkinsonian symptoms (Moro et al. 2002). Despite the ability of microstimulation to manipulate behavior, the efficacy of stimulation parameter settings has not been thoroughly evaluated.

We asked whether the relative impact of current intensity vs. frequency on behavioral output could be assessed by testing the sensitivity of the superior colliculus (SC) with different patterns of microstimulation. The SC plays a critical role in the generation of saccades. Movements are preceded by a high-frequency burst across an ensemble of neurons, where the locus of the response dictates the desired vector and the vigor of the burst contributes to the speed of the movement (Sparks and Mays 1990). The subcortical brain structure serves as a promising candidate site in which to test different patterns of current intensity and frequency since stimulation of the SC evokes saccades whose metrics and kinematics are related analogously to the site and parameters (Groh 2011; Stanford et al. 1996; Van Opstal et al. 1990). Furthermore, the laminar structure, topographical map, and direct role in motor control of the SC provide a highly constrained model system, in which the impact and interactions of stimulation, to yield output, can be understood with limited complexity. With well-defined spatial and temporal population responses that reflect sensitivity to differential changes in stimulation, the impact of each parameter may be discriminated well by the oculomotor system.

We systematically varied both current intensity and frequency within a large range, from subthreshold to suprathreshold levels, and correlated the different input combinations to resulting behavioral properties to obtain well-defined response surfaces. We examined saccade features under two conditions: *1*) one stimulation parameter was fixed while the other varied, or *2*) both parameters varied but the total charge input was held constant. Our results demonstrated that frequency was more influential in inducing optimal behavioral properties, whereas current intensity was less effective and produced more gradual changes. Nevertheless, we found that by preserving total charge, different combinations of current intensity and frequency generated similar outputs, thus providing a sense of interchangeability. Overall, the evaluation on the relative impact of different regions of parameter space on evoked behavioral properties provides a useful measure for scientific and clinical work, where stimulation parameters are customized to achieve desired behavioral responses within specified boundaries.

## METHODS

All procedures were approved by the Institutional Animal Care and Use Committee at the University of Pittsburgh and complied with the guidelines of the Public Health Service policy on Humane Care and Use of Laboratory Animals.

### Subjects and Surgical Procedures

Two juvenile, male rhesus monkeys (*Macaca mulatta*) underwent one or more surgeries in a sterile environment and under isoflurane anesthesia. The initial procedure consisted of placing a Teflon-coated stainless steel wire (Baird Industries, Hohokus, NJ) under the conjunctiva of one eye and securing a head-restraint post to the skull. In the second procedure, one cylinder was cemented over a craniotomy. The chamber was placed stereotactically on the skull, slanted posteriorly at an angle of 38° in the sagittal plane. This approach allowed access to both colliculi and permitted electrode penetrations normal to the SC surface. After each surgery, the monkey was returned to its home cage and allowed to fully recover from surgery. Postoperatively, antibiotics and analgesics were administered as indicated in the protocol.

### Experimental Procedures and Behavioral Tasks

Visual stimuli, behavioral control, and data acquisition were controlled by a custom-built program that uses LabVIEW architecture on a real-time operating system supported by National Instruments (Austin, TX) (Bryant and Gandhi 2005). Each animal was trained to sit in a primate chair with its head restrained and a sipper tube was placed near the mouth for reward delivery. The animal sat inside a dome surrounded by two alternating magnetic fields that induced voltages in the eye coil and thus permitted measurement of horizontal and vertical eye positions (Robinson 1963). The animal fixated targets on the isoluminant wall of the dome that were rear reflected from a circular mirror. Antiwarping software obtained from Paul Bourke, University of Western Australia, allowed reflections from the mirror to appear undistorted and for distances to be properly transferred onto a curved surface. The monkey sat in the center of the dome, which had a radius of 1 m and spans ±150° horizontally and ±30° vertically of the visual field. A photodetector, positioned outside the animal's field of view, detected the actual time of appearance of visual objects, which was then used to correct for time shifts induced by the projector's refresh rate.

Both animals were trained to perform the oculomotor gap task. Every trial began with directing the line of sight to a fixation point for 300–500 ms before it was extinguished. The fixation point was always at the straight-ahead location and was always kept constant across all trials for each stimulation site. Following a 200- to 400-ms “gap” interval, during which the animal was required to maintain the same eye position, another stimulus was illuminated in a random location in the visual periphery. Incorporation of the gap interval permitted fixation to become disengaged from a visual stimulus prior to saccade preparation, allowing the oculomotor system to be more responsive to incoming visual and/or stimulation input (Sparks and Mays 1983). Each animal was permitted 500 ms to redirect its visual axis on the saccade target and hold gaze steady for 300–500 ms to earn a liquid reward. As the animal performed this task, a platinum iridium microelectrode (1.0–1.5 MΩ; MicroProbes for Life Science, Gaithersburg, MD) was advanced with a hydraulic microdrive (Narashigie, Tokyo, Japan). The superficial layer of the SC was first identified by the presence of distinctive bursting of background activity associated with flashes of room lights. The electrode was then driven deeper into the SC until saccadic motor bursts were identified. At this stage, stimulation (40 μA, 400 Hz) was delivered during the gap interval, 100 ms after fixation offset, to determine the vector coordinates. The depth of the electrode was then minimally adjusted to obtain the shortest possible latency of the stimulation evoked saccade (20—40 ms). Train duration was manually set (range: 100—300 ms) and always long enough to allow for completion of the stimulation evoked movements. Moreover, the peripheral saccade target was illuminated only after stimulation offset. Stimulation was delivered through the electrode on 20% of trials.

### Microstimulation Permutation

The objective of the experiment was to assess the effects of a large range of stimulation parameters on the elicited behavioral response. Thus, both current intensity and frequency were incremented within a substantial range at each stimulation site. The minimum of the range was set to either evoke a very small movement or no movement (never lower than 10 μA and 100 Hz), and the maximum was always set in the suprathreshold range (never higher than 70 μA and 700 Hz). Differences in the chosen parameter ranges for each site could be attributed to factors such as depth of electrode and/or electrode tip distance from neurons and fibers (Ranck 1975; Tehovnik 1996). The sampling increment for each range was typically 5 μA and 50 Hz or 10 μA and 100 Hz, but was always chosen such that data sets had a total of 25, 36, or 49 possible combinations of current intensity and frequency (∼5 trials of each combination was collected). Example parameter ranges consisted of current intensity (20, 30, …, 60, 70 μA) and frequency (200, 300, …, 600, 700 Hz) or current intensity (10, 15, …, 45, 50 μA) and frequency (100, 150, …, 450, 500 Hz). Note that the parameters were intentionally set to be tenfold apart from one another. The intuition behind this design was to collect two types of conditions: one in which only one parameter varied relative to the other and the other in which both parameters varied but the total charge input was held constant (i.e., 10 μA × 200 Hz = 20 μA × 100 Hz). Total charge, computed from different combinations of the parameters, allowed for a better definition of the relative impact of current intensity vs. frequency on behavioral output.

### Electrical Stimulation

Constant current stimulation trains were generated using a Grass S88X stimulator in combination with Grass PSIU6 isolation units. Trains consisted of cathodic phase leading, symmetric biphasic pulses (0.25 ms). Frequency and current intensity were dictated by each experimental design, typically spanning in the range of 10—70 μA and 100—700 Hz. In all cases, stimulation duration was always long enough to ensure that it outlasted the eye movement.

### Data Analyses

Each trial was digitized and stored for offline analysis. We used a combination of in-house software and Matlab 7.10.0 (R2011a). Horizontal and vertical eye positions along with onset and offset times of the stimulation train were stored with a resolution of 1 ms. Component velocities were obtained by differentiating the eye position signal. Onset and offset of stimulation-evoked saccades were then detected using a standard 30°/s velocity criteria, respectively.

#### Response surface methodology.

Response surface methodology uses multivariate regression to model and analyze the relationships between a response and the variables that influence that response. Evoked movement properties (amplitude, peak velocity, and latency) were defined as a function of current intensity and frequency and quantified with a regressed quadratic model
*Z* is either saccade amplitude or peak velocity; *x* is the current intensity predictor; *y* is the frequency predictor; *a* and *d* denote the coefficients describing the quadratic and linear effects attributed to current intensity, respectively; *b* and *e* denote the coefficients describing the quadratic and linear effects attributed to frequency, respectively; *c* denotes the coefficient for the interaction effect; and *f* denotes the intercept coefficient.

We found that the response surface for saccade latency (stimulation onset to movement onset) was better approximated by an inverse multilinear model. Therefore, we took the inverse of latency values to approximate the data with a simple multilinear model
*Z* is the inverse of saccade latency. All other terms are as described earlier with the quadratic terms set to zero (*a = b* = 0).

The optimal number of parameters for multivariate regression was determined by using the *F*-test, the coefficients were determined using a nonlinear least-squares optimization routine (Levenberg—Marquardt algorithm), the significance of each coefficient was determined by a *t*-score, and the coefficient of determination (*R*^{2}) was used to assess goodness of fit.

#### Canonical analysis of quadratic model.

Canonical analysis was a technique used to examine fitted quadratic models by assessing second-order equations (*Eq. 1*) in canonical or standard form
*Z*_{s} is a stationary point that describes a region of parameter space in which little change occurs in response to changing the controllable factors; λ_{1} and λ_{2} are eigenvalues that describe the shape of the modeled surface; and *E*_{1} and *E*_{2} are eigenvectors or principal axes that convey the direction of change in yield for the modeled surface. *Equation 3* may be interpreted as follows. Instead of measuring the effects of each factor from the origin and in the direction of the original axes (*x* and *y*; *Eq. 1*), the effects are measured from the stationary point and in the direction of the principal axes (*E*_{1} and *E*_{2}; *Eq. 3*). Therefore, a change in yield when moving away from the stationary point is defined by the expression on the right-hand side of *Eq. 3*. The magnitude of the coefficients, λ_{1} and λ_{2}, defines the amount of change in yield in the direction of *E*_{1} and *E*_{2}, respectively. Moreover, the sign of the eigenvalues describes the nature of the stationary point, in which all negative defines a maximum, all positive defines a minimum, and mixed defines a saddle point. For example, if both λ_{1} and λ_{2} are negative, the stationary point is a maximum. This can be better interpreted by rewriting *Eq. 3* as *Z* = *Z*_{s} − λ_{1} *E* _{1}^{2} − λ_{2} *E* _{2}^{2}, which shows that any point chosen in the new coordinate system (*E*_{1}, *E*_{2}) will result in a loss of yield from *Z*_{s}. Therefore, the stationary point is a maximum with λ_{1} and λ_{2} measuring the fall off in yield in the direction of *E*_{1} and *E*_{2}, respectively.

To derive the components of *Eq. 3*, all second-order equations (*Eq. 1*) were rewritten using matrix notation
*b*_{1} denotes a vector of the first-order parameter estimates; *w* denotes a vector of the controllable factors; and *B* denotes a matrix of the second-order parameter estimates (note the off-diagonal elements of *B* equal half the interaction coefficient). Equating the partial derivatives (d*Z*/d*x* and d*Z*/d*y*) of *Eq. 4* to zero and solving the resulting system of equations, provided values for *w* used to reevaluate the equation and solve for a stationary point. The eigenvalues and corresponding eigenvectors were determined by eigen-decomposition of matrix *B*. For a full account of the usage and mathematics behind canonical analysis on surface responses see Box (1954) and Box et al. (1953).

#### Calculating total charge.

The total charge (TC) elicited by different combinations of current intensity and frequency was calculated using the following equation
*T* is the set period of pulses (1/Hz) in the stimulation train. In biphasic pulsing the first phase, also known as the stimulation phase, has been shown to elicit the physiologic effect and primarily influences activation, whereas the second phase, also known as the reversal phase, is used to reverse electrochemical reactions to minimize extracellular tissue damage, and is believed to have minimal effect on activation (for review, see Merrill et al. 2005). As a result, only the phase leading pulse was used for TC calculation.

The allotted duration used to calculate TC was chosen as a period of time when stimulation was significantly contributing to the saccade feature being analyzed. When analyzing radial amplitude and peak velocity, duration was chosen to be 20 ms before saccade onset until saccade offset (Goossens and Van Opstal 2006; Munoz and Wurtz 1995; Ottes et al. 1986). When analyzing movement latency, duration was chosen to be the time from stimulation onset to saccade onset. We note that additional to the TC equation above, when calculating TC for latency the equation was further divided by duration to compute the TC relative to the rate of incoming pulses. The calculation stems from the notion that the rate of incoming charge is significant to the initiation of a movement (Carpenter et al. 2009; Hanes and Schall 1996).

## RESULTS

We report on a total of 30 microstimulation sites from two monkeys (monkeys 1 and 2: 13 and 17 sites, respectively), sampling a range of the SC motor map that spanned approximately 3° to 40° in amplitude and approximately −75° to 50° in direction. Twenty-seven of the 30 sites were 6 × 6 permutation sets, two were 5 × 5 sets, and one was a 7 × 7 set.

### Defining Response Surface Characteristics

Figure 1*A* illustrates a response surface that describes stimulation-evoked movement amplitude as a function of current intensity and frequency. Qualitative assessment was obtained by generating a surface that connected the mean behavioral response for each current intensity and frequency combination. The surface highlights a concave structure in which amplitude increases monotonically from a minimum to an apex. The individual slices that constitute the surface plot are shown in Fig. 1, *B* and *C*, where qualitative differences can be observed when varying frequency vs. current intensity. Variation in frequency seems to generate nonlinear effects as opposed to the linear-like changes seen when varying current intensity. A similar relationship is observed with peak velocity (not shown). Latency, in contrast, is high and more variable at low stimulation parameters and asymptotes to a minimum with suprathreshold stimulation (not shown). To implement conformity in the structure of multivariate regressions, we correlated the inverse of latency values to current intensity and frequency, seen in Fig. 2*A*. In this representation, the inverse of latency also increased monotonically, although the overall shape of the data was much more linear. This can be better appreciated in the individual slices of the surface shown in Fig. 2, *B* and *C*. A qualitative assessment suggests that variation in frequency seems to cause much steeper slopes than variation in current intensity.

Response surface methodology was used to quantify the surface features across all stimulation sites (see methods). Each data set was first normalized in all three dimensions (*x*-, *y*-, and *z*-axis were each divided by their respective maximum) to make comparisons across predictors and data sets, and then fit with a multivariate model. Figure 3, *A* and *B*, illustrates fits for radial amplitude and peak velocity from the same stimulation site used in Fig. 1*A*. As above, the concave structure and monotonic increase in output variables as a function of stimulation parameters are preserved. A quadratic model provided the proper approximation of the surface characteristics (*Eq. 1*). The *F*-statistic revealed that 48 of the 60 surface fits were significantly (*P* < 0.05) improved when fitted with a quadratic rather than multilinear model. Figure 3*C* illustrates a multilinear model fit to the surface characteristics of the inverse of latency (*Eq. 3*). The *F*-statistic revealed that only 1 of the 30 surface fits was significantly (*P* < 0.05) improved by using a quadratic model. Figure 3*D* illustrates a bar graph of the coefficients of determination (*R*^{2}) derived from all three behavioral output surface fits. Note, *R*^{2} values are always >0.75, with the majority >0.85, demonstrating that the regressed models are significantly capturing the characteristics within the behavioral response data.

To summarize the response surface shape across all stimulation sites and assess how each stimulation parameter contributed to the shape, we examined the coefficients derived from the multivariate regressions. Figure 4*A* illustrates the resulting distributions for each coefficient across all response surfaces. Note that all quadratic terms were found to be significant (*t*-score, *P* < 0.05) and negative, signifying a common concave shape for all amplitude and peak velocity surfaces. All linear terms were found to be significant (*t*-score, *P* < 0.05) and positive, signifying a common monotonically increasing direction across all responses surfaces. To better visualize the relative magnitude of the derived effects, we compare the coefficient associated with the current intensity and frequency terms. Figure 4, *B* and *C*, demonstrates that almost all quadratic and linear coefficients derived from radial amplitude (black) and peak velocity (light gray) lie above the line of unity. The same trend exists for linear coefficients derived from inverse of latency (Fig. 4*D*). Overall, the comparisons signify that the frequency predictor is the more prominent factor in contributing to the shape of a behavioral response surface.

The interaction term for each model was also assessed. For amplitude and peak velocity analyses, the interaction coefficient in the quadratic model was statistically significant (*t*-score, *P* < 0.05); however, Fig. 4*A* illustrates that values derived for the coefficient were close to zero. Consistent with this result, there was no significant (*F*-test, *P* < 0.05) difference when fitting response surfaces for amplitude or peak velocity with a pure quadratic (no interaction term) vs. full quadratic model. In contrast, the interaction effect derived from inverse of latency surfaces was both significant (*t*-score, *P* < 0.05) and large in magnitude (Fig. 4*A*). In addition, all interaction coefficients were positive, implying a synergistic effect of the two predictors on latency.

### Relative Impact of Stimulation Parameter

To evaluate how each stimulation parameter specifically affected behavioral output properties, the total charge (TC) of each parameter combination in a data set was calculated (see methods). A noteworthy piece of information for performing such a calculation is that it allows a three-dimensional (3D) response surface to be collapsed into two dimensions by scaling current intensity and frequency into common units (μA/ms), allowing for a better definition of the relative effects of each parameter. TC values were parsed into two groups: current series, calculated from combinations that have the same frequency but different current intensities, and frequency series, calculated from combinations that have the same current intensity but different frequencies. As a result, each series contained slices of the 3D response surface that encompassed how variation in one parameter affected behavioral output while the other remained fixed. Slices derived from radial amplitude and peak velocity surfaces were regressed with a one-phase exponential (1 − e^{−xβ}; β: slope coefficient), whereas slices from the inverse of latency surfaces were fit with a linear equation (β_{1}*x* + β_{2}; β_{1} and β_{2}: slope and intercept coefficient). Any regression resulting in a coefficient of determination <0.6 was removed from analysis (<10% of data removed).

Figure 5*A* illustrates a histogram of the exponential slope coefficients derived from current series (blue) and frequency series (red) taken from radial amplitude surfaces. The distribution is a result of regression across all stimulation sites. Note that the time constant for a one-phase exponential function is inversely proportional to the regressed slope; therefore, the results demonstrate that variation in frequency generates significantly (Wilcoxon, *P* < 0.001) larger slopes or shorter rise times in behavior compared with variation in current intensity. To visualize the impact of having different time constants, Fig. 5*B* illustrates slices from a frequency series (*left*) and current series (*right*) taken from one stimulation site. The color of a slice indicates the different level at which the noniterating parameter was fixed, whereas each slice is a result of the iterating parameter. We observed that within the frequency series the iteration of frequency causes concave nonlinearity by quickly driving the behavioral output to a plateau (thus, the shorter rise times). In contrast, within the current series the iteration of current intensity causes more gradual linear-like effects. The same analyses performed on peak velocity surfaces produced the same comparable result (not shown).

Figure 6*A* illustrates a histogram of the linear slope coefficient derived from current series (blue) and frequency series (red) taken from inverse of latency surfaces. The results demonstrate that variation in frequency generates significantly (Wilcoxon, *P* < 0.001) larger slopes than variation in current intensity. Figure 6*B* illustrates slices from a frequency series (*left*) and current series (*right*) taken from one stimulation site. The slopes of slices in the frequency series are noticeably steeper, indicating frequency as the more prominent factor in dictating the onset of movements.

### Canonical Analysis

The characteristics captured in fitted quadratic models (amplitude and peak velocity) were further interpreted by assessing their canonical form, where the effects of varying current intensity and frequency were measured from a stationary point and in the direction of the derived principal axes of each modeled surface (see methods). We demonstrate an example of this analysis in Fig. 7*A*, which plots the contours of the example response surface in Fig. 1*A*. The color bar of Fig. 7*A* conveys normalized saccade amplitude, indicating that the maximum of the surface is in the top right corner of the figure and the minimum in the bottom left. Superimposed atop the contours are the derived principal axes. Traditionally, the axes are centered at the stationary point, but for a better visualization of the results we arbitrarily shifted them to the center of the figure. The relative magnitude of the eigenvalues is illustrated by the thick blue portions on each axis. The eigenvalues were both found to be negative (λ_{1} = −1.16, λ_{2} = −0.51), as were the eigenvalues obtained across all stimulation sites (λ_{1}, amplitude: mean = −1.64, median = −1.65, SD = 0.44; λ_{2}: mean = −0.68, median = −0.64, SD = 0.40; λ_{1}, peak velocity: mean = −1.50, median = −1.46, SD = 0.48; λ_{2}: mean = −0.59, median = −0.61, SD = 0.34). As interpreted from *Eq. 3* (see methods), negative eigenvalues define the stationary point as a maximum. Therefore, any direction taken from the stationary point will result in a loss of yield, with the fall off being greatest in the direction of the principal axis corresponding to the largest eigenvalue. Since λ_{1} is >λ_{2} in the example, we find that moving across the surface (*right* to *left*) in the direction of the first principal axis (*E*_{1}) results in a rapid loss in yield (as seen by the change of color in contour ridges); in contrast, moving across the surface in the direction of the second principal axis (*E*_{2}) results in considerably less change in yield. This result was consistent across all surfaces as the magnitude of the first eigenvalue was always found to be greater than the second (λ_{1} > λ_{2}). The slope (ratio of frequency divided by current intensity components) of the principal axes *E*_{1} and *E*_{2} was also assessed to gauge how much each stimulation parameter contributed to the directions of the axes. For the first principal axis, *E*_{1}, in Fig. 7*A*, the ratio (0.57/0.82) indicates that current intensity contributes most to the direction, signifying that changes in the level of current intensity are associated with greater changes in yield. Accordingly, frequency contributed more to the direction of the orthogonal axis, *E*_{2} (slope = 0.82/0.57), signifying that changes in the level of frequency are associated with smaller changes in yield. These results were consistent because current intensity was always found to contribute most to the direction of the first principal axis (across all sites, ratio < 1; amplitude: mean = 0.47, median = 0.45, SD = 0.29; peak velocity: mean = 0.61, median = 0.71, SD = 0.29), and frequency to the second principal axis (across all sites, ratio > 1; amplitude: mean = 4.53, median = 2.25, SD = 5.32; peak velocity: mean = 2.47, median = 1.39, SD = 2.56).

In addition, we arbitrarily divided each contour response into a low and high region and reevaluated how each parameter dictated the direction of the principal axes in each subregion. Figure 7*B* illustrates the principal axes corresponding to the largest eigenvalue derived from the low and high region of the example surface. The reevaluation demonstrates that in lower parameter space the ratio of frequency divided by current intensity is far more equivalent (slope = 0.71/0.70), indicating similar effects on yield. It is not until higher parameter space that the ratio shifts in favor of current intensity (slope = 0.51/0.86), thus demonstrating the same result as the canonical analysis above. Current intensity contributing more to the change in yield at higher levels agrees well with the longer time constants derived from regression of current series slices vs. frequency series (see *Relative Impact of Stimulation Parameter*).

### Interchangeability of Stimulation Parameters

Having established the different effects of current intensity and frequency on behavioral output properties, we next assessed whether the stimulation parameters had compensatory effects. Observing contour responses like that in Fig. 7, we find that a wide range of stimulation parameters correspond to each isoamplitude ridge. Therefore, as a first level approximation we can appreciate that the parameters compensate for one another in the sense that departure from one ridge, by changing one stimulation parameter, can be compensated for by suitably changing the other parameter (returning to that ridge). Figure 8, *A* and *B*, illustrates a qualitative assessment of the contours from another example radial amplitude surface and an inverse of latency surface. The arrows in each figure demonstrate examples of how each factor must change to return to a previously departed ridge. A qualitative assessment of this manner can be made across the entire contour response, where the direction of the ridges indicates the amount of change that will be needed for one parameter to compensate for the other. Accordingly, we find that the direction of ridges for radial amplitude is sensitive to the region of parameter space (see *Canonical Analysis* above); therefore, the amount of change for one parameter to compensate for the other varies. In contrast, the direction of ridges for the inverse of latency is consistently dictated by frequency; therefore, to compensate for small changes in frequency, there must always be large changes in current intensity.

Along with qualitatively assessing the contours of surfaces for compensatory sensitivity, we also wanted to quantify whether current intensity and frequency were interchangeable factors. Interchangeability was gauged by matching total charge (TC) values within each data set that were not significantly (Wilcoxon, *P* > 0.05) different, and then correlating the generated behavioral outputs from each matched pair. Figure 9, *A* and *B*, illustrates the correspondence of paired radial amplitude and paired peak velocity outputs. Within each pair, the output evoked by TC calculated from parameter combinations where frequency (divided by 10) was greater than current intensity was always plotted on the *x*-axis; combinations in which current intensity was greater than frequency (divided by 10) was always plotted on the *y*-axis. Comparing paired outputs in this manner provided a better evaluation of the relative impact of each parameter. Correlation was assessed by linear regression in which the slope and correlation coefficient for each output distribution was: slope = 1.06, *R*^{2} = 0.965 and slope = 0.941, *R*^{2} = 0.899, respectively. To quantify the amount of shift from perfect correspondence we calculated the residual error (Euclidean distance from orthogonal projection) of each data point from the line of unity. The inset of each figure illustrates a histogram of the magnitude of error (radial amplitude: median = 0.001, SD = 0.039; peak velocity: median = −0.006, SD = 0.043). Each distribution was not significantly (Wilcoxon, *P* > 0.05) different from a distribution centered at zero. Overall, the insignificant shift from unity and high degree of correlation demonstrates that the relative impact of current intensity and frequency are comparable when TC is preserved. Despite the different relative impact observed when one parameter is fixed and the other varied, the result conveys a strong sense of interchangeability, in which the parameters can compensate for different relative effects as long as the TC remains constant. In contrast, current intensity and frequency did not show the same interchangeability for paired latency outputs. Figure 9*C* illustrates the correspondence of paired inverse of latency outputs. Linear regression of the distribution derived a slope of 0.804 and correlation coefficient of 0.899, whereas the distribution of residual error (median = −0.06, SD = 0.065) was found to be significantly (Wilcoxon, *P* < 0.05) different from a distribution centered at zero. The relatively shallow slope and significant shift below the line of unity demonstrates that in preserving TC the relative impact of frequency, on driving the onset of behavior, is greater than current intensity.

## DISCUSSION

We examined the relative impact of current intensity and frequency on the behavioral properties of saccades evoked by microstimulation of the SC. Our results demonstrate that activity induced for motor control in the saccadic system reflects a different sensitivity to each parameter. We found that the rate of incoming pulses was the more prominent factor in generating optimal responses. Varying frequency from low to high levels caused rapid saturation in saccade amplitude and peak velocity, and steep reduction in reaction time. In contrast, saccade properties displayed more gradual, linear-like trends as current intensity increased. Interestingly, despite the different effects of each parameter, high correlations were conveyed when amplitude and peak velocity properties were paired by matching total charge values. The result demonstrates that different combinations of current intensity and frequency can generate similar output properties provided the total charge, inducing the activity, is equivalent. No such result, however, was found when pairing latency values because frequency was always the strongest factor in dictating the onset of movements.

### Neural Mechanisms of Action

Prior to generating visually guided saccades, low-level discharge in an ensemble of collicular neurons accumulates gradually toward an activation threshold, at which point it transitions abruptly into a high-frequency burst to initiate a saccade. The population of active neurons can be envisioned as a Gaussian mound, whereby the metrics and kinematics encoded in the activity are defined by the spatiotemporal features (Sparks et al. 1976). It is generally assumed that SC microstimulation also activates a Gaussian distribution as a result of (suprathreshold) stimulation-evoked saccades being nearly equivalent to visually guided movements (Van Opstal et al. 1990). However, previous studies have demonstrated that suboptimal microstimulation can attenuate both kinematics and metrics (Groh 2011; Katnani et al. 2012; Stanford et al. 1996; Van Opstal et al. 1990); moreover, our current study provides a systematic investigation of the differential effects of current intensity and frequency on movement-evoked properties. In this section, we explain our results by speculating on the spatiotemporal features of the overall activity produced during microstimulation. With frequency and current intensity defining the rate and magnitude of incoming pulses, respectively, and being the major factors that constitute the total extracellular current flux to which neural elements respond, we will discuss the collicular networks responsive nature to each factor, reflected through the evoked behavioral trends.

Typically, a very small number of neural elements are directly activated with each pulse delivered through the microelectrode (Histed et al. 2009), and a substantially larger fraction of the population response is recruited through synaptic activation (McIlwain 1982). The frequency of microstimulation contributes to the ensemble response predominantly by driving the temporal integration of action potentials, whereas current intensity regulates mainly spatial summation (Ranck 1975; Tehovnik 1996; Vokoun et al. 2010). During the epoch of time when microstimulation induces activity to trigger a movement, our data demonstrate (Fig. 6) that saccade initiation is determined primarily by frequency. Consistent with the framework that activity must meet a threshold to trigger a movement (Hanes and Schall 1996), temporal summation could be the key integrative mechanism that accumulates the population activity toward that threshold. Once threshold is attained, the low-level discharge transitions to a high-frequency burst, which we suspect is driven by biophysical properties and/or intracollicular interactions (Isa and Hall 2009). Due to the nonlinear scaling and quick saturation of amplitude and peak velocity with stimulation frequency, we speculate that frequency provides the temporal precision for the kinetics and interactions, excitation and inhibition of the collicular network, to operate in a manner that allows for the production of the high-frequency burst mode needed during a movement. In contrast, we expect that the collicular network is less responsive to current intensity. One interpretation is that spatial recruitment of collicular neurons scales approximately linearly with current intensity. A smaller distribution of activity could weaken the intracollicular network and therefore reduce the duration for which the high-frequency burst can be sustained. Alternatively, the weaker drive from the SC to the burst generator could prematurely resume omnipause neuron activity, which inhibits the saccadic system. Either or both of these mechanisms could account for the gradual linear-like trends generated across different levels of current intensity.

Further insight on the relative impact of current intensity and frequency is provided when assessing the compensatory effects between parameters. Our results demonstrate that different combinations of current intensity and frequency delivered during the intrasaccadic period can produce similar amplitudes and peak velocities on average. Given that each parameter is thought to influence neural activity differently, this result might not be expected. However, the finding supports the idea that past the point of triggering a movement, spatial and temporal patterns of activity can begin to compensate for one another, where the degree of compensation to produce similar patterns of depolarization is dependent on how responsive the neural system is to each parameter.

Additional microstimulation parameters can be manipulated to influence neural activation and perhaps behavior. Aside from current intensity and frequency, the pulse width, shape of pulse train, and interphase delay of biphasic pulsing, can be systematically varied. These factors were not incorporated in our investigation, because every additional component causes a significant increase in the number of permutations. Large permutation sets are difficult to obtain due to electrode longevity, animal fatigue, and the desire to preserve tissue health. Nevertheless, we speculate from the relationships quantified in this investigation that, similar to our interpretation of varying current intensity and frequency, the responsive nature of the kinetics and interaction connections in the neural network can explain the differential effects of microstimulation parameters on behavior. Increasing or decreasing the width of a pulse or interphase interval, which can effectively trigger a movement, is essentially equivalent to respectively decreasing or increasing the frequency of pulses. Therefore, we expect the manipulation of these parameters to influence temporal activity patterns in a manner that will generate behavioral trends similar to what was observed when varying frequency (Figs. 5 and 6). However, the impact of scaling these trends (equivalent to shifting the red distribution in Figs. 5*A* and 6*A*) may differ depending on the ability of each parameter to engage the neural network with a more natural/physiologic temporal activation pattern (Kimmel and Moore 2007). Modeling efforts have demonstrated that the spatial selectivity (proximity from electrode tip) of neurons (Grill and Mortimer 1995) along with selectivity of cell bodies vs. axons (McIntyre and Grill 2000) can be affected by pulse width. Thus, manipulating the parameter can influence the activation pattern of distributed neural activity in different ways. However, the models generally use a lattice arrangement with an organized layout of fibers and cell bodies, making it difficult to interpret how different activation patterns would truly influence behavior. A systematic study of pulse width or interphase delay on evoked behavior would be an interesting area for future research.

Finally, we show that saccade kinematics evoked from the SC is nonlinearly correlated with stimulation frequency. In contrast, stimulation in the pontine and reticular formations produces movements that continue for the duration of the stimulation (until biomechanical limits). In addition, the kinematics of these movements has been shown to linearly scale with frequency (Cohen and Komatsuzaki 1972; Gandhi et al. 2008; Quessy and Freedman 2004). This discrepancy in stimulation effects suggests a general neural network property, in which direct or indirect recruitment of neurons within a topographically organized structure activates a network level response that can dominate the stimulation-induced activation, whereas the network level response is weaker or nonexistent in regions that operate primarily with a rate code (see note added in proof).

### Significance in the Scientific Setting

Microstimulation is a pillar in neuroscience. The technique has been quintessential in testing hypotheses, revealing emerging properties and defining brain function. A primary motivation behind this study was to provide a robust evaluation of the different factors in stimulation to help solidify the methodology and extend the experimental capability. Understanding the relative impact of varying current intensity vs. frequency can help to appropriately select the parameter space necessary to achieve desired behavior. Data of this nature can be beneficial when designing an experiment to test specific properties of a neural system or to test the predictions of a neural model. For example, within the oculomotor field, the input to output relationship characterized in this investigation provides a resource to challenge existing concepts and models that describe spatiotemporal motor decoding in the superior colliculus (Groh 2001; Katnani et al. 2012). In addition, a well-developed neural field model of the collicular network trained on these microstimulation results could reveal emerging properties in the SC and provide insightful information on the mechanisms underlying saccade generation. Nevertheless, the characterized relationship goes beyond the oculomotor system and provides insight on the general understanding of properties in neural coding. Although effects of stimulation can be specific to the morphology of a neural network (i.e., type of neurons, density of ion channels, etc.), the underlying concepts revealed in this investigation align well with converging evidence in other subcortical and cortical areas that demonstrate an overall higher sensitivity to the temporal component of microstimulation (Butovas and Schwarz 2003; Kimmel and Moore 2007). Furthermore, as a result of tissue impedance and fiber structure being variable across brain region, we propose a unifying dimension of total charge, in which the impact of current intensity and frequency, as a first level approximation, can generalize to other areas of the brain.

### Significance in the Clinical Setting

Therapeutic stimulation in the central nervous system looks to customize parameters to optimize benefit while minimizing side effects and power consumption. Our results provide a useful measure to achieving these goals by demonstrating the effectiveness of current intensity and frequency within different regions of parameter space. For example, trends like those shown in Figs. 7 and 8 provide a first-level assessment of the trade-off between desired yield and power consumption. The data demonstrate that when varying current intensity or frequency behavioral effects occur within two regimes: low parameter space, in which the effectiveness values of current and frequency are comparable, and high parameter space, in which variation in frequency, in contrast to current intensity, generates steady optimal yield (this may be better visualized by observing the saturation effects at high levels of frequency illustrated in the slices of Figs. 1*B* and 5*B*). As a result, we can conclude that frequency has the potential to generate more stable response output. In correspondence to optimize therapeutic stimulation, it may seem suitable then to fix frequency at high levels to stabilize desired clinical effects, and then vary current intensity to minimize observed side effects or meet specified design constraints.

The work from this investigation also draws many parallels to mechanisms of action involved in deep brain stimulation (DBS). Note that a majority of DBS experimentation uses constant voltage stimulation, where varying voltage is analogous to varying current intensity. Previous work in the DBS field has demonstrated that lowering voltage, as opposed to frequency, causes a steep fall off in alleviating symptoms of neurological disorders (Moro et al. 2002). In parallel we have shown that lowering current intensity causes a rapid drop off in output properties (see *Canonical Analysis*); therefore, similar to current intensity, the relative impact of voltage may be explained by the less responsive nature of neural networks to the parameter. Nevertheless, the impact of voltage on symptom alleviation is confounding. Modeling work has demonstrated that even at low levels of voltage, neural elements are still stimulus driven by the set frequency (Hashimoto et al. 2003; McIntyre et al. 2004). Because parkinsonian symptoms are believed to occur due to corrupted temporal firing patterns within relay loops (Guo et al. 2008), one would not expect the impact of voltage to be very critical. Nevertheless, as discussed earlier (see *Neural Mechanisms of Action*) induced temporal activity may be constrained by spatial recruitment, in which frequency may be critical in defining the proper neural command to drive a system but current intensity influences the strength of that command. Therefore, in context to DBS, frequency may be able to effectively drive relay loops in a stable fashion, but the fidelity of that signal is ineffective until the strength is appropriately set, thus alleviating symptoms.

## GRANTS

This work was supported by a National Institute of General Medical Science Institutional Training Grant T32 GM-081760 (to H.A.K.), a National Science Foundation Institutional Training Grant DGE-0549352, National Eye Institute Grants R01 EY-015485 and P30 EY-008098a (to N.J.G.), and National Institute of Deafness and Communication Disorders Grant P30 DC-0025205.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

H.A.K. and N.J.G jointly conceived and designed the experiment. H.A.K performed experiments, wrote code, analyzed data and wrote the paper. N.J.G wrote the paper.

## ACKNOWLEDGMENTS

We thank A. D. Katnani, A. J. van Opstal, and U. Jagadisan for scientific discussions; J. McFerron for software maintenance; and G. Foster for general assistance.

## NOTE ADDED IN PROOF

Importantly, a suggestion of the nature described at the end of *Neural Mechanisms of Action* is not unique to the SC, as recurrent microcircuit architectures have been found in other areas of the brain (**Keller A.** Intrinsic synaptic organization of the motor cortex. *Cereb Cortex* 3: 430–441, 1993). Furthermore, recent work has demonstrated network level responses in cortex (**Logothetis NK, Augath M, Murayama Y, Rauch A, Sultan F, Goense J, Oeltermann A, Merkle H.** The effects of electrical microstimulation on cortical signal propagation. *Nat Neurosci* 13: 1283–1291, 2010). Therefore, the effective influence of microstimulation in the SC may generalize to other cortical and subcortical regions.

- Copyright © 2012 the American Physiological Society