We used two-photon microscopy to image larval zebrafish with panneuronal expression of the genetically encoded calcium indicator GCaMP5G (Experimental Procedures and (Akerboom et al., 2012)). Larvae were restrained in agarose, and the agarose removed from the eyes and tail to allow for unimpeded movement. To elicit the OKR, we projected sinusoidally rotating, whole-field visual stimuli on a screen below the fish (Figure 1A), which reliably evoked conjugate horizontal eye movements (Beck et al., 2004; Huang and Neuhauss, 2008). These optokinetic responses consisted mainly of slow tracking movements phase-locked to the stimulus and were accompanied by occasional fast reset saccades and swims (Figure 1B and Figure S1A; Movie S1). Eye and tail movements in response to visual stimulation were recorded using an infrared-sensitive camera (Figure 1A). Since the tracking responses persisted robustly for many hours (Figure S1B), we were able to systematically image, under comparable behavioral conditions, the majority of the brain (Figure 1C) of individual fish with a voxel spacing of less than one micron in all three dimensions.

Stimulus-locked, phasic activity was evident in areas throughout the brain, including both cell somata and neuropil regions (Movie S2). To analyze the activity dynamics across different brain areas, we designed an automated algorithm to segment the image into regions of interest (ROIs) by assigning voxels to ROIs based on their correlations with neighboring voxels (see Supplemental Experimental Procedures). Since, for most ROIs, fluorescence changes occurred at the stimulus frequency (Movie S2), we first focused our analysis on the temporal relationship between ROI activity and stimulus. To visualize the whole-brain activity dynamics, we assigned each ROI a hue based on the timing of its response relative to the stimulus (Figure S2), with intensity representing the relative power at the stimulus frequency (dorsal view Figure 1D, other views Movie S3). We observed that the neural activity pattern showed a complex, but markedly ordered, distribution, which is repeated symmetrically on both sides of the brain, with corresponding regions on opposite sides active in antiphase (Figure 1D). The consistent structure of this pattern is evident when the map is flipped left to right, and phase shifted by 180 degrees (Figure 1D, bottom).

Activity locked to the stimulus was found throughout the brain, particularly in the optic tectum and pretectum, cerebellum, habenulae and throughout the caudal brainstem. Dividing up the ROIs by brain area (Movie S4), we found roughly equal numbers on the left and right in each region. A notable exception was the dorsal habenula, a region known to have pronounced anatomical asymmetry (Bianco and Wilson, 2009), where we found responses predominantly on the left side (Figure 1D,E). The activity in the habenula, as well as in the torus longitudinalis, shows a distinct temporal pattern from the rest of the network, and may relate to the light/dark structure of the projected stripe pattern, rather than the visual motion signals that drive the OKR.

Many ROIs align with individual cell somas (Figure 1D, panel i), while some are localized to synaptic regions, such as thin stripes in the superficial layers of the tectum (Figure 1D, arrowhead in panel ii), which are likely to be the axon terminals of direction-selective retinal ganglion cells (Gabriel et al., 2012; Nikolaou et al., 2012). Large domains of activity in neuropil are seen in the pretectal/thalamic retinal arborization areas (Figure 1D, panel i and iii), and two areas in the ventral hindbrain (Figure 1D, panel iv). Visual input from the retina reaches the brain through 10 retinal arborization fields (AFs) (Burrill and Easter, 1994); however, it is not known which of these AFs are part of the OKR circuit. AF6 has been proposed to correspond to the mammalian nucleus of the optic tract (Roeser and Baier, 2003), but the dark half mutant, which lacks AF6 and other ventral tract AFs has an intact OKR (Muto et al., 2005). We observed a large domain of activity spanning from AF4 to AF6 along the rostrocaudal axis, and extending dorsally to the most rostral portion of AF9, suggesting the possibility that multiple redundant pathways exist to carry visual information for the OKR. In the dorsal medulla, we observed that activity has a columnar structure (Figure 1D, panel v), with interleaved regions active at opposite stimulus phases, which is interesting in light of the fact that this region is known to be organized into stripes of neurons with different projection patterns and neurotransmitter types (Kinkhabwala et al., 2011). It would be interesting to determine if the functional segregation we observe coincides with the anatomical organization described.

