The main anatomical and physiological features of the OI are summarized in Figures 1A–C. The OI is located in the hindbrain and composed of two bilaterally symmetric neuronal nuclei (Figure 1A). In goldfish, these nuclei consist of around N = 40 neurons (Pastor et al., 1994; Aksay et al., 2000), and are referred to as Area I. Quick horizontal eye movements (saccades) are generated by motoneurons in the abducens and oculomotor cranial nuclei (ABN, OMN, Figure 1A), which receive a velocity command from excitatory burst cells. The OI receives a collateral of the velocity signal and integrates it in the mathematical sense to create a position signal. The firing rates of the OI neurons are therefore linearly related to eye position (Figures 1B,C). The slope of their tuning curves is higher for neurons with more central eye position firing thresholds and lower for neurons with contraversive, peripheral eye position thresholds (Aksay et al., 2000), a property that has been called “recruitment order” (see Materials and Methods). Neurons from the same side are usually positively correlated, whereas neurons from opposite sides are negatively correlated, suggesting that the two sides are coupled by mutual inhibition and self-excitation (Aksay et al., 2003) as shown in Figure 1A. The firing rates of these “position” neurons remain stable in the absence of visual feedback and provide the signal that controls fixation of the horizontal eye position (Mensh et al., 2004).

These observations indicate that the network of position neurons can maintain a continuum of persistent firing rates over several seconds, i.e., the time scale of a typical eye fixation. Since only these persistent firing patterns are observed, the population activity of the two sides must be highly constrained. We use these constraints to reduce the network dynamics to the dynamics of the two interacting populations, which we describe by their summed activity (Figures 1D,F – see also Figures A1D–F). We will write X_L and X_R for the left- and right-side population output, measured as the resulting post-synaptic conductances (Figures 1E,F) and will refer to this population output as “population activity.”

As further discussed below, our modeling framework requires that the synaptic currents saturate to balance the progressive recruitment of neurons. Although such synaptic non-linearities are yet to be found in the OI, we here assume their presence. As discussed in the Materials and Methods, we choose step input-output functions to simplify the model tuning but can relax that assumption by using smoother sigmoidal functions (see Appendix). Given this choice, the synaptic output tuning curves are idealized, saturated versions of the firing rate tuning curves (Figures 1C,E). For simplicity of the model tuning, we also assume that the single-cell synaptic outputs have thresholds spanning the whole eye position range. This assumption can be reconciled with the data if the synaptic currents have different thresholds (Aksay et al., 2007). In consequence, the left- and right-side synaptic population outputs X_L and X_R resemble a staircase function where each step is caused by the synaptic input-output function of a single neuron (Figures 1E,F—see also Materials and Methods). Since there are many neurons involved, these population outputs approximate linear functions of the eye position.

The essence of our network models is shown in Figure 2 (see also Materials and Methods; a Matlab-based implementation of the models is available in the online supplementary information). In Figure 2A, we plot the population activities of Figure 1F against each other. Given the staircase shape of the population outputs, the resulting relationship is composed of multiple points on a line, which we here approximate by the red-blue line. This line illustrates the persistent firing states of the system, which range from virtual silence on the right and strong activity on the left side (X_R = 0, X_L large, red color) to virtual silence on the left and strong activity on the right (X_R large, X_L = 0, blue color). Different eye positions correspond to different points on this line (Figure 2B; see also below), and each point corresponds to a stable mode of firing for the network. The line is therefore often called a “line attractor” (Seung, 1996).

In the absence of synaptic input, neurons cease to fire given the neuronal leak. Therefore, for neural activities to be stable in a network in the absence of an external input, the recurrent network input to each neuron has to exactly match the neuronal leak (see Materials and Methods). An imbalance between those would cause either runaway excitation or inhibition. As the eye position moves from left to right, this balance of inputs and neuronal leaks needs to be maintained despite the non-linear distortions (such as thresholds or saturations) introduced by the biophysics of neurons and synapses. The resulting distortions can be compensated by the successive recruitment of neurons from one side and the suppression of neurons from the opposite side, as prescribed by the recruitment order (Figure 1C, see Materials and Methods, for details) (Seung, 1996; Seung et al., 2000; Aksay et al., 2007). It is important to note however that to obtain a line attractor in our network models, fine-tuning of the parameters is necessary, since changes in the parameters as small as 1% disrupt the line attractor dynamics.

Given the connectivity of the two sides of the system, the response of a neuron (on the right) is modeled as x_R = H(aX_R − cX_L + h), where parameter “a” determines the weight of the excitatory input from the right population X_R, parameter “c” determines the weight of the inhibitory input from the left population X_L, and parameter “h” models constant external inputs to the network. The function H(.) models the neuron's input-output function, and is either a threshold-linear function (in the case of a firing rate response) or a Heaviside function (in the case of a synaptic output response). The threshold of this input–output function corresponds to the set of points for which aX_R − cX_L + h = 0 (orange threshold line, Figure 2A). Any (X_R, X_L) combination that is below the orange line will make the right side neuron shown in Figure 2A fire.

The recruitment order fixes the thresholds of the neurons on the red-blue line as illustrated by the blue circles for a few example neurons (Figure 2A). However, the data do not specify how a neuron would respond to population activities X_R and X_L outside of this line, leaving a degree of freedom that is related to the relative strength of the self-excitatory input “a” and cross-inhibitory input “c” into each neuron (Figure 2A). Depending on how this threshold line is chosen for each neuron, the dynamics of the population activities outside of the line attractor change accordingly.

The dynamics of four exemplary models are illustrated in Figures 2C–F. Here, the arrows point in the direction in which the population activities evolve from different starting points. In Figure 2C, the line attractor is generated through mutual inhibition of the two sides. This configuration corresponds to one of the oldest models proposed for the OI (Cannon et al., 1983). In Figure 2D, the external excitatory connections are weakened, yet the resulting detrimental effect is compensated by weak self-excitation. This model was tested by Aksay et al. (2007), and is an extrapolation from the model in Seung et al. (2000): the authors here included mutual inhibition, while keeping orthogonal relaxation dynamics to the line attractor. We note that for both the mutual inhibition model and the weak self-excitation model, the dynamics outside of the attractor are orthogonal to the line, although the mutual inhibition model suggests faster dynamics (as indicated by the longer arrows).

