Microstimulation in a spiking neural network model of the midbrain superior colliculus

The midbrain superior colliculus (SC) generates a rapid saccadic eye movement to a sensory stimulus by recruiting a population of cells in its topographically organized motor map. Supra-threshold electrical microstimulation in the SC reveals that the site of stimulation produces a normometric saccade vector with little effect of the stimulation parameters. Moreover, electrically evoked saccades (E-saccades) have kinematic properties that strongly resemble natural, visual-evoked saccades (V-saccades). These findings support models in which the saccade vector is determined by a center-of-gravity computation of activated neurons, while its trajectory and kinematics arise from downstream feedback circuits in the brainstem. Recent single-unit recordings, however, have indicated that the SC population also specifies instantaneous kinematics. These results support an alternative model, in which the desired saccade trajectory, including its kinematics, follows from instantaneous summation of movement effects of all SC spike trains. But how to reconcile this model with microstimulation results? Although it is thought that microstimulation activates a large population of SC neurons, the mechanism through which it arises is unknown. We developed a spiking neural network model of the SC, in which microstimulation directly activates a relatively small set of neurons around the electrode tip, which subsequently sets up a large population response through lateral synaptic interactions. We show that through this mechanism the population drives an E-saccade with near-normal kinematics that are largely independent of the stimulation parameters. Only at very low stimulus intensities the network recruits a population with low firing rates, resulting in abnormally slow saccades.


Introduction
High-resolution foveal vision covers only 2% of the visual field. Thus, the visual system has to gather detailed information about the environment through rapid goal-directed eye movements, called saccades. Saccades reach peak eye velocities well over �1000 deg/s in monkey, and last for only 40-100 ms, depending on their size. The stereotyped relationships between saccade amplitude and duration (described by a straight line) and peak eye velocity (a saturating function) are termed the 'saccade main sequence' [1]. The acceleration phase of saccades has a nearly constant duration for all amplitudes, leading to positively skewed velocity profiles [2]. In addition, the horizontal and vertical velocity profiles of oblique saccades are coupled, such that they are scaled versions of each other (through component stretching), and the resulting saccade trajectories are approximately straight [3]. These kinematic properties all imply that the saccadic system contains a nonlinearity in its control [3][4][5]. More recent theories hold that this nonlinearity reflects an optimization strategy for speed-accuracy trade-off, which copes with the spatial uncertainty in the retinal periphery, and internal noise in the sensorimotor pathways [6][7][8][9].
The neural circuitry responsible for saccade programming and execution extends from the cerebral cortex to the pons in the brainstem. The midbrain superior colliculus (SC) is the final common terminal and a major point of convergence of descending saccade related signals, and it has been hypothesized to specify the vectorial eye-displacement command for downstream oculomotor circuitry [10][11][12]. The SC contains an eye-centered topographic map of visuomotor space, in which the saccade amplitude is mapped logarithmically along its rostralcaudal anatomical axis (u, in mm) and saccade direction maps roughly linearly along the medial-lateral axis (v, in mm; [10]). The afferent map (Eq 1a) and its efferent inverse (Eq 1b) has been described by [13]: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx þ AÞ 2 þ y 2 > > > > > > = > > > > > > ; Afferent mapping ð1aÞ evidence for an additional anisotropy for upward (v > 0) vs. downward (v < 0) directions, which would lead to slightly different inverse mapping relations than Eq 1b (see Discussion). Each saccade is associated with a translation-invariant Gaussian-shaped population within this map, the center of which corresponds to the saccade vector, (x,y), and a width of σ � 0.5 mm [13,15]. It is generally assumed that each recruited neuron, n, in the population encodes a vectorial movement contribution to the saccade vector, which is determined by both its anatomical location within the motor map, (u n ,v n ), and its activity, F n .

Vector averaging vs. linear summation models
Precisely how individual cells contribute to the saccade is still debated in the literature. Two competing models have been proposed for decoding the SC population: weighted averaging of the cell vector contributions ( [16][17][18]; Eq 2a) vs. linear summation ( [3,9,19]; Eq 2b), respectively, which can be formally described as follows: S SUM ðtÞ ¼ N is the number of active neurons in the population, K n < t the number of spikes in the burst of neuron n up to time t, F n its mean (or peak) firing rate, and M n = (x n , y n ) is the saccade vector in the motor map encoded at SC site (u n ,v n ) (Eq 1b). m n = zM n is the small, fixed vectorial contribution of cell n in the direction of M n , for each of its spikes, with z a fixed, small scaling constant that depends on the adopted cell density in the map and the population size, and δ(t − τ k,n ) is the k'th spike of neuron n, fired at time τ k,n .
The vector-averaging scheme of Eq 2a only specifies the amplitude and direction of the saccade vector, and thus puts the motor map of the SC outside the kinematic control loop of its trajectory. It assumes that the nonlinear saccade kinematics are generated by the operation of horizontal and vertical dynamic feedback circuits in the brainstem [16,20,21], or cerebellum [22,23]. Note also that vector averaging is a nonlinear operation because of the division by the total population activity.
In contrast, the linear dynamic ensemble-coding model of Eq 2b encodes the full kinematics of the desired saccade trajectory at the level of the SC motor map through the temporal distribution of spikes by all cells in the population [9,19,24]. As a result, the instantaneous firing rates of all neurons in the population, usually estimated by their instantaneous spike-density functions, f n (t), together encode the desired vectorial saccadic velocity profile: where S n is the number of spikes of cell n, with the spikes occurring at times t k,n . The Gaussian acts as a linear kernel that smooths the discrete spike into a continuous function (e.g., [25]). Although the models of Eqs 2a and 2b cannot both be right, each is supported by different lines of evidence. For example, electrical microstimulation produces fixed-vector (E-)saccades with normal main-sequence kinematics that are insensitive to a large range of stimulation parameters [10,15,26,27]. If one supposes that electrical stimulation directly activates a large population of SC cells, and that the firing rates follow the (typically rectangular) stimulation profile, a vector-averaging scheme with downstream dynamic feedback circuitry readily explains why E-saccades are normal main-sequence, since the center of gravity of the population specifies the desired saccade vector only, regardless the firing rates.
In addition, reversible inactivation of a small part of the SC motor map produces particular deficits in the metrics of visually-evoked (V-)saccades that may not be readily explained by the linear summation model of Eq 2b [16]. As the amplitude and direction of a V-saccade to the center of the lesioned site remain unaffected, saccades to locations around that site are directed away from the lesion. For example, V-saccades for sites rostral to the lesion undershoot the target, while V-saccades for sites caudal to the lesion will overshoot the target.
The simple vector-summation model of Eq 2b yields saccades that would always undershoot targets, as the lesioned population produces fewer output spikes than under normal control conditions. However, [9,19] observed that their estimate of the total number of spikes from the SC population, was remarkably constant, regardless of saccade amplitude, direction, or speed. Yet, they also observed that many cells in the normal SC fire some post-saccadic spikes. They therefore assumed that saccades are actively terminated by a downstream mechanism, whenever the criterion of a fixed number of spikes, N TOT , is reached: They demonstrated, by simulating the summation model of Eq 2b with actual recordings from �150 cells, that by including the criterion of Eq 4 (which constitutes a cut-off nonlinearity in the model), the pattern of saccadic over-and undershoots to a focal SC lesion can be fully explained. In addition, the extended summation model of Eqs 2b and 4 also accounts for weighted averaging of double-target stimulation in the motor map [10,28,29]. Moreover, although the vector-averaging model (Eq 2a) correctly predicts the pattern of saccadic dysmetrias, it fails to explain the substantial slowing of the lesioned saccades [16]. As this latter observation is also accounted for by Eqs 2b and 4 [9], it further supports the hypothesis that the SC population encodes both the saccade-vector, and its kinematics [24].