To further characterize the response dynamics, we looked at the distribution of phases at which the activity peaked for all active voxels (Figure 1F). This distribution shows four major peaks, with responses on both sides of the brain occurring mainly during ipsilateral nasal-to-temporal (ipsiversive) stimulus motion (Figure 1F). Activity was clustered in a narrow range of stimulus phases encompassing the time between maximum stimulus velocity and maximum eye position (Figure 1F). Responses in areas associated with sensory processing, such as the tectal neuropil, pretectum and inferior olive tended to precede activity in motor structures in the hindbrain (Figure 1G). Most areas, including the pretectum and inferior olive, contained regions active in phase with ipsiversive rotation, although the cerebellum showed the opposite pattern, consistent with crossed projections from the inferior olive (Straka, 2006). The activity observed in the pretectum is consistent with recordings from the homologous region in mammals, the NOT, that show selectivity for ipsiversive movement (Collewijn, 1975a).

Clustering of the raw traces, based on their temporal profiles, confirms the division of activity into four laterally symmetric, early and late clusters (Figure 2A). Four clusters of activity with similar profiles were found in every fish with a robust OKR (Figure S3). Based on their anatomical localization (Figure 2B–C) and their temporal order (Figure 2A), these clusters likely represent a division of the circuit into distinct sensory and motor domains, with the two structures repeated bilaterally.

In order to identify regions that were consistently activated during the OKR, and to quantify the stereotypy of the spatiotemporal pattern of activity, we compared the whole-brain activity maps across multiple fish. To that end, data from each fish was registered in three dimensions to a reference brain (Figure 3A) (Rohlfing and Maurer, 2003; Jefferis et al., 2007) with a registration accuracy on the order of one cell soma diameter (4.5 ± 0.5 μm; see Experimental Procedures). Movie S5 shows how the ROI maps across two or more fish were aligned and averaged to reveal a core set of three-dimensional structures that were reliably activated. We quantified the degree of similarity across fish first by calculating the probability of finding an ROI in each registered voxel. Many voxels showed activity in nearly every fish tested (Figure 3B). To assess the stereotypy of the temporal structure of the responses, we asked how far, on average, you have to travel from any active region to find a region with similar activity in another fish (see Experimental Procedures). For most ROIs, this distance was less than 15 μm, and, in an extensive network throughout the brain, especially in the hindbrain and neuropil regions around the retinal arborization fields, it was less than 2 μm (Figure 3C,E), which is on the scale of a single cell soma in the larval fish brain, and at the limit of our registration accuracy. The stereotypy and bilateral symmetry are particularly evident when we compare each fish’s brain to a left-right flipped version of the same brain, allowing for a 180 degree phase shift (Figure 1D, bottom panel). In this case, most ROIs are less than 10 μm away from a matching ROI on the opposite side (Figure 3D,E). Some cells in the left habenula still find matches, on a slightly longer spatial scale, due to small numbers of similarly responding cells in the right habenula. In the stratum periventriculare (spv) of the optic tectum, where a few active neurons are scattered within a dense layer of cell bodies, the precise location of the active somas was more variable from individual to individual (Movie S5).

Taking advantage of this stereotypy, we averaged the raw pixel fluorescence data across fish to reveal structure in weakly active areas that lie below the detection threshold for individual fish. Combining data from three fish with closely matched behavior we detected activity in 65% more voxels than the combination of the three individuals alone. In many places, the alignment is sufficiently precise to allow the identification of fiber tracts, based on functional signals, between active cell bodies and co-active neuropil regions (arrowhead, Movie S6). The resulting average activity map reveals a distributed network, consisting of fewer than 5% of all imaged voxels, but which nevertheless spans a large volume of the brain (Figure 3F, Movie S6). As expected, activity was evident in regions such as the oculomotor nucleus and multiple retinal projection fields, confirming the validity of our activity maps (Figure 4, Movie S7). Surprising features are also revealed, such as left-right alternating activity in neuropil in the vicinity of the interpeduncular nucleus or the closely associated median raphe (Figure 4), both of which receive bilateral innervation from the habenula (Aizawa et al., 2005; Amo et al., 2010).