The population dynamics of the model in Figure 2E were introduced by Aksay et al. (2007) to account for unilateral inactivation experiments in the goldfish. To obtain such dynamics in our modeling framework, the inhibitory and excitatory connections are set up so that each half of the oculomotor range is stabilized by an independent line attractor (ILA model). As a consequence, the population with the high activity (e.g., X_R) does not change its activity when the other side's population activity (X_L) is reduced. This situation is given when the left half of the system is silenced, which is equivalent to setting X_L = 0. Although the dynamics above the line attractor are unconstrained by experimental data, this model proposes that the dynamics above and below the line attractor are antiparallel to each other. We note that the ILA model captures the same population state space dynamics as the model suggested in Aksay et al. (2007), although the detailed implementation differs from the one in Aksay et al. (2007): the model does not incorporate input-output functions with high synaptic thresholds and uses a different distribution of tuning curves and cross-inhibition (see Materials and Methods).

A last example model is shown in Figure 2F. Here, the line attractor is stabilized through self-excitation only, and the inhibitory connections are non-existent. We note that this model is an extension of the model in Seung et al. (2000) from one population to two populations of excitatory neurons. In this case, any point outside of the line will be a potential stable fixed point as well. The system may still be confined to the line in practice, if the burst input during saccades always moves the population activities back onto the line.

The network models allow us to predict precisely how a perturbation would affect the system. Most importantly, these perturbations can be observed at the level of the eye position which makes the predictions experimentally accessible. To link the population activities to the eye position, we note that the iso-eye-position curves are likely passing through a non-linear transform introduced in the abducens nucleus (Figure 2B; discrepancy of position cell and motoneuron tuning curves—for details see Figure A3 and Materials and Methods). This bending of the curves also provides a simple explanation for the results of unilateral silencing of the OI in which the stabilization of eye position remains functional in only half the motor range and for roughly half the range of population activities (Aksay et al., 2007).

We can then simulate the response of the models to unilateral instantaneous inhibition and excitation, mimicking optogenetic stimulations with NpHR and ChR2, respectively. This idea is illustrated in Figure 3, where we focus on one of the network models, the ILA model (Figure 2E).

Figures 3A–D shows the effect of inhibiting the left half of the OI in the ILA model. Due to the extra inhibition, the left population activity X_L decreases immediately, as indicated by the arrows in Figure 3A. If the eye position before stimulation (initial eye position) is in the left range (black point), the inhibition causes the system state to cross θ isolines transversely so that the eye makes large movements to the right (orange point). If the initial eye position is in the right range, the system state shifts mostly in parallel to θ isolines, and the eye movements are small or null. After switching-off the inhibition stimulus, the system relaxes back to the line according to the dynamics of the intact system (Figure 3B). We note that if the initial eye position is in the left range, X_R increases, thereby moving the eye further to the midline. If the initial position is in the right range, X_L increases and X_R does not change after the stimulation is turned-off, so that the system returns to its initial state and the net eye movement is null (Figures 3A,B). We can extend this perturbation analysis to all initial eye positions, i.e., all points on the line attractor. Naturally, the results will depend on both the length and intensity of the stimulation. We used a brief stimulation (200 ms) and varying stimulation intensities. The net eye movement resulting from the combination of stimulation and relaxation was measured as the difference Δθ in eye positions just before the stimulation and 1 s after the stimulation (see simulations in Figure 3C). Large Δθ are observed when the initial eye position is in the left range, and negligible Δθ when the initial position is in the right range, (Figure 3D).

With similar reasoning, we can explore the system's response to excitatory perturbations (Figures 3E–H). During left excitation of the ILA model (Figure 3E), the left population activity, X_L, increases. After the stimulation, the system state relaxes back to the line attractor with the dynamics of the intact system (Figure 3F). If the initial eye position is in the left range, the value of X_L stays constant, and the right population activity, X_R, decreases (Figures 3F,G). Altogether, the eye makes therefore large movements to the left. If the initial position is in the right range, the system state moves mostly in parallel to θ isolines and the eye movements are small or null. Consequently, the perturbations Δθ of eye position occur mostly on the ipsilateral side to the stimulation, similar to the inhibitory perturbations, but with opposite sign (Figures 3G,H).

In both cases, the perturbations Δθ reflect the relaxation dynamics of the system, i.e., the dynamics of the intact system. By measuring these simulated perturbations for the full range of initial eye positions, we can cover all points of the line attractor. While we have illustrated these perturbations for the ILA model, we can perform similar predictions for the whole range of models. Conversely, we can measure the system's response to perturbations in optogenetic experiments, and then simply infer the dynamics of the system around the line attractor that are consistent with the experiments.

We note that we here modeled NpHR stimulations as divisive and ChR2 stimulations as additive. This distinction is based on electrophysiological recordings from the caudal zebrafish hindbrain (not limited to OI cells) which showed that NpHR stimulations induce a change in firing rate that is dependent on the initial firing rate, while for ChR2 stimulations no such effect was observed in the (small) range of firing rates tested (Figure A9). Therefore, for simplicity, we modeled the effect of ChR2 stimulation as being additive. However, assuming a subtractive influence of NpHR on population activity, or a multiplicative influence of ChR2 yields qualitatively the same results (Gonçalves, 2012) (simulation data not shown), and does not impede our ability to infer the overall dynamics from measurements.