Electrical microstimulation in SC
Interestingly, electrical microstimulation experiments have also shown that at low current strengths, just around the threshold, the evoked saccade vectors become smaller and slower than main sequence [15,30]. These results do not follow from vector averaging (Eq 2a, which would always generate the same saccade, but might be predicted by dynamic summation (Eqs 2b and 4), if low-amplitude electrical stimulation were to recruit a smaller number of neurons at lower firing rates. However, if supra-threshold microstimulation would produce a large square-pulse population profile around the electrode tip (mimicking the profile of the imposed current pulses, as is typically assumed), the summation model would generate severely distorted saccade-velocity profiles, which are not observed in experiments. Yet, little is known about the actual activity profiles in the motor map evoked by electrical microstimulation, as simultaneous multi-electrode recordings in the SC during microstimulation are not available and would be obscured by the large stimulation artefacts [31].
Under microstimulation, two factors contribute to neuronal activation: (1) direct (feedforward) current stimulation of cell bodies and axons by the stimulation pulses of the electrode, and (2) synaptic activation through lateral (feedback) interactions among neurons in the motor map. How each of these factors contributes to the population activity in the SC is unknown. It is conceivable, however, that current strength falls off rapidly with distance from the electrode tip (at least by �1/r 2 ), and that hence a relatively small number of SC neurons would be directly stimulated by the electric field of the electrode.
Indeed, a two-photon imaging study, carried out in cortical tissue from rodents and cat are V1, showed that microstimulation at physiological current strengths directly activates only a sparse set of neurons directly around the immediate vicinity of the stimulation site [32]. These considerations therefore suggest that the major factor in explaining the effects of microstimulation in the SC motor map may be synaptic transmission through lateral excitatory-inhibitory connections among the cells. Such a functional organization in the SC is supported by anatomical studies [33,34], by electrophysiological evidence [35][36][37], and by pharmacological studies [38].

Spiking neural network model
We recently constructed a biologically plausible, yet simple, spiking neural network model for ocular gaze-shifts by the SC population to visual targets [39]. This minimalistic (one-dimensional) model with lateral interactions can account for the experimentally observed firing properties of saccade-related cells in the gaze-motor map [9,19], by assuming an invariant spiking input pattern from sources upstream from the motor map (e.g., FEF).
We here extended that simple network model to the full two-dimensional network map that accounts for microstimulation results over a wide range of stimulation parameters. To simplify the analysis of the network properties, and to limit the number of independent parameters that describe the electrical stimulation pulses, we used rectangular current profiles with different heights (current intensities) and durations. In line with the evidence from previous work, the network was tuned such that microstimulation provides an initial seed that directly activates only a small set of SC neurons, which subsequently sets up a large SC population activity through lateral synaptic interactions. Our results show that stimulating the network indeed sets up a near-normal population activity profile that generates appropriate saccadic command signals across the two-dimensional oculomotor range through the linear dynamic summation mechanism of Eq 2b.

Log-polar afferent mapping
The afferent mapping function (Eq 1a) translates a target point in visual space to the anatomical position of the center of the corresponding Gaussian-shaped population in the SC motor map. It follows a log-polar projection of retinal coordinates onto Cartesian collicular coordinates (Eq 1a; [13]). To allow for a simple 2D matrix representation of the map in our network model, we simplified the afferent motor map to the complex logarithm: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi with B u = 1 mm and B v = 1 mm/rad (isotropic map). Thus, the contribution, m, of a single spike at site (u, v) to the eye movement is computed from the efferent mapping function as: We thus constructed a spiking neural network model as a rectangular grid of 201 x 201 neurons. The network represents the gaze motor-map with 0 < u < 5 mm (i.e., up to amplitudes of 148 deg), and −π/2 < v < π/2 mm. The network generates saccadic motor commands of different directions and amplitudes into the contralateral visual hemispace through a spatial-temporal population activity profile. The location of the population in the motor map determines the direction and amplitude of the saccade target, whereas the temporal activity profile encodes the eye-movement kinematics, through Eq 2b. As described below, and in our previous study [39], the eye-movement main-sequence kinematics result from location-dependent biophysical properties of the neurons within the map, together with their lateral interconnections.

AdEx neuron model
We investigated the dynamics of the network model numerically through simulations developed in C++/CUDA [40]. The motor map is represented as a rectangular grid of neurons with a Mexican hat-type pattern of lateral interactions. The neural activities were simulated by custom code utilizing dynamic parallelism to accelerate spike propagation on a GPU [41]. The code was developed and tested on a Tesla K40 with CUDA Toolkit 7.0, Linux Ubuntu 16.04 LTS (repository under https://bitbucket.org/bkasap/sc_microstimulation). Simulations ran with a time resolution of 0.01 ms. Brute-force search and genetic algorithms, described below, were used for parameter identification and network tuning since there exists no analytical solution for the system.
The neurons in the network were described by the adaptive exponential integrate-and-fire (AdEx) neuron model [42], which accommodates for a variety of bursting dynamics with a minimum set of free parameters. The AdEx model is a conductance-based integrate-and-fire model with an exponential membrane potential dependence. It reduces Hodgkin-Huxley's model to only two state variables: the membrane potential, V, and an adaptation current, q. The temporal dynamics of the system are given by the following differential equations for neuron n: where C is the membrane capacitance, g L is the leak conductance,E L is the leak reversal potential,η is a slope factor, V T is the neural spiking threshold, q n is the adaptation time constant,a is the sub-threshold adaptation constant, and I inp,n is the total synaptic input current. In our previous paper [39] the input-layer of Frontal Eye Field (FEF) neurons had identical biophysical properties, and only received a fixed external input current, I inp,n = I ext . In the present simulations, we did not include a FEF input layer, as the electrical stimulation was applied within the SC motor map as an external current. Two parameters specify the biophysical properties of the SC neurons: the adaptation time constant, τ q,n (which is assumed to be location dependent), and the synaptic input current, I inp,n = I syn,n + I E (where I syn,n is a location-and activity-dependent synaptic current, and I E is the applied microstimulation current). Both variables change systematically with the spatial location of the cells within the network (rostral to causal). The remaining parameters, C, g L , E L , η, V T and a, were tuned such that the cells showed neural bursting behavior (see Table 1 for the list and values of all parameters used in the simulations, and The AdEx neuron model employs a smooth spike initiation zone between V T and V peak , instead of a strict spiking threshold. Once the membrane potential crosses V T , the exponential term in Eq 7a starts to dominate and the membrane potential can in principle increase without bound. We applied a practical spiking ceiling threshold at V peak = −30 mV for the time-driven simulations. For each spiking event at time τ, the membrane potential is reset to its resting potential, V rst , and the adaptation current, q, is increased by b to implement the spike-triggered adaptation: After rescaling the equations, the neuron model has four free parameters (plus the input current) [43]. Two of these parameters characterize the sub-threshold dynamics: the ratio of time constants, τ q /τ m (with the membrane time constant τ m = C/g L ) and the ratio of conductances, a/g L (a can be interpreted as the stationary adaptation conductance). Furthermore, the resting potential V rst and the spike-triggered adaptation parameter b characterize the emerging spiking patterns of the model neurons (regular/irregular spiking, fast/slow spiking, tonic/phasic bursting, etc.).

Current spread function
We applied electrical stimulation by the input current, centered around the site at [u E , v E ], according to Eq 5. We incorporated an exponential spatial decay of the electric field from the tip of the electrode: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi with λ (mm −1 ) a spatial decay constant, I 0 the current intensity (in pA), and a rectangular stimulation pulse given by P(t) = 1 for 0 < t < D S , and 0 elsewhere. Thus, only a small set of neurons around the stimulation site will be directly activated with this input current (see Results). Throughout this paper, we used a fixed input current profile (I 0 = 150 pA), λ = 10 mm −1 and D S = 100 ms) except for the final section, where we explore the effect of changing the microstimulation parameters on the resulting saccade. These parameters were determined by the neural tuning of the AdEx neurons in their bursting regime (see Neural tuning and bursting mechanism section in Results).  For simplicity, we incorporated a single rectangular stimulation pulse, P(t), rather than a train of narrowly spaced stimulation pulses. A train of pulses would introduce additional parameters, like pulse height, pulse duration, pulse intervals, pulse polarity, and number of pulses (stimulus duration), each of which would affect the network response. We have shown before that the spiking neural network model with AdEx neurons and lateral interactions can deal with such spiking input patterns [39]. However, varying these different stimulations parameters would complicate the analysis, and is deemed a topic for future work (see Discussion). Note also that the AdEx neurons act as 'leaky integrators' for membrane potentials below V T . Therefore, a sequence of pulses and a single rectangular pulse yield qualitatively similar membrane responses.
Remark on the current scale. In SC microstimulation experiments, one typically applies extracellular currents in the micro-Ampère range (10-50 μA) to evoke a saccade. In our simulations, we instead take the effective intracellularly applied current, which amounts to only a tiny fraction of the total extracellular current leaving the electrode.

The SC model: Synapses and lateral connections
The total input current for an SC neuron, n, located at (u n , v n ), is governed by the spiking activity of surrounding neurons, through conductance-based synapses, and by the externally applied electrical stimulation input (Eq 9): where g exc n and g inh n are excitatory and inhibitory synaptic conductances acting upon neuron n, E e and E i are excitatory and inhibitory reversal potentials respectively. These conductances increase instantaneously for each presynaptic spike by a factor determined by the synaptic strength between neurons, and they decay exponentially otherwise, according to: with τ exc and τ inh , the excitatory and inhibitory time constants; w exc i;n and w inh i;n are the intracollicular excitatory and inhibitory lateral connection strengths between neuron i and n, respectively (Eqs 12a and 12b) and τ i,s is the spike timing of the presynaptic SC neurons that project to neuron n. With conductance-based synaptic connections, spike propagation occurs in a biologically realistic way, since the postsynaptic projection of a presynaptic spike depends on the instantaneous membrane potential of the postsynaptic neuron. In this way, the state of a neuron determines its susceptibility to presynaptic spikes.
We incorporated a Mexican hat-type lateral connection scheme in the model, where the net synaptic effect is given by the difference between two Gaussians [44]. Accordingly, neurons were connected with strong short-range excitatory and weak long-range inhibitory synapses, which implements a dynamic soft winner-take-all (WTA) mechanism: not only one neuron remains active, but the "winner" affects the temporal activity patterns of the other active neurons. The central neuron governs the population activity, since it is the most active one in the recruited population. As a result, all recruited neurons exhibit similarly-shaped bursting profiles as the central neuron, leading to synchronization of the spike trains within the population [39]. Two Gaussians describe the excitatory and inhibitory connection strengths between collicular neurons as function of their spatial separation: where � w exc > � w inh and � s inh > � s exc , and s n is a location-dependent synaptic weight-scaling parameter, which accounts for the location-dependent change in sensitivity of the neurons due to the variation in adaptation time constants.

Network tuning
Electrophysiological experiments have indicated that the neural responses are well characterized by four principles: (i) a fixed number of spikes for each neuron associated with its preferred saccade vector N u,v ffi 20 spikes, (ii) a systematic dependence of the neuron's cumulative spike count on the saccade vector (dynamic movement field), N u,v (R, ϕ, t), (iii) scaled and synchronized burst profiles of the neurons in the population, resulting in a high cross-correlation, C pop (f n (t), f m (t)) � δ nm , between the firing rates of recruited neurons, and (iv) a systematic decrease of the peak firing rate of central neurons in the population, F peak , along the rostralcaudal axis, together with an increase of burst duration, T burst , and burst skewness, S burst . [19] argued that these properties follow from a systematic tuning of the gaze-motor map, and that they are responsible for the observed saccade kinematics. Here we applied these principles to determine a similarity measure between our simulated responses, and the experimentally recorded gaze motor-map features. In our network model, these features emerge from the interplay between intrinsic biophysical properties of the SC neurons, and the lateral interactions between them.
Distinct biophysical properties. The intrinsic biophysical properties of the neurons were enforced by systematically varying the adaptation time constant, τ q,n , and the synaptic weightscaling parameter, s n , in the motor map. Changes in the adaptive properties of the neurons result in a varying susceptibility to synaptic input. The synaptic weight-scaling parameter corrects for the total input activity. These distinct biophysical properties capture the systematically changing firing properties of SC cells along the rostral-caudal axis of the motor map, while keeping a fixed number of spikes for the neurons' preferred saccades N u,v (R, ϕ). Following the brute-force algorithm from our recent paper [39], the location-dependent [τ q,n , s n ] value pairs for the neurons were fitted to ensure a fixed number of spikes per neuron under a given microstimulation condition, and the subsequent excitation through lateral interactions (see below, Eqs 15 and 16). These parameters were first tuned for isolated neurons. The lateral interactions ensured that the bursting profiles in the population remained scaled versions of each other and had their peaks synchronized (evidenced from a high cross-correlation, C pop , between the burst profiles across the population). The s n values of Eqs 12a and 12b were scaled by the number of neurons in the population.
Lateral connectivity. The single-unit recordings also suggested that for each saccade the recruited population size, and hence its total number of spikes, is invariant across the motor map. The widths of the Mexican-hat connectivity (σ exc and σ inh ) govern the spatial range of a neuron's spike influence in the network, and directly affect the size of the neural population. In our model, these widths were fixed, such that they yielded local excitation and global inhibition. The connection strengths ( � w exc and � w inh ), on the other hand, affect the spiking behavior and local network dynamics, as they control how much excitation and inhibition will be received by each single neuron, and transmitted to others, based on the ongoing activity. Strong excitation would result in an expansion of the population, whereas a strong inhibition would fade out the neural activity altogether. Thus, balanced intra-collicular excitation and inhibition would be required to establish a large, but confined, Gaussian population. The parameters for the lateral connection strengths were found by a genetic algorithm, as described in our previous paper (Kasap and Van Opstal, 2017). In the current model we used eight saccade amplitudes for each generation to calculate the fitness of each selection (selected as R = [2,3,5,8,13,21,33,55] deg, and ϕ = 0 deg, to cover equidistant locations on the rostral-to-caudal plane: u = [0.69, 1.08, 1.60, 2.07, 2.56, 3.04, 3.49, 4.00] mm, and v = 0 mm, respectively).
The genetic algorithm minimized the root-mean squared errors (RMSE) between the spiking network responses and the rate-based model of [45]: from the fitness evaluation for each generation, we calculated the RMSE between the peak firing rates, F peak ; the number of elicited spikes from the central cells in the population, N u,v (R, ϕ); burst durations, T burst ; and burst skewness, S burst . Furthermore, the cross-correlations, C pop , between all active neurons and the central cell were included too to ensure that the experimentally observed gaze-motor map characteristics were taken into account for parameter identification. The fitness function was defined by a weighted RMSE summation: where the weights (0.1, 10, 103) were empirically chosen to cover similar ranges, since the F peaks vary from roughly 430-750 spikes/s, the number of spikes varies between 18 and 22, and the cross-correlation values are < 1. Peak firing rates of the central neurons from each population were calculated by convolving the spike trains with a Gaussian kernel (Eq 3; 8 ms kernel width), to determine spike-density functions of instantaneous firing rate. RMSE values for F peak along the rostral-caudal axis of the motor map were subsequently tuned by approximating the following relation: F peak ðrÞ ¼ F 0 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where F 0 = 800 spikes/s and β = 0.07 ms/deg (taken from [45]. The RMSE of the total spike counts during the burst from the central cells in the population were tuned to N u,v = 20 spikes, and was required to be independent of the neuron's position in the map. Synchrony of the neural activity within the recruited population was quantified by the RMSE of deviations for the cross-correlations between the central cell and all other active cells in the recruited population.

Generating eye movements
Eye movements were generated by the population activity following the linear ensemble-coding model of Eqs 2b and 3. We applied the two-dimensional efferent motor map of Eq 5. For any network configuration throughout this paper, the unique scaling factor of the efferent motor map (z) was calibrated for a horizontal saccade at (x, y) = (21,0) deg. The resulting eyedisplacement vector,SðtÞ, was calculated from the spike trains by interpolation with a first-order spline to obtain equidistant time samples. The interpolated data were further smoothed with a Savitzky-Golay filter, to obtain smooth velocity profiles.

Neural tuning and bursting mechanism
Fig 1 shows the membrane potential traces for three model neurons, differing in their adaptation time constants, τ q , which were stimulated under different microstimulation paradigms. The electrical stimulus strength increased from a low amplitude (I 0 = 50 pA; light blue traces) to a high intensity (I 0 = 250 pA, dark-blue traces), for stimulation durations between 25 and 225 ms. Note that for these different microstimulation regimes, the burst onsets and burst shapes (i.e., the instantaneous firing rates) could differ, even when the number of elicited spikes would be the same. These responses illustrate how the biophysical properties of the neurons affected their bursting behavior. First, the neuron could respond after the stimulation had terminated. Such a feature, as well as the bursting behavior, is only captured by more complex spiking neuron models. Even when the input current amplitude cannot drive a neuron rapidly to its first spike to initialize the burst (light traces), it suffices if the neuron's membrane potential crosses a certain threshold (V T in the AdEx neuron). The neuron can then elicit a spike after the stimulation is over (visible for stimulation durations < 75 ms).
Second, the stimulation amplitude determines the response onset: as the amplitude increases, the first spike occurs earlier. Such a behavior is to be expected, since the neuron model acts as an integrator [30]; higher input currents thus drive a neuron faster to its spiking threshold.
Third, the different neurons respond differently to long stimulation trains (> 175 ms). While the neuron with a longer adaptation time constant (τ q = 84.6 ms; Fig 1A) responds with repetitive bursts of 4 to 5 spikes, separated by a silent period, the faster recovering neuron (τ q = 52.4 ms; Fig 1C) elicits more and more spikes after the initial burst, especially for the higher current amplitudes (dark traces).
Interestingly, the neurons with the intermediate (Fig 1B) and short (Fig 1C) adaptation time constants switch between different bursting behaviors as the current amplitude increases along with longer stimulation durations. Regular short bursts with silent periods in between result from the slow decay of the adaptation current, which acts on the membrane potential as an inhibitory current. Hence, the adaptation time constant determines how fast a neuron will recover after each spike in a burst. Therefore, the strongly adapting neuron with a long will require more input current to elicit another spike (Fig 1A and 1B for stimulation duration >175 ms), and thus after the fourth spike in the burst, the adaptation current is already high enough to break the bursting cycle. The fast recovering neuron (Fig 1C, short τ q ) continues its burst with more spikes (dark traces at longer durations (B, C).
A phase plot of the instantaneous adaptation current vs. the membrane potential provides a graphical analysis of the effects of changing the neural parameters, the current input, and the initial state, on the evolution of the dynamical system. Fig 2 shows a number of phase-trajectories for the Adex model, for the parameters used in the simulations of the SC motor map. Nullclines illustrate the boundaries of the vector fields in the AdEx neuron's phase plane. The Vnullcline (V null ; i.e., dV/dt = 0 for Eq 7a) and the q-nullcline (q null ; i.e., dq/dt = 0 for Eq 7b) are shown as gray lines. Fixed points of the system lie at the intersections of these nullclines. A stable fixed point of the system is found at [-53 mV, 0 nA]. In all subfigures that is the starting point of the trajectories, and the state variables of the neurons will converge to this stable fixed point in the absence of input.
The q-nullcline follows a linear trajectory, whereas the V-nullcline represents a convex function because of the superposition of two V-dependent parts. For V < V T , the exponential term can be omitted and the linear V dependence will have a slope of g L . For V > V T , the exponential term will dominate with a sharp increase as V increases. When a neuron receives input, the V-nullcline shifts upward by as much as the current density, and the response of the neuron follows a trajectory on the phase plane toward the spiking threshold. The blue trajectories show the evolution of the state variables for three neurons with different τ q values, and stimulated at different current strengths. The horizontal arrows show the membrane potential in the spike initiation zone, V > V T . Spikes occur when the membrane potential overcomes the spiking threshold, V > V thr . After a spike, the membrane potential is reset, and the adaptation current is increased by b (Eq 7). The spiking threshold, V thr , and the reset potential, V rst , are indicated by the vertical dashed lines. With each spike, the adaptive current increases more and once it reaches values above the V-nullcline, the adaptive current is high enough to suppress the neuron from continued bursting, and hyperpolarizes.
In Fig 2A, the phase trajectory crosses values over V null = 150 pA after 5 spikes. Due to the hyperpolarization, the membrane potential starts to drop. The phase plot shows that the microstimulation is finished when the membrane potential decreases to -58 mV, and the smooth trajectory is seen disrupted. In Fig 2B, there is a second burst cycle since the microstimulation duration is much longer. After the first burst cycle crosses V null + 200 pA with 6 spikes (arrows are placed closer to V thr ), neuron follows the trajectory to the spike initiation zone for a second burst cycle with 5 spikes. The end of the microstimulation coincides with the second burst cycle and afterwards the membrane potential decreases fast due to the high adaptive current acting on the neuron. In Fig 2C, the neuron gets stuck in its first cycle and continues spiking repetitively. This pattern is due to the fast decay of the adaptive current, which drops by more than b after each spike. Therefore, the neuron would continue spiking repetitively, as long as the current is applied.
The neurons in the network were tuned to respond with a fixed number of spikes in a burst cycle (as in Fig 2A). This initial burst sets up a large population activity through the lateral connections. V null fluctuates for each neuron with the network dynamics, depending on the input from other neurons in the population. Microstimulation parameters were chosen such that the central neuron of the population would respond with a burst cycle of 4-5 spikes (typically, D S = 100 ms, and I 0 = 150 pA), independent of the biophysical properties of the neuron. To that end, the adaptation time constant, τ q,n , and the synaptic weight-scaling parameter, s n , for each neuron were determined by applying a fifth order polynomial fit to produce a fixed number of spikes (N = 20) for self-exciting neurons: s n ¼ ð8:808 � 10 À 9 � t 5 q;n À 3:280 � 10 À 6 � t 4 q;n þ4:855 � 10 À 4 � t 3 q;n À 3:607 � 10 À 2 � t 2 q;n þ1:383 � t q;n À 8:396Þ � 10 À 3 The self-excitation mimics the population activity, since the central cell's burst profile is representative for the entire population activity, due to burst synchronization across the active neurons. The adaptive time constant, τ q.n , varied from 100-30 ms in a linear way with the anatomical rostral-caudal location of the neurons, according to:

Microstimulation without lateral interactions
The current density drops rapidly with distance from the microelectrode tip, as given by the current spread function (Eq 9, with λ = 10 mm −1 , D S = 100 ms, and I 0 = 150 pA). Fig 3A illustrates this decay of current density on the motor map surface. The pulsed input current is presented onto the collicular surface at a site corresponding to the visual image point (u(R), v(ϕ) in Eq 5; Fig 3B and 3C). Microstimulation directly activated only a small set of neurons within a 250 μm radius. Fig 3B and 3C shows the number of spikes elicited by the activated neurons in the absence of intra-collicular lateral interactions. Each activated neuron elicited only 4-6 spikes within a given input duration range, regardless the electrode's location. These spikes arose from the initial bursting regime of the neurons until the adaptation current built up with repetitive spikes that canceled the microstimulation input (see Fig 2). The input amplitude affected the response delay of the neurons between stimulation onset and their first spike. Thus, in the model these small neuronal subsets generated only a brief pulse signal that is supposed to set up the entire population activity through lateral connections.

Including lateral interactions
We next tested the collicular network response to the same microstimulation parameters as in Fig 3, while including the lateral interactions. Fig 4A-4C shows the recruited neural population The peak firing rate of the central cells was close to 700 spikes/s and dropped in a regular fashion with distance from the population center. Note also that the cells near the fringes of the population were recruited slightly later than the central cells, but that their peak firing rates were reached nearly simultaneously. Moreover, the bursts all appeared to have the same shape. Fig 4C shows the saccade that was elicited by this neural population, together with its velocity profile. The saccade had an amplitude of 5 deg, reaching a peak velocity of about 200 deg/s. Fig 4D-4F shows the results for stimulation at the more caudal location in the motor map, yielding an oblique saccade with an amplitude of 31 deg. The size of the resulting population activity is very similar to that of the rostral population, and also the number of spikes elicited by the cells is the same. The peak firing rates of the neurons, however, were markedly lower, reaching a maximum of about 450 spikes/s. As a result, the burst durations increased accordingly, from about 50 ms at the rostral site, to more than 70 ms at the caudal site. Note that the saccade reached a much higher peak velocity (about 900 deg/s) than the smaller saccade in Fig  4C, but its duration was prolonged. Note also that the horizontal and vertical velocity profiles were scaled versions, indicating a straight saccade trajectory. These burst properties, which are due to a precise tuning of the biophysical cell parameters, underlie the kinematic main-sequence properties of saccadic eye movements [19,39,45]. Fig 6A shows the amplitudes and directions of 45 elicited saccades across the 2D oculomotor range (stimulation parameters: I 0 = 120 pA, D S = 100 ms). We avoided stimulating near the vertical meridian, as our model included only the left SC motor map (e.g., [15]), and stimulation at very caudal sites (R > 40 deg), where edge effects of the finite motor map would lead to truncation of the elicited population at the caudal end. Crosses indicate the coordinates of the corresponding motor map locations where stimulation took place; blue dots give the coordinates of the evoked saccade vectors. There is a close correspondence between the motor map coordinates and the elicited saccade vectors. Only for the most caudal sites the saccade vectors tended to show a slight undershoot. We have not attempted to compensate for these minor effects, e.g. by including heuristic changes to the efferent mapping function. The panels of Fig  6B and 6C show the evoked saccades for the nine stimulation sites along the horizontal meridian. Note that the saccade duration increased with the saccade amplitude, and that the peak eye velocity showed a less than linear increase with saccade size.  The main-sequence behavior of the model's E-saccades is quantified in Fig 8. Fig 8A shows the nonlinear amplitude vs. peak eye-velocity relationship, described by the following saturating exponential function: These main-sequence relations were combined into a single, characteristic linear relationship that captures all saccades, normal and slow ( Fig 8C)

Properties of electrically evoked eye movements
All three relations correspond well to the normal main-sequence properties, as have been reported for monkey and human saccades (e.g., [2]). Importantly, the main-sequence behavior of E-saccades was largely insensitive to the applied current strength as soon as it exceeded the stimulation threshold. This feature of the model is illustrated in Fig 9, which shows E-saccade peak eye-velocity as function of current strength for a fixed stimulation duration of D S = 100 ms (Fig 9A). The stimulation was applied at three different sites on the horizontal meridian (corresponding to R = 15, 21 and 31 deg). Below I 0 = 80 pA no movement was elicited, but around the threshold, between 90-120 pA, stimulation evoked slow eye movements, which eventually yielded the final amplitude ( Fig  9B). Immediately above the threshold at 130-140 pA, the evoked movement amplitudes and velocities reached their final, site-specific size (Fig 9A and 9B), which did not change with current strength over the full range between 140-220 pA. The associated peak eye velocity followed a similar current-dependent behavior for changes in stimulus duration (at a fixed current strength of 150 pA; Fig 9C). Thus, the quantity that determines evoked saccade initiation is the total amount of current (current amplitude times duration; e.g., [30]).

Summary
The simple linear ensemble-coding model of Eq 2b [9,45,46] seems inconsistent with the results of microstimulation, when it is assumed that (i) the rectangular stimulation input profile directly dictates the firing patterns of the neural population in the motor map, and (ii) that the neurons are independent, without synaptic interactions.
We here argued that these assumptions are neither supported by experimental observation, nor do they incorporate the possibility that a major factor determining the recruitment of SC neurons is caused by synaptic transmission within the motor map, rather than by direct activation through the electrode's electric field. We implemented circular-symmetric, Mexican-hat like interactions in a spiking neural network model of the SC motor map and assumed that the current profile from the electrode rapidly decreased with distance from the electrode tip ( Fig  3A). As a consequence, only neurons in the direct vicinity of the electrode were activated by the external electric field (Fig 3B and 3C; [31,32]).
Once neurons were recruited by the stimulation pulse, however, local excitatory synaptic transmission among nearby cells rapidly spread the activation to create a neural activity pattern which, within 10-15 ms, was dictated by the bursting dynamics of the most active central cells in the population (Fig 4). As a result, all cells yielded their peak firing rates at the same time, and the burst shapes of the cells within the population were highly correlated. Similar response features have been reported for natural, sensory-evoked saccadic eye movements [19], and it was argued this high level of neuronal synchronization ensures an optimally strong input to the brainstem saccadic burst generator to accelerate the eye with the maximally possible innervation.
Note that the evoked population activity does not grow without limit, but ceases automatically, both in its spatial extent, and in its bursting behavior, while the inhibitory currents acting on the neurons accumulate during the stimulation pulse. These currents are due to the synaptic far-range lateral inhibition, and to each neuron's own adaptive current. Thus, once the network is perturbed by an excitatory input current, the SC will set up a bursting population activity, without the need of an external comparator, or external feedback by a resettable integrator. Indeed, the adaptive current functionally acts as a putative 'spike counter' at the single neuron level. With this spiking neural network model, we thus offer an alternative framework for the oculomotor system, in which the SC motor map not only provides a spatial signal for the saccade vector, but also the instantaneous eye-movement kinematics, through the temporal organization of its burst profiles.

Network tuning
The site-dependent tuning of the biophysical parameters of the AdEx neurons, in particular their adaptive time constants and lateral-interaction weightings specified by Eqs 15 and 16, caused the peak firing rates of the cells to drop systematically along the rostral-to-caudal axis, while keeping the total number of spikes constant (Fig 5). As a result, the saccade kinematics followed the nonlinear main-sequence properties that are observed for normal (visuallyevoked) saccadic eye movements (Figs 6-8). In addition, the long-range weak inhibition ensured that the size of the population remained fixed to about 1.0 mm in diameter, and resulted to be largely independent of the applied current strength and the current-pulse duration (Fig 9).
The lateral excitatory-inhibitory synaptic interactions ensured three important aspects of collicular firing patterns that underlie the saccade trajectories and their kinematics: (i) they set up a large, but limited, population of cells in which the total activity (quantified by the number of spikes elicited by the recruited cells) can be described by a circular-symmetric Gaussian with a width (standard deviation) of approximately 0.5 mm (Fig 4A and 4D), (ii) the temporal firing patterns of the central cells (their peak firing rate, burst shape, and burst duration) solely depend on the location in the motor map (Eq 14), but the number of evoked spikes remains invariant across the map, and for a wide range of electrical stimulation parameters (Fig 5), and (iii) already within the first couple of spikes, the recruited neurons all became synchronized throughout the population, in which the most active cells (those in the center) determined the spike-density profiles of all the others (Fig 4B and 4E).
Here we described the consequences of this model on the ensuing kinematics and metrics of E-saccades as function of the electrical stimulation parameters. We showed that the network could be tuned such that stimulation at an intensity of 150 pA and a total input current duration of D S = 100 ms, sets up a large population of activated neurons, in which the firing rates resembled the activity patterns as measured under natural visual stimulation conditions. As a result, the kinematics of the evoked saccades faithfully followed the nonlinear main-sequence relations of normal, visually evoked saccades (Fig 8). Importantly, above threshold the saccade properties were unaffected by the electrical stimulation parameters (Fig 9).

Network normalization
Only close to the stimulation threshold, the evoked activity remained much lower than for supra-threshold stimulation currents, leading to excessively slow eye movements, that started at a longer latency with respect to stimulation onset. Similar results have been demonstrated in microstimulation experiments (e.g. [15,30]. The saccade peak eye velocity of the model saccades followed a psychometric curve as function of the amount of applied current (Fig 9). We found that the kinematics of the evoked eye movements at near-threshold microstimulation were much slower than main sequence (Fig 9). Although this property is readily predicted by the linear summation model (Eq 2b), it does not follow from center-of-gravity computational schemes (like Eq 2a), in which the activity patterns themselves are immaterial for the evoked saccade kinematics.
Conceptually, the lateral interactions serve to normalize the population activity. Therefore, the total number of spikes emanating from the SC population remains invariant across the motor map, and to a large range of (sensory or electrical) stimulation parameters at any given site. The nonlinear saturation criterion of Eq 4 is thus automatically implemented through the intrinsic organization of the SC network dynamics, and do not seem to require an additional downstream 'spike-counting' mechanism in order to terminate the saccade response, e.g. during synchronous double stimulation at different collicular sites (see, e.g. [28]).
Although other network architectures, relying e.g. on presynaptic inhibition across the dendritic tree, have been proposed to accomplish normalization of the population activity and vector averaging [28,45,[47][48][49], substantial anatomical evidence in the oculomotor system to support such nonlinear mechanisms is lacking. We here showed, however, that simple linear summation of the effective synaptic inputs at the cell's membrane, which is a well-recognized physiological mechanism of basic neuronal functioning, can implement the normalization when it is combined with excitatory-inhibitory communication among the neurons within the same, topographically organized structure. Such a simple mechanism could suffice to ensure (nearly) invariant gaze-motor commands across a wide range of competing neuronal inputs.

Further supporting evidence
Our model predicts near-normal activity profiles within the SC during microstimulation (Figs 4-6), and hence near-normal recruitment of the downstream brainstem circuits. Although simultaneous recordings in the SC during microstimulation are lacking, [50] described recordings from neural populations in the downstream brainstem burst generators (EBNs) and omnipause neurons (OPNs) during SC microstimulation. Their results indicated normal discharge patterns for OPNs and EBNs, and indistinguishable movement kinematics for stimulation-evoked and volitional saccades [51]. These results are nicely in line with the predictions or our model (Figs 8 and 9), at least for suprathreshold stimulation levels [26].

Future work
The two-dimensional extension of our model is a substantial improvement over our earlier one-dimensional spiking neural network model [39]. It can account for a much wider variety of neurophysiological phenomena. Yet, we have not attempted to mimic every experimental result of microstimulation. A few aspects in our model have not been incorporated yet, or some of its results seem to deviate slightly from experimental observations, which we briefly summarize here.
First, although the network output is invariant across a wide variety of stimulation parameters, and evoked saccade kinematics drop markedly around the threshold (Fig 9), the present model did not produce small-amplitude, slow movements near the stimulation threshold. This behavior has sometimes been observed for near-threshold stimulation intensities [15,30]. In our model, the saccade amplitude behaved as an all-or-nothing phenomenon (Fig 9B), which is caused by the strong intrinsic mechanisms that keep the number of spikes of the central cells fixed. Although we have not tested different parameter sets at length, we conjecture that a major factor that is lacking in the current model is the presence of intrinsic noise in the parameters and neuronal dynamics that would allow some variability of the evoked responses for small inputs. When near the threshold the elicited number of spikes starts to fluctuate, and becomes less than the cell's maximum, the evoked saccades will become smaller (and slower) too. Such near-threshold responses would also explain the truncated saccades generated when stimulation train durations are shortened [26].
Second, although the main-sequence relations of the model's E-saccades (Eqs 17 and 19) faithfully capture the major kinematic properties of normal eye movements, the shape of the evoked saccade velocity profiles were not as skewed as seen for visually-evoked saccades. As a result, the peak velocity is not reached at a fixed acceleration period, but at a moment that slightly increased with the evoked saccade amplitude (Fig 6C). We have not attempted to remediate this slight discrepancy, which in part depends on the applied spike-density kernels (here: Gaussian, with width σ = 8 ms, Eq 3), and in part on the biophysical tuning parameters of the AdEx neurons. However, it should also be noted that a detailed quantification of E-saccade velocity profiles, beyond the regular main-sequence parametrizations [15,30], is not available in the published literature. It is therefore not known to what extent E-saccade velocity profiles and V-saccade velocity profiles are really the same or might slightly differ in particular details.
Third, as explained in Methods, the electrical stimulation inputs were described by simple rectangular pulses, rather than by a train of short-duration stimulation spikes, in which case also the pulse intervals, pulse durations, pulse heights, and the stimulation frequency would all play a role in the evoked E-saccades [26,30]. We deemed exploring the potential results corresponding to these different current patterns as falling beyond the scope of this study, which merely concentrated on the proof-of-principle that large changes in the input for the proposed architecture of a spiking neural network led to largely invariant results. Note, however, that in our previous paper [39] the presumed input from FEF cells to the SC motor map did indeed provide individual spike trains to affect the SC-cells. We there demonstrated that the optimal network parameters could be found with the same genetic algorithm for such spiky input patterns, as applied here (Eq 13). The small differences in neuronal tuning parameters for the 1D model with FEF input, compared to the 2D model tuned to electrical pulse input, are mostly due to these fundamentally different input dynamics.
Fourth, [14] recently reported an asymmetric, anisotropic representation in the afferent mapping for the upper vs. lower visual hemi-fields, that would explain kinematic differences between upward vs. downward saccades. The underlying mechanism for this anisotropy is not yet clear. For example, it could result from (i) differences in lateral interaction strengths for up vs. down, thus creating different population profiles in the SC; (ii) differences in cell density along the medial-lateral SC coordinate, or (iii) systematic differences in the efferent projection strengths from medial-lateral SC neurons to the up-and down burst generators. In principle, our model could accommodate an anisotropic organization for upward vs. downward saccades by incorporating parametric changes at any of these levels. Here, we focused on a simple scheme, in which the SC was taken fully isotropic (Eqs 5 and 6), and the horizontal/vertical burst-generating circuits in the brainstem, including the horizontal/vertical ocular plants, were taken identical [9]. This ensured perfectly straight saccade trajectories in all directions, with homogeneous main-sequence properties, due to a full cross-coupling between the horizontal and vertical movement components ('component stretching'; see Fig 7).
Any change in this organization (e.g. more realistic eye-position related differences in the oculomotor plants, or different gains and delays in the up-vs. down vs. horizontal burst generators) will cause saccade trajectories to become curved, and direction and eye-position dependent, and may be made to resemble more closely the idiosyncratic differences observed in measured oblique saccades (e.g. [5]). Although an interesting topic, working out these many different factors, however, falls beyond the scope of this paper.
Fifth, double-stimulation experiments at different sites within the SC motor map have shown that the resulting saccade vector appears a weighted average between the saccades evoked at the individual sites [10,27]. In the present paper, we have not implemented double stimulation, although an earlier study had indicated that Mexican-hat connectivity profiles in the motor map effectively embed the necessary competition between sites to result in effective weighted averaging [28]. In a follow-up study, we recently explored the spatial-temporal dynamics of our model to double stimulation at different sites, and at different stimulus strengths [52]. Indeed, double stimulation results in weighted-averaged saccade responses, even when the SC activity is decoded by a dynamic linear-ensemble coding scheme, and without the need to implement an explicit cut-off on the total spike count, like in Eq 4. Thus, our SC scheme with excitatory-inhibitory interactions results to automatically normalize the total activity within the SC motor map (see also above). Hence, double stimulation results do not support the vector averaging scheme per se, as they can be explained by linear summation, in combination with intracollicular interactions, as well.
Finally, close inspection of the burst profiles in Fig 1 (showing stimulation results for single, isolated neurons) suggests that prolonged stimulation at sufficient current intensities could in principle generate multiple bursts of activity in the SC cells. For example, the top-left trace (I 0 = 250 pA, D S = 225 ms) shows a burst of 6 spikes, followed by a second burst of 5 spikes about 150 ms later. In principle, each of these bursts could be part of its own saccade, provided that the total network dynamics (including the lateral interactions) would preserve these properties. Indeed, the literature has shown that prolonged stimulation can lead to a series of eye movements of decreasing amplitude in the same direction (a so-called 'staircase' of saccades; [10,50,51]). Here we haven't tested our network for its potential to generate staircases, as we limited the stimulation durations to 250 ms. We suspect that the inhibitory currents and neural recovery may have to be balanced better to allow the prolonged input current to overcome the dynamic inhibition. Yet, although our network was not a priori designed for these staircases, their occurrence would be an interesting emerging property of the model.