While activity in each region is concentrated around particular response phases, within most areas we observed smooth spatial gradients of activity timing, which propagate either along different anatomical axes, or radiate from the center of active structures outwards (Figure 5, Movie S8). Prominent examples included one group of pretectal neurons in which responses progress from rostral to caudal, while a distinct rostral sweep is seen in the cerebellar neuropil. Several such gradients are evident in structures throughout the hindbrain. These could result from gradients of physiological or tuning properties, as have been described in neurons in both the oculomotor integrator and spinal motor networks in zebrafish (Mclean et al., 2007; Miri et al., 2011a), or from topography of local or inter-region connectivity. These gradients could serve several functional purposes. Some may reflect signal transformations, such as integration of activity within areas. Smooth progressions of responses times could also permit specific sequencing of actions, or provide a substrate for learning properly timed and coordinated responses (Mauk and Buonomano, 2004).

We hypothesized that the early and late response peaks in Figure 1F reflected processing of sensory and motor signals. However, when presented with motion throughout the visual field, fish moved both eyes together with similar amplitude and phase relative to the stimulus (Figures 1B and S1A), making it hard to distinguish contributions of individual signals to the observed neural activity. To dissociate left and right eye, and sensory and motor signals in different areas, we presented stimuli on the left or right side of the fish only, which causes different amplitude of motion of the two eyes, as well as in-phase and out-of-phase stimulation to both eyes which approximates translational motion and results in either conjugate eye rotation or vergence eye movements and optomotor swimming (Figures 6A and S1C–D, Movie S9). In order to identify the sensorimotor signals represented in the activity of the ROIs, we correlated the measured fluorescence of every ROI with the predicted fluorescence time course arising from an array of different stimulus and behavioral parameters (Figure 6B and Figure S4) (Miri et al., 2011b).

Many active regions showed strong correlations with behaviorally relevant variables (Figure 6C and S6). Some could be classified as carrying sensory signals relating to visual motion, such as stimulus velocity (Figures 6C,D and S6), while others carried sensory information that was not directly related to the OKR. For example, many cells in the left dorsal habenula and torus longitudinalis (TL) were strongly modulated by the presence of the visual stimulus on one side or the other, probably encoding differential light levels between the two eyes (Figures 6C,D and S5) (Northmore et al., 1983). Other regions, which we classified as motor, were correlated with eye and tail variables (Figure 6C,D). For instance, activity in ROIs clustered near the dorsal surface in the caudal hindbrain was strongly correlated with nasal-temporal positions of the ipsilateral eye (Figures 6C, S5 and S6), which is consistent with responses found in Area I cells in fish, which form part of the horizontal velocity-to-position integrator (Pastor et al., 1994). Other regions were correlated with more complex combinations of sensorimotor variables. Responses of many tectal neurons were inhibited in a direction-selective manner by stimuli presented to the ipsilateral side (Figures 6C,D, S5 and S6). They responded much more strongly when motion was presented in opposite rotational senses to both eyes, which is consistent with translational movement through the water (Frost et al., 1990), than when whole-field rotational motion was shown. Thus, these neurons are unlikely to be directly involved in the OKR, which is in agreement with previous ablation studies (Roeser and Baier, 2003). Conversely, neurons in some hindbrain areas were selective for whole-field rotating stimuli and unresponsive to whole-field translation motion stimuli (Figures 6C, S5 and S6), suggesting that these neurons may be involved in the OKR.

We next analyzed the three-dimensional arrangement of the voxels correlated with these sensorimotor variables, in order to identify different functional networks across the whole brain. To identify these structures in an unbiased way, we took fluorescence traces from overlapping cubes of approximately 5 μm, tiling the entire imaging volume, correlated each of these cubes with 191 different sensorimotor variables, and then aligned these correlation maps to the reference brain. Figure 7 shows the distributions of cubes whose activity best correlated with different parameters, summed across seven fish (Figure S7). Each signal localizes to one or a few distinct areas, which are similar between the left and right sides, and align very well across different individuals (Figure 7). One exception to this pattern is the response to dimming on one side during hemifield stimulation. In this stimulus condition, cells in the left habenula are selective for both left and right side stimulation, while there is little activity in the right habenula. In addition, there is a striking asymmetry in the ventral cerebellar area of the hindbrain, where a small area on the left is responsive during left stimulus presentation, whereas a more diffuse bilateral activation is observed during presentation on the right (Figure 7), suggesting different processing of information from the two eyes.