To measure the effects of such instantaneous perturbations, and in turn infer the dynamics around the line attractor, we used fiber optic stimulations (Arrenberg et al., 2009) in behaving transgenic zebrafish (Figures 4A,B, Figure A7). Zebrafish are likely to have the same basic oculomotor circuit architecture and physiology as adult goldfish. The zebrafish larvae (5–8 days post-fertilization, dpf) were immobilized in agarose, and the agarose surrounding the eyes was removed. An optic fiber was positioned above the hindbrain to stimulate halorhodopsin (NpHR) or channelrhodopsin (ChR2). The transgenes Tg(UAS:ChR2(H134R)-mCherry)s1986t or Tg(UAS:NpHR-mCherry)s1989t were driven by the enhancer trap line Et(E1b:Gal4)s1101t, resulting in broad expression of ChR2 or NpHR in neurons (Scott et al., 2007; Arrenberg et al., 2009). To localize the zebrafish OI (Miri et al., 2011a) we performed unilateral light stimulations on different rostro-caudal positions of the hindbrain and measured the resulting eye drift magnitudes (Figures 4B–D). We chose short 200 ms stimulations to use the same time scale as burst signals which the position neurons typically integrate and to allow the NpHR-induced hyperpolarization to saturate. Two previous reports (Miri et al., 2011a,b) localized the integrator neurons in rhombomere 7 and 8 of the larval zebrafish, based on the corresponding location in the goldfish, two-photon laser ablations and on optogenetic loss-of-function experiments. In these laser ablation and optogenetic experiments, a single location was tested (rhombomere 7 and 8) and found to have an effect on the integrator performance. Here, we used small diameter optic fiber stimulations (50 μm) to test multiple positions (Figure 4C). The photoactivated volumes were columns of tissue approximately 10–15 cells wide and protruding through the entire dorsoventral extent of the hindbrain, as judged by Kaede photoconversion experiments (Figure A7). The maximal effect was observed around 50–150 μm caudal of the Mauthner cells, somewhat more rostral (rhombomere 5, 6, and 7) than in the previous reports (Miri et al., 2011a,b) (Figure 4D, Supplementary Movie 1).

We first focused on inhibitory (NpHR) perturbations. Unilateral stimulations about 100 μm caudal of the Mauthner cells (rhombomere 7 and 6) resulted in a drift of the eye position following a saccade (Figure 4E). As in the numerical simulations protocol, we computed the changes Δθ in eye position from just before the stimulation to 1 s after the stimulation. For simplicity, all unilateral stimulation results are plotted as left stimulations. Figures 5A,B show the results for a single fish: the eye positions are perturbed strongly in the left range and only weakly in the right range. In the left range, the more eccentric the initial eye position, the higher the elicited change in the eye position toward the right, as indicated by the positive Δθ. Additionally, the magnitude of this change increased with the stimulation intensity. These results become more distinct when averaging over all fish (n = 24) (see Figure 5C).

The results of the halorhodopsin activation experiments therefore agree with the ILA model that we used as an example network model in Figures 3A–D. We note two small differences to the predictions of the ILA model. First, even in the absence of NpHR stimulation, the eyes move slowly toward the null position from the whole motor range. Accordingly, the zebrafish eyes are slightly mistuned, whereas the ILA model is not. Second, the eyes show positive perturbation Δθ at the null position, unlike the model predictions (with Δθ = 0 for θ = 0). This difference could be due to uncertainties in the eye position: since the model only predicts positive perturbations, Δθ ≥ 0, any ambiguity in eye position, be it due to uncertainty of the null position (see Materials and Methods), measurement errors, or hysteresis in the system itself, will shift the average perturbation at θ = 0 to a positive value. Altogether, the dynamics of the OI below the line attractor are therefore well-captured by the ILA model (compare Figure 2E). The halorhodopsin inactivation experiments thereby confirm the pharmacological inactivation experiments of Aksay et al. (2007).

Following the same protocol as in the NpHR experiments, we next performed unilateral instantaneous stimulations (100 ms) in animals expressing ChR2 so that one side of the OI was excited (Figures 5D–F). Again, all unilateral stimulations were pooled and plotted as left stimulations. When the initial eye position was in the right range, the perturbations generally caused an eye movement toward the left, i.e., toward the null position, as indicated by the negative Δθ. In this case, we furthermore observed that more eccentric initial eye positions or higher intensity stimulations induced stronger changes in the eye position. On the other hand, when the initial eye position was in the left range, stimulations elicited very small or no changes toward the right after 1 s, when compared to the control case (Figures 5D–F).

For very high ChR2 stimulation intensities, we additionally noted eye movements toward the null position when the initial eye position was on the same side as the stimulation (ipsiversive eye positions; Figure A6). Since we are interested in inferring the dynamics in the immediate vicinity of the line attractor, we linearly regressed the magnitude of the drift at very low stimulation intensities (Figure 6A). In agreement with the data for medium intensities (Figure 5F), low intensity stimulation induces eye movements mostly during contraversive eye positions (when the initial eye position is on the other side as the stimulation; Figures 6B,C).

The results of the channelrhodopsin activation experiments therefore invalidate the ILA model, as we find that weak excitatory unilateral perturbations of the left OI induce centripetal (toward the null position) and not centrifugal eye movements. Furthermore, we find the strongest effect on the side contralateral to the perturbation, and almost no effect on the ipsilateral side. We note that the perturbations resulting from experimental excitation and inhibition are approximately point-symmetric with respect to the central eye position, (θ = 0, Δθ = 0), and not reflection-symmetric with respect to Δθ = 0, as proposed by the ILA model.

The results from the NpHR and ChR2 experiments provide us with the means to infer the OI dynamics around the line attractor. Since the inhibitory stimulations yield results similar to those predicted by the ILA model (compare Figure 3D and Figure 5C), we conclude that the dynamics below the line attractor are similar to those of the ILA model (Figure 2E). To obtain the observed point-symmetry of the perturbations (Figures 5C,F), however, we need to assume that the dynamics above the line attractor are roughly orthogonal to the dynamics below. Consequently, our experimental results suggest a model with dynamics as shown in Figure 7A. Here, the dynamics outside the line attractor are organized around a central point, corresponding to the null position of the eyes (the “null position” or NP model). Although Figure 7A illustrates the overall flow of the trajectories toward the line, we note that our line attractors are composed of individual fixed points, such that the fine-scale dynamics change in the vicinity of the line (see Figure A4C). Although the fine-scale dynamics predict small changes toward the left in the left range after ChR2 stimulation, these small changes do not appear in the predictions if we assume that the system is noisy.

These conclusions are further supported when we analyze the transient dynamics of the experimental eye movements around the NpHR or ChR2 stimulation period (Figures 7C–F). If the initial eye positions are in the left range, the NpHR stimulation generates monotonic eye movements toward the right, as indicated by the transient, positive derivatives of the eye movements (green traces in Figure 7C). If the initial eye positions are in the right range, the NpHR stimulation leaves them essentially unaffected. This contrasts with the ChR2 stimulations which transiently affect all eye positions. For initial eye positions in the left range, the stimulation causes transiently biphasic eye movements, as indicated by the negative, then positive derivatives of the eye movements (green traces in Figure 7E). In other words, the eye is briefly moved to a more eccentric position, before relaxing back to its original position. For initial eye positions in the right range, we observe an overall negative derivative of the eye movement, generating a centripetal eye position drift (yellow traces in Figure 7E).

These transient dynamics match with the population activity dynamics predicted by the NP model (Figures 7D,F). Let us focus on the excitatory perturbations. When the left OI is excited, the left population activity X_L increases in the whole motor range (Figure 7F). In the relaxation phase, if the initial eye position is in the left range, the system counterbalances the left stimulation by decreasing X_L without changing X_R. Accordingly, eye movements are biphasic and altogether eye position changes only marginally. In the right motor range, the right population activity decreases in the relaxation phase, while the left population activity stays intact, thus causing centripetal (to the null position) eye movements (Figure 7F).

Besides explaining the perturbation data well (Figures 7G,H), the resulting NP model is also in better agreement with the natural leakiness of the OI, i.e., the slow drift in eye positions that we observe even in the absence of perturbations (Figures 5C,F). First, we note that many (random) mistunings of the synaptic parameters will automatically lead to leaky eye positions. Second, even in the perfectly tuned NP model, the slow centripetal drift will emerge due to the presence of noise in the neuron's firing. Since any internally generated noise in the OI will have similar (if smaller) effects than our perturbations, such noise will cause slow centripetal, and thereby leaky eye movements (see Figure A4D).

While the dynamics that we inferred through the experiments seem perfectly reasonable, the applicability of the NP model to the OI hinges on the correctness of our experiments. One particular concern is that the limits of the integrator cells in the larval zebrafish are not well-established. Hence, we cannot rule out that we may be stimulating efferent motor neurons or afferent inputs (for instance, vestibular neurons) to the integrator. Since stimulation of motor neurons should not lead to persistent eye position changes, we ignore the potential stimulation of those neurons. However, stimulating afferent inputs to the integrator, in addition to the integrator cells remains a possibility. Using the modeling framework, we show below that accidental stimulation of these afferent inputs does not invalidate our conclusions.

We distinguish four possible scenarios for afferent stimulation (Figure 8A): in the first two scenarios, afferent stimulation leads to either ipsilateral excitation or inhibition of the integrator cells (violet and green arrows), and in that case our overall stimulation direction does not change in the synaptic population output space. In the other two scenarios, the afferent stimulation leads to either contralateral excitation or inhibition of integrator cells (red and blue arrows), therefore introducing a component to the stimulation which is orthogonal to the pure integrator stimulation.

For a cross-inhibitory architecture, the dynamics above the line attractor are either as in the ILA model (Figure 2E), as in the NP model (Figure 7A), or somewhere in between (in that case, the dynamics are roughly perpendicular to the line attractor, as in Figure 2D). We simulated all four afferent stimulation scenarios for each of the three models, and assumed that afferent stimulation has half the magnitude of the direct integrator stimulation. As can be seen in Figures 8B–E, neither the ILA model nor the intermediate model produce results similar to the data (regardless of the type of afferent stimulation). Only the NP model leads to the results observed in the data.

Furthermore, the NP model provides an explanation for the effect of the high intensity stimulation mentioned earlier (compare Figure A6): if at high intensity, unilateral ChR2 stimulation also excited the contralateral side via light scattering, then the NP model would indeed predict centripetal eye drifts on the ipsilateral side (Figure 8E, red curve).

In conclusion, the NP model is the only candidate model within our theoretical framework that agrees with all aspects of the data. Most importantly, these dynamics suggest a different connectivity than implied by the previously proposed models. The respective, effective connectivities for the NP model are shown in Figure 7B. The dynamics are mediated by weak mutual inhibition and a self-excitation whose effective strength depends on eye position, suggesting that both sides are almost, albeit not completely, decoupled, in line with previous findings in the goldfish (Debowy and Baker, 2011). To keep the self-excitation in check, external (e.g. vestibular) inputs to this model are inhibitory. A detailed description of the implementation of the network models can be found in the Materials and Methods.

We showed that a whole range of network models could account for the set of electrophysiological features that have been measured in OIs across animals of different species. These models differ from each other by the relative strength of self-excitation and mutual inhibition, and the respective dynamics prevailing in the population activity state space outside of the line attractor. To test these models and experimentally unveil the dynamics of the OI, we performed optogenetic perturbations in the larval zebrafish. Using the silencer NpHR, we found that unilateral light stimulation induced eye movements back to the midpoint, if the eye position prior to the stimulation was on the same side as the stimulation. The experimental results for NpHR in larval zebrafish were in accordance with results from goldfish using inactivation with lidocaine (Aksay et al., 2007), therefore corroborating the previously proposed ILA population dynamics. However, we found that unilateral ChR2 stimulations did not have the opposite effect to NpHR stimulations predicted by the ILA dynamics, i.e., centrifugal (away from the null position) eye movements on the stimulation side (ipsiversive eye positions). Instead, ChR2 stimulations had a centripetal (toward the null positions) effect on eye positions on the side opposite to the stimulation (contraversive eye positions).