Individual regressors may not fully capture the functional properties of a given neuron. We hypothesized that a more accurate functional classification could result from analyzing the correlations with all relevant regressors, rather than the best individual fit. We therefore performed a k-means clustering analysis on the complete vectors of correlation coefficients for voxels from seven fish brains. Several symmetric pairs of functional clusters were found, some of which match very well to the single parameter maps in Figure 7 (Figure 8A–C, S8 and Movie S10). Others were not so clearly linked to a single variable, but were broadly correlated with, for example, the left and right eye velocity and position signals that were consistent with clockwise or counterclockwise conjugate rotation. In addition, clear and consistent symmetrical structures associated with swimming activity could be extracted from individual fish, even though the average activity trace from this cluster was rather flat, since this behavior was poorly correlated across individuals, or across stimulus repetitions, under these conditions (Figure 8B,C). Some areas correlated with swimming also show strong correlations with other parameters, such as eye acceleration due to rapid eye movements that occur when swimming, and therefore may not all be grouped together based only on peak correlation value. The multidimensional clustering therefore gives a much more comprehensive map of this hindbrain locomotor network (compare Figure 7, 8B and S8).

In our experiments using a sinusoidally-modulated OKR stimulus, a large cluster of neurons in the pretectum showed responses in phase with temporal to nasal (TN) motion in the contralateral eye (e.g. Figures 1D(i,ii),1G , 4, Movie S6). Experiments in many different species have shown that sensory inputs to the horizontal OKR originate from ON-type, direction-selective ganglion cells, which project to pretectal and accessory optic nuclei, in particular the nucleus of the optic tract and dorsal terminal nucleus (NOT-DTN) in mammals, and the nucleus lentiformis mesencephali (LM) in other tetrapods (Oyster et al., 1972; Collewijn, 1975b; Soodak and Simpson, 1988; Büttner-Ennever and Horn, 1997; Masseck and Hoffmann, 2009). Similar to our observations, neurons in both these areas respond mostly to TN motion in the contralateral eye. In fish, two candidate nuclei were originally suggested as homologs of the LM and the accessory optic nucleus of the basal optic root (nBOR) (Finger and Karten, 1978), based on features of their projection patterns, but, more recently, others have proposed a single pretectal nucleus in fish, which contains neurons responding to all directions of motion (Klar and Hoffmann, 2002; Masseck and Hoffmann, 2009). The pronounced asymmetry in directional sensitivity that we find in our data may account for the observation that, in larval zebrafish, as in many species, TN motion drives stronger OKR than nasal to temporal (NT) motion (Roeser and Baier, 2003; Burgess et al., 2009), but does raise the question of how monocularly presented OKR stimuli can drive a response in both directions. Burgess et al. observed that robo3 mutant zebrafish, which have impaired formation of commissures, have a specific deficit for NT motion in the monocular OKR, and propose that NT motion acts by suppressing crossed inhibitory connections between pretectal nuclei (Burgess et al., 2009; Prochnow et al., 2007). Interestingly, here we find functional evidence for a crossed inhibitory pathway with a preference for TN motion, in the responses of translation-sensitive neurons in the optic tectum (Figures 6C and 8). These respond in phase with contralateral NT motion, except when it occurs simultaneously with ipsilateral TN motion.

Responses to whole-field motion stimuli have previously been found in retinal ganglion cell terminals and neurons in the optic tectum of zebrafish larvae (Niell and Smith, 2005; Del Bene et al., 2010; Nikolaou et al., 2012). Similarly, we found signals in distinct, superficial layers of the tectum responsive to each of the two directions of motion (e.g. Figures 1Dii and 4), and, although most neurons in the large periventricular layer did not show stimulus-locked responses, a substantial scattered population did, and could, in principle, be involved in driving the OKR. However, ablation studies in both zebrafish larvae and adult goldfish have indicated that the optic tectum is not required for this behavior (Springer et al., 1977; Roeser and Baier, 2003). It is striking that up to half of the tectal neurons that responded to visual stimuli in our experiments appear to integrate information from both eyes (Figure 6C–D, 8), since their response to binocularly presented motion differed from the response to contralateral motion alone. In particular, responses to rotational motion stimuli were strongly suppressed. While we are not aware of previous functional evidence for such neurons integrating binocular visual information in the tectum, anatomical tracing experiments have identified bilateral projections in mammals from the pretectum to the superior colliculus, the optic tectum homolog, and in vitro recordings have confirmed functional inhibitory connections (Baldauf et al., 2003; Born and Schmidt, 2004). The preference for NT motion suggests that they would respond strongly to re-afferent motion signals during forward locomotion, and could be involved in gating visual processing, or detecting the proximity of objects, during swimming. This scattered population represents a very small percentage of the neurons in the tectum, which, in the zebrafish larva, contains on the order of 10,000 cells, so these cells would only be rarely sampled in electrophysiological recordings, highlighting the potential of whole-brain imaging methods such as the one used in this study.