Hence, perturbations always tend to drive the eye positions toward the midpoint, indicating that this point is the focus of the OI dynamics. Consequently, we inferred the dynamics around the line attractor from these experiments, and named the resulting model the null-position or NP model. This model suggests an OI architecture with strong self-excitation and weak cross-inhibition. Only a specific combination of excitation on one side and inhibition on the other will lead to actual changes in the position signal as required in saccadic eye movements (e.g., movements from one range to the opposite range) (Gonçalves, 2012). Interestingly, our results can at least in part explain the findings in pharmacological experiments where glutamate or GABA agonists were injected in the OI. Both an increase in excitation as well as an increase in inhibition resulted in centripetal eye movements, much as in our experiments (Arnold et al., 1999).

In the presented study we show how optogenetic experiments can be combined with modeling to infer the dynamics of a neural circuit module for integration. In the last few years, the field of optogenetics has provided a powerful set of techniques to perform gain- and loss-of-function experiments (reviewed in Zhang et al., 2007a; Luo et al., 2008; Fenno et al., 2011) and has been applied to zebrafish (Szobota et al., 2007; Douglass et al., 2008; Arrenberg et al., 2009; Baier and Scott, 2009; Zhu et al., 2009; Schoonheim et al., 2010). A fundamental problem in interpreting the effects of optogenetic stimulations is that a system's response is a combination of the stimulation magnitude and the intrinsic network dynamics. Indeed, the eye movements induced through the optogenetic perturbations depended on both the light-intensity, i.e., the strength of stimulation, and on the eye position, i.e., the internal state of the system prior to stimulation. To understand these interdependencies, we relied on network modeling (Seung, 1996; Seung et al., 2000; Goldman et al., 2003; Eliasmith, 2005; Aksay et al., 2007). In turn, the mismatches between the model predictions and the experimental results allowed us to constrain the class of feasible network models and thereby improve our understanding of the OI. This general approach illustrates the importance of the internal state of a system during a perturbation. Wherever this internal state is at least partially known, optogenetic perturbations can provide useful clues toward the underlying network dynamics.

In our experiments we made use of a Gal4 driver line that drove strong expression broadly in neurons. While local stimulation in this line led to marked behavioral changes, the manipulation was not specific to the neural integrator. The neural integrator in larval zebrafish is distributed across approximately 150 μm in the hindbrain of the larval zebrafish (Miri et al., 2011a), with non-integrator neurons interspersed between the cells with position signals. However, our results are not easily explained by a stimulation of these non-integrator neurons. First, both NpHR and ChR2 stimulation induce stable and persistent changes in eye position. This persistent change makes an influence of the motoneurons that lie in close rostral proximity to the integrator unlikely. Exclusive motoneuron stimulation should cause the eyes to move back to the original position immediately after stimulation offset, an effect we did not observe. Second, we may have stimulated cells which project to the integrator such as the saccade-generating neurons. However, saccade generating neurons are only active during saccades and unilateral stimulation is therefore expected to only change the saccade frequency (Schoonheim et al., 2010) and have no effect in-between saccades. Nevertheless, a low level stimulation of the saccade-generating neurons could have occurred without the generation of a measurable saccade: in that case, given that saccade-generating neurons excite the ipsilateral OI and inhibit the contralateral OI, a perturbation of these neurons would cause an indirect stimulation of the integrator neurons roughly in the same direction as the direct integrator stimulation, therefore not invalidating the interpretation of our results. Third, we may have affected some of the vestibular inputs to the OI. However, even in this scenario, our conclusions about the integrator dynamics hold up. Since vestibular inputs are included in the network models, we can simply simulate their accidental stimulation. As shown in Figure 8, accidental stimulation offsets the magnitude of induced eye drifts, but overall does not alter their eye position dependence. Within the range of models considered, the data can therefore only be explained by the NP model, but not by the other models.

We have stimulated excitatory and inhibitory integrator cells at the same time, which could potentially lead to unexpected network effects, e.g., due to induced imbalances of excitation and inhibition within the network. However, several observations support our interpretation of the data. First, electrophysiological recordings in the hindbrain of the same zebrafish lines used in this study suggest that more than 80% of NpHR expressing cells were significantly silenced during illumination and more than 90% of ChR2 expressing cells showed an increase in firing rate upon illumination (Arrenberg et al., 2009). Second, NpHR stimulation results are in agreement with the pharmacological inactivations in the goldfish integrator (Aksay et al., 2007), therefore confirming the inhibitory nature of the NpHR stimulations on the integrator. Third, ChR2 stimulation leads to different results than NpHR stimulation, which is proper inhibition. Fourth, both ChR2 and NpHR experiments were performed with the same Gal4 driver line. Therefore, if ChR2 stimulation led to net inhibition of the integrator, then NpHR stimulation would lead to excitation of the integrator which is ruled out based on the second observation. While these results suggest that optogenetic manipulations of the neural integrator changed the network activity in the expected direction, future zebrafish lines, e.g., with specificity for excitatory or inhibitory neurons, will facilitate the dissection of this circuit.

One could hypothesize that the effect of ChR2 stimulation saturates or reverses (depolarization block, Kleinlogel et al., 2011) with increasing stimulation intensities or for highly active cells. Our previous electrophysiological recordings provided no evidence for such an effect (Arrenberg et al., 2009, Figure A9). More importantly, this possibility is not supported by the eye movement data in the range of intensities analyzed, since the eye movements scale linearly with the stimulation light intensity (Figure 6B). Also, the effect does not reverse at high light intensities (Figure A6). As a final note, the modeling framework already assumes that the synapses of highly active neurons are saturated, so that, at least within the model, ChR2 stimulation does not affect these neurons.

In this study, we deliberately focused on the dynamics in the neighborhood of the stable eye position states. Consequently, our analysis was restricted to stimulation with low or medium light intensities. The NP model, however, makes predictions for any stimulation intensity, opening the question of what happens when the stimulation intensity is increased. As shown in the Figure A6, for high NpHR stimulation intensities, we additionally found small centripetal movements when the eye position was on the side opposite to the stimulation. For high ChR2 stimulations, we found centripetal eye movements when the eyes were on the same side as the stimulation. In the NP model, this would require a change in the dynamics far away from the line attractor, requiring the arrows to bend further toward the midline. While it seems unlikely that the effects at higher light intensities can be explained through scattering of light into the other hemisphere (data not shown), the effect could potentially be explained through strong stimulation of vestibular inputs. Furthermore, we notice that strong ChR2 stimulation could synchronize the activities of cells, which may have a range of effects, including complete shutting down of persistent activity (Dipoppa, 2012). We therefore refrained from including these observations in the model.

Integrators are ubiquitous in the brain and are involved in several important computations. For instance, in decision-making tasks requiring sensory integrations, neurons in the lateral intraparietal cortex behave similar to integrators (Gold and Shadlen, 2001; Wong et al., 2007). In working memory tasks, neurons in the prefrontal cortex exhibit almost linear dynamics during the times in which an animal needs to remember a stimulus, similar to integrators operating in several dimensions (Singh and Eliasmith, 2006; Machens et al., 2010). In the head direction system, a head velocity signal is integrated into head position (Zhang, 1996).

In previous line attractor models, it has generally been assumed that noise causes random drift along the line (Seung, 1996). While this is true in models with orthogonal dynamics around the line (such as the model illustrated in Figure 2D), in the case of the NP model the relaxation to the line has a preferred direction, therefore causing a systematic drift toward the null position. The term “line attractor” for the NP model is therefore strictly only valid in the limit of vanishing noise. For large noise levels, the model shows flow toward the central eye position with equivalent speed from every point in the state space. Hence, one could interpret the NP model as suggesting that the OI operates like a single fixed point, and not a line attractor, as extensively suggested in previous literature. However, for large noise levels, we can re-tune the NP model to recover the NP dynamics in the proximity to the line by implementing stronger cross-inhibition (simulations not shown). In any case, random perturbations of the NP model (such as noise) are unlikely to cause a centrifugal drift of the eye position.

Given the centripetal drift suggested by the NP model, we hypothesize that the OI features dynamics with a higher degree of built-in “safety” than previously thought. The OI has been observed to be leaky on longer time scales, both in goldfish and zebrafish (Mensh et al., 2004; Miri et al., 2011a). This leakiness may be a behaviorally advantageous feature, since, by bringing the system to the central position by default, it enables the relaxation of the eye muscles. Yet even higher brain systems may rely on such a built-in leakiness. In working memory tasks that employ graded persistent activity (Machens et al., 2005), for instance, a tendency to drift toward the central point while memorizing a sensory stimulus could explain the psychophysical errors that are known as contraction bias (Ashourian and Loewenstein, 2011).

While this built-in “safety” may help against noise in the system, it does not solve the fine-tuning problem, i.e., the instability of the line attractor against perturbations in the synaptic weights in the network. Indeed, this fine-tuning problem is a separate problem, somewhat orthogonal to the problems that we have investigated here, for which several solutions have been proposed (Koulakov et al., 2002; Goldman et al., 2003; Moreau and Sontag, 2003).

While we here have assumed that neural integration in the oculomotor system is generated through precise recurrent feedback in a neural circuit, in previous literature single-cell mechanisms have been put forward to explain neural integration observed in multiple areas in the brain. In particular, following an experimental demonstration of integration in individual cells from the entorhinal cortex (Egorov et al., 2002), a body of theoretical work has proposed several biophysical mechanisms which could underlie single cell integration, dispensing synaptic feedback (Loewenstein and Sompolinsky, 2003; Fransén et al., 2006). In the OI, unilateral disruption of the connectivity leads to neural activity drifts with time constants which are typically above 1 s (Pastor et al., 1994; Aksay et al., 2007), suggesting that single-cell mechanisms possibly play a role in the process of integration. However, single-cell mechanisms remain largely uncharacterized in the integrator, and therefore we here have followed the network mechanisms hypothesis as in previous studies of this system (Seung, 1996; Aksay et al., 2007). The contribution of single-cell mechanisms to the slow dynamics in the integrator is a challenge for future research.

We also note that our network model is a rate-model, in which the activities of individual cells are described by rates rather than precise spike times. Although we lose biophysical realism with this type of model, we gain analytical tractability, a very useful asset in interpreting experimental results within a theoretical framework, and in constructing models in accordance with data. Since the position cells exhibit persistent activity with regular firing (Aksay et al., 2003), temporal averages of spiking events are a good qualitative description of the system. Nevertheless, an equivalent model with spiking neurons could be built as shown in Seung et al. (2000), Eliasmith (2005).

Given the weak mutual-inhibition, the NP model behaves close to a system with a plane of stable fixed points (see Figure 2F), and shows slow dynamics around the line. We note that the relaxation of the eye positions after stimulation is indeed slow (on the order of 200 ms). However, these slow dynamics could be reflecting the dynamics of the muscle physics rather than the slow dynamics of the integrator. Future work should show whether such slow dynamics can indeed be observed.

Our network model features multiple stable fixed points, which suggests that the eye positions corresponding to these fixed points should be held comparatively longer than eye positions in-between the fixed points. However, in our data, the system seems to visit a range of eye positions in a homogeneous way (both during NpHR stimulation and during spontaneous, slow eye position decay), which contrasts with the prediction of our model. Such homogeneity could be due to external factors to the integrator, such as small saccadic commands causing smooth eye movement fluctuations, or the dynamics of the motor neurons and muscles. Nevertheless, the homogeneity found in the data challenges our hypothesis of discrete fixed points in the integrator and suggests further studies to elucidate this question.

Our modeling framework assumes a specific mapping from network activity to eye position, based on the difference between the tuning curves of position neurons and motor neurons (see Materials and Methods). Although this assumption is essential in our framework to account for the unilateral inactivation results in Aksay et al. (2007), one could relax it and assume linearity between population synaptic outputs and eye position by introducing high synaptic thresholds in the same fashion as in Aksay et al. (2007). However, we believe that the non-linear mapping is most likely present in the system and should be included in future modeling studies of the oculomotor system.