In many species, retinal motion signals pass through two integration processes which generate both a stored representation of eye velocity (velocity storage mechanism, VSM, Figure 9), which can serve to make the OKR robust to interruptions of sensory input, and a velocity to position neural integrator (VPNI), which generates eye position signals that are used to counteract elastic forces which pull the eye towards its resting position, and to trigger fast reset saccades during optokinetic nystagmus (OKN). Although early experiments argued against the existence of velocity storage in larval zebrafish (Beck et al., 2004), more recently a VSM has been shown to exist in fish as early as five days post fertilization, but its presence was masked in previous experiments by the rapid decay of the position integrator at this stage of development (Chen et al., 2013a). The VPNI has been investigated in larval zebrafish in the context of maintenance of eye position following spontaneous saccadic movements in the dark (Miri et al., 2011a; 2011b). Curiously, the distribution of neurons we observe with the strongest eye position correlations are distributed differently from those previously observed, occupying a more restricted volume at a more dorsal level. Nevertheless, the region we describe does still fall within the volume used for optogenetic localization of the VPNI (Miri et al., 2011a). This apparent discrepancy could be explained by the fact that the two experimental paradigms differ in both motor dynamics (saccades vs. smooth movements) and sensory inputs (continuously varying motion vs. darkness), and it will be interesting to investigate the activity across these populations during different behaviors.

The OKR also shows learned adaptation driven by visual feedback, for which the cerebellum has been shown to be important (Nagao, 1983). The cerebellum receives OKR-related sensory and motor signals through mossy fibers, and in particular, retinal image motion signals via climbing fibers from the contralateral inferior olive, which in turn receive direct input from the pretectum. We observe a notable temporal gradient of activity timing in this olivocerebellar pathway during OKR (Figure 5 and Movie S8), which could arise from spatial gradients of dynamics or velocity tuning, and may provide a substrate for fine-tuning of the dynamics of adaptation. Although adaptation of the OKR has not been observed in larval zebrafish, rapid adaptation of the optomotor response, which depends on both cerebellar function and on visual motion signals from the inferior olive, has been shown as early as six days post fertilization, suggesting that this circuit is functional even at early developmental stages (Portugues and Engert, 2011; Ahrens et al., 2012)

Current attempts to map functional organization, or connectivity, of the brains of entire animals at cellular resolution, raise the question of how similar such maps are across individuals. The correct approach to, and success of, proposed large-scale projects to assemble brain activity maps (Alivisatos et al., 2012) depends critically on the answer to this type of question. Previous studies in fruitflies have revealed anatomical stereotypy on the order of a few microns (Jefferis et al., 2007; Datta et al., 2008), but studies of neuromuscular connectivity in vertebrates found substantial variability (Lu et al., 2009). Meanwhile, in humans, spatiotemporal activity patterns from fMRI data have revealed similarities in processing of visual stimuli across individuals (Hasson, 2004). In this study we have compared, across individual animals, neural circuit activity in the whole brain during behavior, with micrometer resolution. We find that a large part of the functional network active during the OKR is extremely stereotyped, such that it is possible to reliably predict the activity in individual voxels from their anatomical position alone. This has important practical consequences, since these areas can now be targeted, without further functional mapping, for ablations, electrophysiological recordings, optogenetics or photoactivation, and also aligned with the scaffolds for existing anatomical atlases to find molecular markers for specific loci of activity (Ronneberger et al., 2012).