In a recent study (Miri et al., 2011a), zebrafish position neurons were shown to have variable timescales of integration, so that the associated relaxation time constants varied across neurons over one order of magnitude. This suggests that the dynamics of the OI could be high-dimensional (on the order of the number of neurons), in contrast with our line attractor model, which is implemented with homogeneous time constants across neurons and has low (2D) dimensional dynamics. Given that in this study we are interested in the dynamics of the population activities, the details of single-cell time constants do not affect our conclusions. In the future, it will be interesting to perform optogenetic stimulations and at the same time measure the activities of the position cells to explore the full state space, and realize the dimensionality of the system's dynamics.

Our models belong to a series of works suggesting that the OI builds up a line attractor by a balance between neural saturation and progressive recruitment of neurons to compensate such saturation (Seung, 1996; Seung et al., 2000; Aksay et al., 2007). Future work should specifically target the validity of this assumption as it is crucial for the whole modeling framework. Specifically, this assumption predicts that neurons have no responsibility or influence on eye positions that are below their firing threshold or that are sufficiently above that threshold (when they run into saturation). This prediction could e.g., be tested with single-cell ChR2 stimulation. Excitatory stimulation of different cells would then lead to movements in different eye ranges, enhancing the fact that different neurons are responsible for different stable activities on the motor range (Gonçalves, 2012). Consequently, combining single-cell optogenetics with the framework here designed has the potential to provide even deeper insights into the detailed structure of the OI in the future.

The files containing eye position trace and stimulation data were analyzed with custom algorithms using MATLAB. We only analyzed the first 3 h of each experiment, since in some recordings (without dish lid) water evaporation affected the camera focus and thereby the correct detection of eye position angles. The induced eye movements were measured as the change in eye position from stimulation start to 1 s after stimulation start. Most induced eye movements consisted of a single monotonic drift that was completed within 400 ms after stimulation. However, a 1 s interval was chosen to account for potential variations of the duration of the post-stimulation drift. The eye position change, averaged across left and right eye, was then plotted against the initial eye position, averaged across eyes. The results from each fish were fitted by a cubic smoothing spline with boundary second derivatives equal to zero, using the spline fit-package developed for MATLAB by Jonas Lundgren. This fitting procedure was repeated for several stimulation intensity bins. A population plot was generated by averaging the spline fits of all fish, for each intensity bin. We observed some variability of motor ranges across fish, especially in animals expressing ChR2, which tended to have a smaller motor range (average motor range and respective standard deviations in NpHR animals: [−20.0° ± 5°; 18.7° ± 4.4°]; in ChR2 animals [−13.7° ± 3.9°; 13° ± 3.5°]). Therefore, before averaging across animals, we normalized all motor ranges to generate a smooth population plot. In some fish, there was an undersampling of events at the eccentric eye positions, and therefore we removed the 5% most eccentric events on the right side and on the left side before averaging across animals. In several fish, some intensity bins did not span the full motor range, and we extrapolated the corresponding splines. To account for the uncertainty of the extrapolation, we performed a weighted average over fish for eccentric eye positions, with weights that decreased as the extrapolation distance increased.

The null-position of the eyes was defined as the average of the extreme eye positions. Note that this definition only works well for eyes which explore eye positions in both directions equally well, which was not the case for every fish. However, in our hands this was the best definition since other definitions (median, average, fixed at +8°/−8°) were more subjective or decreased the across-fish consistency of the results. An animal's distribution of post-saccadic eye positions was often bimodal, since the animal explored the peripheral positions more frequently. In some experiments, we stimulated multiple times (1, 3, 5, and 7 s) after a single post-saccadic fixation in order to access the full scope of eye positions.

The location of stimulation was in caudal proximity to the previously identified region containing saccade generating neurons (Schoonheim et al., 2010), and in a fraction of the events, the stimulation induced a saccade (more frequently for ChR2 stimulations than for NpHR). Since this study focuses on the performance of the integrator, we did not include these events in our analysis. In Figure 5, these were excluded by a velocity threshold (>20°/s). In Figure 4D, the data points result from eye velocity vs. eye position fits (linear regression through the origin, see Figure A8), and correspond to the inverse of the eye position drift time constant at a specific hindbrain location. For each fit, only the middle 70% of the eye velocity/eye position data points were fitted, thereby excluding outliers caused by saccades. For each eye, the eye drift was normalized across hindbrain positions. Each data point in Figure 4D corresponds to one fish (average of the normalized left and right eye, for individual eyes see Figure A8) and hindbrain position. The data points for Figure 4D (eye drift vs. hindbrain position) were then fitted by gaussian functions. For the Gaussian fit functions in Figure 4D, the mean was restricted to the interval [0; 200]. This modification was made because for approximately flat distributions (control data points ipsi ctrl, contra stim, contra ctrl) peripheral data points sometimes caused Gaussian fits with means far away from the tested hindbrain region.

In previous work (Gonçalves, 2012), we show how to build a line attractor model accounting for the main findings in the oculomotor integrator, using basic design principles that have previously been suggested in the literature (Seung, 1996; Eliasmith, 2005; Machens and Brody, 2008). We use these principles to clarify how the network connectivity is related to the tuning curve properties. Here, we summarize the main features of the model developed in Gonçalves (2012).

We start by describing the assumptions and simplifications underlying the network model. Most of these assumptions are adapted from the previous literature and are based either on arguments of plausibility and simplicity or on observations about the oculomotor integrator (Seung, 1996; Aksay et al., 2000; Seung et al., 2000; Aksay et al., 2007; Machens and Brody, 2008).

Assumption 1 (Two opposing populations). We assume that the dynamics are controlled by two population variables, X_R and X_L, one for each side. These variables shall represent the effect one population has on the postsynaptic currents of neurons in the other population; a specific definition will follow below. We assume that in its normal working regime, these two “synaptic” population activities oppose each other, so that, as the activity in one population grows, that of the other drops. For simplicity we therefore assume that (1)β=XR+XL where β is a constant value.