While two-photon laser scanning microscopy offers many advantages for rapid, long-term volumetric functional imaging, especially during visual behavior (Renninger and Orger, 2013), alternative methodologies are advancing which offer different trade-offs. Two recent studies in larval zebrafish used light-sheet imaging to image simultaneously from multiple z-planes (Panier et al., 2013; Ahrens et al., 2013). Scanned light-sheet imaging allows for the collection of many more voxels per unit time than standard two-photon microscopy, since an entire line, rather than a single point, is being exposed at any point in time. A disadvantage of the current plane-illumination methods, however, is that they require a constantly moving, visible laser passing laterally through the brain. This light may be sensed by the eye, which makes the study of visual responses more challenging, but could also produce responses directly through activation of light-sensitive areas in the brain (Fernandes et al., 2012; Kokel et al., 2013). Methods which combine two-photon excitation and plane illumination are a promising direction for the future. For example, wide-field two photon excitation, through temporal focusing, was recently used to image simultaneously from much of the brain of C. elegans (Schrödel et al., 2013).

Whole-brain, single-cell resolution functional maps as described here, especially when combined with connectivity maps and high-resolution anatomical atlases (Ronneberger et al., 2012; Arrenberg and Driever, 2013) have the power to reveal new functionally defined structures in the brain and provide a fundamental framework for understanding the complete neural circuits for behavior (Alivisatos et al., 2012).

Adult fish were maintained at 28°C on a 14:10 hour light cycle following standard methods. Larvae were raised in E3 embryo medium (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl₂ and 0.33 mM MgSO₄). All procedures followed NIH and the Portuguese Direcção Geral de Veterinária guidelines and were approved by the Standing Committees on the Use of Animals in Research of Harvard University and Champalimaud Foundation.

Five and six day-old nacre (mitfa −/−) elavl3:GCaMP5G zebrafish larvae were placed in a drop of 1.5% low-melting-temperature agarose in a Petri dish on a Sylgard 184 (Dow Corning) base and immersed in E3 solution. The agarose around the tail, caudal to the pectoral fins, and eyes was cut away with a fine scalpel to allow for movement. The dish was placed onto a light-diffusing screen and imaged on a custom-built two-photon microscope. A Ti-Sapphire laser (Spectra Physics Mai Tai) tuned to 950 nm was used for excitation. Larval brains were systematically imaged, while they performed the OKR. Frames were acquired at 2 Hz in four, 0.83 μm spaced interlaced scans allowing each area of the brain to be evenly sampled at around 8 Hz. Visual stimuli, generated using LabView, were projected at 60 frames per second using a 3M MPro110 microprojector, and a red long-pass filter (Kodak Wratten No.25) to allow for simultaneous imaging and visual stimulation. The stimulus consisted of radial light and dark stripes with a period of 45 degrees (Movie S1). The stimulus position was rotated sinusoidally, centered on the fish, at 0.1 Hz with amplitude of 40 degrees. The focal plane was shifted by 0.5 μm every 3 stimulus repetitions. The fish were illuminated from above using a ring of infrared LEDs (850 nm wavelength) coupled to the objective. To track eye and tail movements, a small hole was cut in the diffusing screen to allow the fish to be imaged from below at up to 200 frames per second using an infrared-sensitive CCD camera (Pike F032B, Allied Vision Technologies). Tail and eye image data was streamed to hard disk using software written using National Instruments LabView.

Eye angles and tail movements were tracked using custom routines in Matlab (Mathworks, USA). For details see Supplemental Experimental Procedures.

Image analysis was performed with Matlab (Mathworks, USA). To correct for drift or small movements of the fish, each image frame was aligned, using translation only, to the average image of that z-plane, and consecutive z-planes were aligned to each other with subpixel precision (Nestares and Heeger, 2000). Occasionally vigorous swimming movements caused a large motion in a single frame. Such isolated frames that could not be aligned with a correction of less than 2 pixels were removed from analysis. Region of interest (ROI) segmentation was performed using automated algorithms, based on local correlations between pixels (see Supplemental Experimental Procedures).