Assumption 2 (Eye position). We furthermore assume that within this working regime (but not outside of it) the eye position can be read out as (2)θ=XR−XL where right range positions are defined as positive eye positions, and left range positions as negative.

Assumption 3 (Firing rates). We furthermore assume that each neuron's firing rate in the stationary state is determined by its excitatory input from the same population, scaled by a weight a_i, the inhibitory input from the opposite population, scaled by a weight c_i, and some external (e.g., vestibular) input, h_i, so that (3)rR,i=[aiXR−ciXL+hi]+ (4)rL,i=[aiXL−ciXR+hi]+ where i = 1, …, n indexes the neuron, and [·]₊ is a threshold-linear function (Figure A1A). For simplicity, we assume complete symmetry of the two systems, so that neurons in the two populations have exactly the same set of parameters values.

Now that we have fully determined the relationship between the integrator neural activities and the eye position, we can appropriately interpret the pharmacological inactivation results of (Aksay et al., 2007). After unilateral silencing of the oculomotor integrator, the animals manage to fixate the eyes on the contralateral half of the motor range, and the activities of the intact neurons remain stable within this range. In our framework, this suggests a range of stable points on the upper half of the activity axis X_R and X_L. These stable ranges provide an important constraint for the ILA and NP models, limiting the possible choices of λ_i that we can make in Equations (9–11).

We have clarified the dynamics along the axes required by the pharmacological inactivation results of Aksay et al. (2007), so that we can proceed to the final tuning of the network model.

Independent of the choice of λ_i, the proposed network connectivities feature a line attractor in the (X_R, X_L) space that agrees with the recruitment order. However, the desired dynamics outside of the attractor, e.g., the dynamics along the axes required by the pharmacological inactivation results, further constrain the tuning of the parameter λ_i. Given the low number of neurons and simplicity of the tuning, we manually tuned this parameter to change the dynamics around the line attractor. However, a least squares optimization procedure could also be built to solve the problem.

In Figure A4, we present two model solutions presented in the main section, the independent line attractor model (ILA model) and the null position model (NP model) (see Appendix for the detailed parameter tuning procedure). In general, the nullcline of each half OI is the set of points in space where the dynamics in the corresponding direction are zero. In both ILA and NP models, we tuned the nullclines of either side to consist of densely arranged parallel lines, such that the vectors in the flow field are bound to be either horizontal or vertical. For instance, in the NP model, left range, vertical nullclines of the right side (in blue) ensure that the system has weak dynamics on the horizontal direction, and therefore the dynamics are mostly on the vertical direction.

We note that in the NP model, we have tuned the nullclines to end close to the line of fixed points, to ensure that each piece of nullcline intersects only one piece of the opposite nullcline and therefore obtain a line of fixed points (Figure A4B, left). If we relax that assumption and allow for a band attractor for instance, then the nullclines do not have to end as close to the diagonal.

The ILA model with 36 neurons on each half of the integrator is illustrated in Figure A4A. As can be seen, the nullclines resulting from the summation of the input-output functions H(.) of the individual neurons intersect in a line of stable points (Figure A4A, left). Furthermore, after complete left inactivation, the system is able to maintain stable positions in the right range (results not shown), as obtained in the pharmacological experiments (Aksay et al., 2007). Finally, this model shows the experimentally observed recruitment order feature, i.e., the higher the threshold the higher the slope of the tuning curves (Figure A2), since it was imposed into the model. The dynamics of the ILA model are qualitatively similar to the dynamics of the model developed in Aksay et al. (2007) (see Figure A4A, right). The respective network parameters are listed in Table A1.

In Figure A4B, we illustrate the NP model. Just as the ILA model, the system is able to maintain stable eye positions in half the original range after unilateral silencing. The model differs from the ILA model in one important aspect, however. The upper half of the nullcline (Figure A4B, left) is now maintained largely by self-excitation. In fact, the two halves of the oculomotor integrator are almost independent, and the threshold of the neurons is mostly determined by the self-excitatory inputs. However, due to the weak, mutual cross-inhibition, this independence is slightly disturbed, leading to the dynamics shown in Figure A4B, right. Altogether, the NP model has qualitatively the same dynamics below the line attractor as the ILA model but different dynamics above the line attractor. However, the NP model dynamics are slower than the ILA model dynamics in the whole state space, given the low inter-dependence of the two areas in the NP model (the arrow lengths in Figure A4A, right panel, and Figure A4B, right panel, are in different scales). The respective network parameters are listed in Table A1.

Despite the connectivity differences between ILA model and NP model, since these models share the same qualitative dynamics below the line attractor, NpHR left stimulation leads to very similar results in both models, with minor differences (results not shown).

Due to the use of Heaviside input-output functions and the resulting staircase character of the nullclines, the network implementation of the NP model features additional fine-scale dynamics that are not illustrated in Figures 7A, A4B (see detail of the dynamics in Figure A4C). These dynamics do not affect the trajectories in the noise-less case, and disappear when more realistic (e.g., sigmoidal) input-output functions are used for tuning the model (simulations not shown).

As shown in the results section, we experimentally test the range of models by performing instantaneous perturbations with optogenetic tools: inactivation with halorhodopsin (NpHR) and excitation with channelrhodopsin (ChR2). To model the perturbations, it is crucial to understand the effects of optogenetic tools in the activities of stimulated cells.

It has previously been shown (Arrenberg et al., 2009) that inactivation by NpHR stimulation leads to stronger inactivation of high activity cells than low activity cells, so the inactivation effects are approximately divisive. ChR2 stimulation, on the other hand, causes an increase in activity independent of the original activities (Arrenberg et al., 2009). To incorporate optogenetic stimulations in our dynamical systems model, we take these experimental observations into account.

A Matlab package in which the ILA model and the NP model are implemented is provided in the online supplementary information.

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