Image registration was performed using the free Computational Morphometry Toolkit (CMTK - http://www.nitrc.org/projects/cmtk/) (Rohlfing and Maurer, 2003). For each fish, an anatomical stack was computed by summing the GCaMP5G fluorescence. One of the these stacks was chosen as the reference brain, and non-affine volume transformations were computed to align each fish’s anatomical stack to this reference stack using the affine and warp functions. These transformations were then used to transform individual ROIs from each fish into the frame of reference of the reference brain, allowing us to compare the anatomical location of ROIs from different fish. The precision of the registration was evaluated by selecting salient anatomical features, identifying them in every fish and checking how closely they were transformed to the equivalent feature in the reference brain. Before transforming, the average distance between the features in individual and reference brain was 48.3 ± 4.7 μm (mean ± SEM). After transforming the distance between the features in the transformed brains and the reference brain was 4.5 ± 0.5 μm.

The phase of a particular ROI’s activity relative to the stimulus was determined by taking the Fourier transform of the z-score of its fluorescence trace, and determining the phase of the 0.1 Hz component. To make color-coded images, each ROI was assigned a hue, based on its phase, from the Matlab HSV (Hue Saturation Value) colormap, and a value proportional to the amplitude of the 0.1 Hz component, reflecting how much modulation in that ROI was at the stimulus frequency (Figure S2). Different anatomical regions were defined manually in the reference brain (Movie S4). For reasons of symmetry, when analyzing the phase distribution in particular anatomical regions, a phase shift of 180° was applied to all the ROIs from the right brain, before pooling with the ROIs from the left brain. The laterality index was calculated using the formula (voxL-voxR)/(voxL+voxR) where voxL and voxR are the number of voxels found in ROIs detected on the left and right of the brain respectively. Cluster analysis of the traces was performed on individual fish using the k-means method. The optimal number of clusters was selected by inspection of silhouette plots for different even numbers of clusters (Figure S3). Reliability of finding activity in different areas (Figure 3B) was measured by the percentage of fish, following alignment, with ROIs including each voxel. The spatial scale of stereotypy was measured by asking for every voxel in every ROI in every fish, how far you have to go, in all fish, to find a voxel that is active within 18° of the same stimulus phase (Figure 3C). To assess symmetry, we asked how far you need to go in a reference-aligned, left/right-flipped version of the same brain, to find an ROI within 18° of the opposite stimulus phase. To control for the overall spatial distribution of ROIs, the same analysis was performed with the phases of the starting ROIs randomly shuffled. Three-dimensional renderings were performed using the ImageJ 3D Viewer plugin.

Regressors for correlation analysis were constructed from a set of 191 behavioral and stimulus-related variables (e.g. position and velocity of the eyes, stimulus velocity, onset/offset of stimuli). These were convolved with a kernel with an exponential decay based on the measured half-decay time for GCaMP5G (0.667s) (Chen et al., 2013), to produce a set of predicted fluorescence traces (Miri et al., 2011b) (Figure 6B and Figure S4), and were compared with the measured fluorescence traces by correlation. The interlaced scan waveform meant that each frame was scanned in four, slightly offset frames. Therefore, any ROI which spanned at least four adjacent scanned lines was partially imaged in each one of these four frames. Raw fluorescence traces were deinterlaced, by considering the voxels in each of these four subframes independently, to yield fluorescence traces at 8 Hz and smoothed with a 1.125s width boxcar filter prior to correlation. For Figure 6C, correlation analysis was performed for automatically-segmented ROIs to identify the variable that best described the signal of individual ROIs. For Figures 6D, 7 and 8, systematic correlation analysis was performed on cuboidal ROIs of 8 × 8 × 9 voxels (roughly 6.6 × 6.6 × 4.5 microns), tiling the brain in steps of 2 (x,y) or 3 (z) voxels. To assign voxels to a particular functional group, we required a threshold of at least 0.3 for the best absolute correlation coefficient. Based on comparisons with shuffled sets of regressors this level of correlation rarely arises by chance, and is therefore a conservative choice for detecting voxels with behavioral correlation (Figure S7A). At the same time the precise threshold chosen is not critical to the set of areas identified (Figure S7B). Since the maximum r value depends on the original choice of regressors, it remains possible that new functional groups may be identified through better identification of appropriate regressors, as well as by indicators, or fish lines, that offer better signal to noise in particular cell types. For Figure 8, clusters of correlation patterns were found as follows: First the set of correlations with all behavioral regressors was calculated for all cuboidal, tiled ROIs across 7 fish brains, and all those in which one regressor had a correlation of at least 0.3 were pooled (1,283,831 ROIs). Principal component analysis was applied to the data set, and k-means clustering was performed on the data projected onto the first 10 components. Different numbers of clusters were tested, and the number of 25 clusters was selected manually because it gave the clearest separation of functional classes. Clustering with larger numbers produced inconsistent fracturing of anatomical structures or segregation of clusters across fish, suggesting that true functional classes were being artificially divided.