Predicting modular functions and neural coding of behavior from a synaptic wiring diagram

Vishwanathan, Ashwin; Sood, Alex; Wu, Jingpeng; Ramirez, Alexandro D.; Yang, Runzhe; Kemnitz, Nico; Ih, Dodam; Turner, Nicholas; Lee, Kisuk; Tartavull, Ignacio; Silversmith, William M.; Jordan, Chris S.; David, Celia; Bland, Doug; Sterling, Amy; Seung, H. Sebastian; Goldman, Mark S.; Aksay, Emre R. F.

doi:10.1038/s41593-024-01784-3

Article
Open access
Published: 22 November 2024

Predicting modular functions and neural coding of behavior from a synaptic wiring diagram

Ashwin Vishwanathan,
Alex Sood,
Jingpeng Wu,
Alexandro D. Ramirez,
Runzhe Yang^,,
Nico Kemnitz,
Dodam Ih,
Nicholas Turner^,,
Kisuk Lee^,,
Ignacio Tartavull,
William M. Silversmith,
Chris S. Jordan,
Celia David,
Doug Bland,
Amy Sterling,
H. Sebastian Seung^,,
Mark S. Goldman^,,,
Emre R. F. Aksay &
the Eyewirers

Nature Neuroscience volume 27, pages 2443–2454 (2024)Cite this article

9956 Accesses
2 Citations
144 Altmetric
Metrics details

Abstract

A long-standing goal in neuroscience is to understand how a circuit’s form influences its function. Here, we reconstruct and analyze a synaptic wiring diagram of the larval zebrafish brainstem to predict key functional properties and validate them through comparison with physiological data. We identify modules of strongly connected neurons that turn out to be specialized for different behavioral functions, the control of eye and body movements. The eye movement module is further organized into two three-block cycles that support the positive feedback long hypothesized to underlie low-dimensional attractor dynamics in oculomotor control. We construct a neural network model based directly on the reconstructed wiring diagram that makes predictions for the cellular-resolution coding of eye position and neural dynamics. These predictions are verified statistically with calcium imaging-based neural activity recordings. This work demonstrates how connectome-based brain modeling can reveal previously unknown anatomical structure in a neural circuit and provide insights linking network form to function.

Neural signal propagation atlas of Caenorhabditis elegans

Article Open access 01 November 2023

Connectome-constrained networks predict neural activity across the fly visual system

Article Open access 11 September 2024

Neuronal wiring diagram of an adult brain

Article Open access 02 October 2024

Main

A fundamental question is how the neural coding properties of a circuit relate to its connectivity structure. Recent connectomic efforts have enabled the generation of wiring diagrams in the fly and in vertebrate circuits like the retina^1,2, utricular system^3,4, olfactory bulb⁵, visual thalamus⁶, visual cortex^7,8 and somatosensory cortex⁹. In the fly navigation system, models for the organization of the head direction system were literally realized in the ring-like structure of the central complex¹⁰, and the topographically organized sinusoidally varying synaptic connectivity weights directly suggested the trigonometric identities required to encode coordinate transformations^11,12. Likewise, in the early visual system of both invertebrates¹³ and vertebrates¹⁴, connectomics has provided a strong anatomical basis for the computations underlying motion processing. By contrast, within and across most vertebrate brain regions, the connectivity is highly recurrent without obvious topographic structure. In such cases, it remains unclear how the neural coding properties of the system emerge from and relate to the underlying connectivity and cellular properties.

We address this question within a multiregion network of the vertebrate hindbrain responsible for maintaining an animal’s gaze. This network contains the oculomotor neural integrator, which mathematically accumulates (that is, integrates in the sense of calculus) transient eye velocity-encoding command signals into persistent eye position-encoding commands. Persistent activity in this system is thought to be realized by low-dimensional attractor dynamics^15,16. Low-dimensional attractors have been proposed to underlie a large range of neural computations, from storage of working memory to reduction of noise in sensory and cognitive representations^16,17. Although it is generally understood that recurrent interactions in a network play an important role in the generation of attractor dynamics, we lack understanding of what features of the recurrent circuitry lead to the observed low-dimensional patterns of activity. Within the oculomotor integrator, there is broad physiological^18,19,20 and coarse anatomical²¹ evidence implicating recurrent excitatory interactions in generating attractor dynamics. However, even in this relatively simple setting, we do not yet have a wiring diagram for the precise patterns by which neurons are recurrently interconnected, a necessary substrate for any quantitative, predictive model of network dynamics. To address this limitation, we reconstruct a synapse-resolution wiring diagram from a vertebrate brainstem region that includes the oculomotor neural integrator and more broadly serves a broad range of visual and sensorimotor functions. The wiring diagram is based on a previously acquired three-dimensional (3D) electron microscopy (EM) image of a larval zebrafish brainstem volume that is roughly bounded by a 250 × 120 × 80 μm³ box²². We reconstructed ~3,000 neurons as completely as the borders of the volume allow through human proofreading of an automated segmentation.

Classical anatomical studies have suggested that brainstem reticular circuitry within which the studied regions lie is ‘undifferentiated’²³ or ‘diffuse’²⁴. Furthermore, this region contains a tightly interwoven mixture of neurons involved in eye movements and, separately, body movements controlled by the spinal cord²⁵. By applying graph clustering algorithms to the connectome, we reveal that this naively disorganized structure hides a strongly modular hierarchical organization. At the highest level, we find two anatomically defined modules that are specialized for different behaviors: one for eye movements and the other for body movements. The oculomotor module is in turn subdivided into two submodules specialized for movements of each eye. Further, each submodule contains a three-block cyclic substructure. Our linking of structure and behavior by modularity analysis at synaptic resolution is unique for the vertebrate nervous system.

Using this wiring diagram, we generated a computational network model to predict neural coding properties at cellular resolution across three brain regions. When relating the wiring diagram to function, a common approach starts by extracting rules of synaptic connectivity between cell types and then applies these rules to explain function²⁶. This approach has worked well for direction selectivity in the retina^1,14,27 but is not guaranteed to generalize to other brain structures that seem less stereotyped and precisely organized than the retina. Here, we take a more direct approach: use the wiring diagram to estimate a synaptic weight matrix characterizing physiological interactions between neurons and literally insert that matrix into a network model incorporating a minimal number of additional constraints from physiology. Surprisingly, the predictions turn out to be statistically consistent at a population level with neural activity recorded by calcium imaging of larval zebrafish during oculomotor behavior.

Results

Neuronal wiring diagram reconstructed from a vertebrate brainstem

We applied serial section EM to reconstruct synaptic connections between neurons in a larval zebrafish brainstem (Fig. 1). The imaged volume includes neurons involved in eye movements^16,21,22,28. First, the dataset of Vishwanathan et al.²² was extended by additional imaging of the same serial sections. By improving the alignment of EM images relative to our original work, we created a 3D image stack that was amenable to semiautomated reconstruction²⁹ of 100× more neurons than before. We trained a convolutional neural network to detect neuronal boundaries (Methods)²⁹ and automatically segment the volume. To correct errors in segmentation, we repurposed Eyewire, which was originally developed for proofreading mouse retinal neurons¹⁴. Eyewirers proofread ~3,000 objects, which included neurons with cell bodies in the volume as well as ‘orphan’ neurites with no cell body in the volume (Fig. 1a). We will refer to all such objects as ‘neurons’ of the reconstructed network. Convolutional networks were used to automatically detect synaptic clefts and assign presynaptic and postsynaptic partner neurons to each cleft³⁰ (Fig. 1b–d and Methods). The reconstructed dataset contained 2,884 neurons, 44,949 connections between pairs of neurons and 75,163 synapses (additional relevant features are summarized in Supplementary Fig. 1). Most connections (65%) involved just one synapse, but some involved two (19%), three (7.9%), four (3.7%) or more (4.0%) synapses, up to a maximum of 21 (Supplementary Fig. 1a).

**Fig. 1: EM reconstructions of brainstem neurons.**

After registration of our EM volume to the Z-Brain reference atlas³¹ (Supplementary Fig. 2), several important groups of neurons and neurites were identified (Methods). Comparison with transgenic lines and cranial nerves in the atlas identified the abducens (ABD) neurons that control extraocular muscles, secondary vestibular (Ve²) neurons³² involved in sensorimotor transformations such as the vestibulo-ocular reflex (Fig. 1e–h, Supplementary Fig. 3 and Methods), and two populations of reticulospinal (RS) projection neurons involved in escape behaviors and controlling body movements^33,34,35. The neurons of one RS population were large and dorsally located, and the neurons of the other RS population were smaller and ventromedially located (Fig. 1h and Supplementary Fig. 4).

Axial and oculomotor modules in the brainstem

Complex systems, both biological and artificial³⁶, can often be divided into modules such that interactions within modules are stronger than interactions between them. Although such modules are structurally defined, they often turn out to be functionally specialized. To look for modularity in our data, we divided the reconstructed neurons into ‘center’ and ‘periphery’ (Methods). Neurons in the ‘center’ (540 neurons, 283 with soma) are recurrently connected to other reconstructed neurons. Neurons in the ‘periphery’ (2,344 neurons), by contrast, have negligible (62 neurons) or zero recurrent connectivity as quantified by eigenvector centrality (Methods) and are involved in feedforward pathways that supply input to the center or transmit output from the center; for example, the periphery includes ABD neurons and most RS projecting neurons. We then applied graph clustering algorithms to the center population (Methods) to identify recurrently connected modules that may not be easily discerned from peripheral associations. To demonstrate that this graph-theoretic procedure identified biologically relevant modules and to suggest functional roles of these modules, we then examined the patterns of connections between the central and peripheral populations. In the interest of clarity, we present the resultant organization in a hierarchical manner.

At the highest level, the center contains two supermodules: one termed the ‘axial module’, or ‘modA’, which we hypothesize primarily plays a role in movements of the body axis, and a second termed the ‘oculomotor module’, or ‘modO’, which we propose primarily plays a role in eye movements (Fig. 2a,b). RS neurons in the periphery received much stronger connectivity from modA than from modO (Fig. 2a, RS (large), and Extended Data Fig. 1). Of the ten smaller RS neurons contained in the center, all were in modA and none were in modO (Extended Data Fig. 1a, black arrows). By contrast, the 54 peripheral ABD neurons received much stronger connectivity from modO than from modA (Fig. 2a), and all 34 of the Ve² neurons (Fig. 2c) were members of modO. Together, these features suggest that modA is functionally related to the control of body movements, and modO plays a role in eye movements.

**Fig. 2: Modularity and functional specialization of interneurons.**

Although the two populations were quite intermingled, there were notable differences in the spatial distribution of soma positions (Fig. 2c), postsynaptic densities (Fig. 2d) and presynaptic terminals (Fig. 2e). Therefore, we decided to probe to what extent these distributions could be contributing to modularity.

To quantify modularity, we defined an index of wiring specificity as the ratio of the sum of within-module synapse densities to the sum of between-module synapse densities. Here, synapse density is defined as the number of synapses normalized by the product of presynaptic and postsynaptic neuron numbers. To look at the contribution of organization at the level of the neurites, we also defined a ‘potential synapse’ as a presynapse and a postsynapse within some threshold distance of each other (Fig. 2f), similar to the conventional definition as an axo-dendritic apposition within some threshold distance³⁷. Based on actual synapses, the wiring specificity index was roughly 6 for the center neurons (Fig. 2g). The index dropped to less than 3 for potential synapses defined by a distance threshold of 5 μm and close to 2 at a distance threshold of 10 μm (Fig. 2g). Similar trends were observed using a corresponding index for wiring specificity to the periphery (Fig. 2h). Together, these results imply that the division into modA and modO requires synaptic identity to be precisely identified to within a few microns or less⁹.

We also identified, at the population level, statistical differences in other anatomical features of modA and modO neurons and synapses (Extended Data Fig. 2), such as total arbor length and synapse distance from the soma (Extended Data Fig. 2a–c). For completeness, we note that the total number of synapses onto a neuron (regardless of origin) was 310 ± 261 for modA neurons, 174 ± 98 for modO neurons, 1,576 ± 1,555 for peripheral RS neurons and 411 ± 167 for ABD neurons.

The modules in Fig. 2 were obtained with the Louvain algorithm for graph clustering³⁸, which is designed to maximize modularity. Similar modules are obtained when spectral clustering³⁹ is used (Extended Data Fig. 3). As shown next, the binary division into two modules is the first step in a hierarchical procedure that reveals further substructure.

Two major submodules of the oculomotor network

At the next level of organization, the oculomotor center divided into two major submodules, one associated with control of the ipsilateral eye and one associated with control of the contralateral eye (Fig. 3a). The soma, postsynaptic density and presynaptic terminal locations for neurons in each module are shown in Fig. 3b–d. Neurons in the two modules were five times more likely to connect with neurons in the same module than the other module (Fig. 3e). This strong modularity could not be explained solely by coarser-scale proximity of presynaptic terminals (Fig. 3d) and postsynaptic densities (Fig. 3c) as measured by potential synapses, as the wiring specificity ratio dropped to just over 2 when considering potential synapses at a distance of only 2 μm away (Fig. 3e). The relationship of each module to the periphery was assessed by focusing on two neuronal groups within the ABD complex, the ABD motor (ABD_M) neurons that directly drive the lateral rectus muscle of the ipsilateral eye and the ABD internuclear (ABD_I) neurons that indirectly drive the medial rectus muscle of the contralateral eye through a disynaptic pathway; increased activities in both ABD_M and ABD_I neurons drive eye movements toward the side of the brain on which the neurons reside (‘ipsiversive’ movements). Neurons in one submodule of modO preferentially connected to ABD_M neurons, whereas neurons in the other submodule preferentially connected to ABD_I neurons. This is evident from both visual inspection (Fig. 3a) and quantitative analysis (Fig. 3a,f) that revealed a wiring specificity of almost 6 (Fig. 3f). Preferences of individual neurons can be extreme. Many modO cells connect to ABD_M neurons only or ABD_I neurons only (Fig. 3g). We therefore refer to the submodules as motor targeting (modO_M) and internuclear targeting (modO_I) and suggest that they are largely involved in controlling movements of the ipsilateral and contralateral eyes, respectively.

**Fig. 3: Submodules specialized for the two eyes.**

Population-level differences between modO_M and modO_I and their interconnections were also evident in the statistical distributions of synapse size and synapse distance from the soma (Extended Data Fig. 2d,e).

Each oculomotor submodule contains a three-block cycle

The results above suggest the existence of two weakly connected submodules within the oculomotor network. To gain insight into whether these submodules contained substructure, we next examined the eigenvalues of the modO connectivity matrix. The set of eigenvalues of a connectivity matrix, known as the eigenvalue spectrum, contains hallmark signatures of the matrix’s underlying structure. To understand the relationship between the eigenspectrum and network structure, we consider various null models.

When all structure is removed from modO by randomly shuffling connections (Methods), a plot of the real and imaginary parts of the eigenvalues consists of a roughly circular disk centered around 0 plus a single leading eigenvalue along the positive real axis (Fig. 4a, left). Such a pattern of eigenvalues is characteristic of a wide class of matrices with random entries^40,41, with the radius proportional to the standard deviation and the leading eigenvalue proportional to the mean of the matrix entries. If, instead, we keep submodules modO_M and modO_I independent, removing all connections between these submodules, and then shuffle within each submodule, the eigenvalues consist of a roughly circular disk plus two leading eigenvalues (Fig. 4a, middle). Here, the slightly larger leading eigenvalue for modO_I (purple point) relative to modO_M (orange point) indicates that modO_I has slightly stronger overall recurrence. If the connections coupling these submodules are maintained (but separately shuffled), a similar eigenspectrum emerges but with the leading eigenvalues pushed apart from one another along the real axis (Fig. 4a, right). Overall, this null model analysis shows that, if the previously identified submodules modO_M and modO_I contained no further substructure, then we would expect the eigenvalues for the modO connectivity matrix to occupy a roughly circular disk in the complex plane plus have two distinct leading eigenvalues along the positive real axis.

**Fig. 4: ModO contains two distinct cycles with different projections to the ABD.**

The eigenspectrum of the actual modO connectivity matrix is notably richer than expected from the null models (Fig. 4b). Similar to the null model of unstructured submodules (Fig. 4a, right), there is a dense central mass of eigenvalues plus two prominent leading eigenvalues. However, unlike the unstructured submodules, the two leading eigenvalues are each part of a triangle of eigenvalues (Fig. 4b, orange and purple). Such triangles are a hallmark of matrices with a cyclic structure (that is, a set of three blocks connected as 1 → 2 → 3 → 1) and suggest the presence of two weakly coupled cycles of blocks in the modO connectivity matrix. Indeed, our recent analysis of three-cell single cell motifs in cellular connectivity in this system revealed an overrepresentation of three-cell cycles that could support cyclic block structure⁴². We therefore used a stochastic block model (SBM; Methods) to decompose modO into statistically interconnected blocks. We found that the predicted structure of two coupled three-block cycles emerged most clearly when we divided the network into seven blocks, where the seventh block predominantly consisted of Ve² neurons (Fig. 4c).

The first three-block cycle (Fig. 4c, orange: O_M1 → O_M2 → O_M3 → O_M1) connects most strongly to ABD_M neurons, and 83% of its cells come from modO_M. We therefore refer to this cycle as the ‘motor cycle O_M1–O_M3’. The other three-block cycle (Fig. 4c, purple: O_I1 → O_I2 → O_I3 → O_I1) connects predominantly to ABD_I neurons, and 98% of its cells come from modO_I. We therefore refer to this cycle as the ‘internuclear cycle O_I1–O_I3’. Because it predominantly consists of Ve² neurons (32/34), we refer to the seventh block as the vestibular-containing block ‘ ${{\rm{O}}}_{{{\rm{Ve}}}^{2}}$ ’; unlike the three-block cycles, this block projects roughly equally to both the ABD_M and ABD_I populations. This hierarchical division into three-block cycles plus a vestibular-containing block is robust to errors in synapse assignments or definitions for inclusion of neurons in the center population (Extended Data Fig. 4) and leads to a wiring specificity to the periphery of 5.9, similar to that for the division into only two submodules of Fig. 3f.

To better isolate the correspondence between the eigenvalue triangles and the cycles seen in the seven-block organization, we decoupled the motor and internuclear three-block cycles by removing all connections to and from the internuclear block. This preserved the orange and purple triangles of eigenvalues, but, as was seen for the coupled versus uncoupled modO_M and modO_I (Fig. 4a, middle), the biggest effect of removing the coupling is to bring the two leading eigenvalues closer together (Fig. 4d). Additionally shuffling the connections within the internuclear cycle preserves the orange triangle in the eigenspectrum while removing the purple triangle (Fig. 4d), enabling us to identify the purple triangle with the internuclear cycle and the orange triangle with the motor cycle (Fig. 4e). A graphical summary of modO recurrent connectivity and connectivity to the periphery is provided in Fig. 4f.

The coupling between cycles, although relatively weak, may have important functional roles. The effect of this coupling can be seen in the eigenvectors corresponding to the two leading eigenvalues (Fig. 4g, top). These eigenvectors represent the patterns across the neurons of the network that have the first and second strongest recurrent feedback connectivity. When the cycles are decoupled, one modO eigenvector has nonzero entries only for neurons in the O_I1–O_I3 cycle (Fig. 4g, bottom left), whereas the other contains nonzero entries only for neurons in the O_M1–O_M3 cycle and ${{\rm{O}}}_{{{\rm{Ve}}}^{2}}$ block (Fig. 4g, bottom right). Including the ABD neurons in the eigenvector analysis (which does not change the eigenvalues of Fig. 4b,d) further shows that the leading decoupled eigenvectors associate predominantly with ABD_M and ABD_I neurons, respectively (Fig. 4g, bottom). By contrast, when the full coupling is included, the leading eigenvectors approximately combine the eigenvectors found in the decoupled case. The leading eigenvector approximates a scaled sum of the two decoupled eigenvectors, and the second leading eigenvector approximates a scaled difference of the decoupled eigenvectors. Biologically, this suggests that activation of the pattern of modO_M neurons corresponding to the first eigenvector may generate ‘conjugate eye movements’ (eyes moving together in the same direction, which requires coactivation of the ABD_M and ABD_I populations), whereas activation of the second eigenvector may be critical for maintaining either disconjugate eye movements (eyes moving in opposite directions) or deviations from conjugate eye movements.

Two of these graph-theoretically derived blocks would not be easily identifiable as part of modO by naively tracing connections back from the periphery. Block O_M1 of the motor submodule projects extremely weakly to the motor neurons but has the strongest projections to and receipt of projections from the internuclear submodule. Thus, it appears to be involved nearly solely in internal processing. Similarly, block O_I3 of the internuclear cycle has relatively little connection to the ABD. Other blocks, such as block O_M3, appear to divide into two halves with different functions for each half; for block O_M3, only half projects to the ABD, whereas the other half provides the strongest projection to the ${{\rm{O}}}_{{{\rm{Ve}}}^{2}}$ block.

We separately applied an SBM to modA and examined how axial submodules aligned with those from modO. We found three axial submodules, with strong connections within and relatively weak connections between the submodules. Coupling between the oculomotor and axial submodules was weak overall, with the strongest coupling being mediated by connections to and from axial block A1 (Extended Data Fig. 5).

Predicting neural coding and dynamics from the wiring diagram

We hypothesized that modO contains the majority of the ‘neural integrator’ for horizontal eye movements that transforms angular velocity inputs into angular position outputs. This transformation has been suggested to depend on a relatively high degree of recurrent interactions within the velocity-to-position neural integrator (VPNI). However, it is not known how these interactions shape the single-neuron coding properties of neurons within modO and its peripheral outputs in the ABD.

We developed a network model to assess whether functional properties of modO and ABD neurons could be predicted from connectomic information. Within modO, we distinguished between generic modO neurons and the exceptions (Ve²) that received inputs from primary vestibular afferents (Fig. 5a). We first focused on predicting the relative strength of the eye position signal in these different populations during spontaneous fixation behavior and then examined what aspects of the firing rate dynamics could be predicted.

**Fig. 5: Connectome-based model captures functional characteristics of the oculomotor system.**

A minimal set of physiological constraints were imposed on the model. Based on alignment with previous studies (see Methods and Identification of excitatory versus inhibitory neurons), we assumed that all Ve² cells were inhibitory (Supplementary Fig. 3c and Methods)³² and the remaining cells were excitatory²¹ (Methods). We omitted interactions with contralateral neurons based on physiological studies indicating that persistent activity in the oculomotor integrator is generated independently by each half of the brain^20,43,44. Likewise, we omitted interactions between modO and modA based on their separate roles in controlling eye movement versus body movement behaviors^34,35,45 and the body-fixed nature of the oculomotor behavior being investigated.

We determined the model’s weight matrix from our EM reconstruction as follows. Each element W_ij of the weight matrix is the physiological strength of the connection received by neuron i from neuron j. First, we approximated W_ij as proportional to N_ij, the number of synapses received by i from j. This approximation effectively assumes that all individual synapses involved in a connection have the same strength. For each neuron i, we then divided by a normalization factor Σ_kN_ik equal to the total number of synapses received by neuron i, including those from neurons outside modO. This normalizing factor has two motivations. First, normalizing by the total number of incoming synapses provides a form of homeostatic synaptic scaling⁴⁶. Second, the normalization provides a simple way to compensate for truncation of dendritic arbors by the borders of the EM volume. With this normalization, each W_ij is proportional to the fraction of the total inputs that neuron j provides to neuron i. Finally, to set the proportionality constant, we did an overall scaling of the weight matrix by a factor β to make the network response exhibit persistent activity (Methods). Therefore, the final synaptic weight matrix took the form ${W}_{ij}=\pm \beta \frac{{N}_{ij}}{{\sum }_{k}{N}_{ik}}$ .

We then constructed a linear network model with weight matrix W_ij (Methods). Linear model neurons⁴⁷ were used for simplicity because most VPNI cells have firing rates that vary linearly with position for ipsilateral eye positions^48,49.

We first used this model to predict the position sensitivities (that is, changes in firing rates for a given change in eye position) of all modO and ABD neurons. To form these predictions, we simulated the response to saccadic burst input and recorded the firing rates of model neurons during the subsequent fixation (Extended Data Fig. 6). We then tested these model predictions against prior calcium imaging experiments⁵⁰ (Methods). We assessed experimental eye position sensitivity during spontaneous fixations in the dark (Fig. 5a–c). Because the imaged neurons come from different fish than the reconstructed one, they cannot be placed in one-to-one correspondence with neurons in our EM volume. Therefore, we resorted to a population-level comparison.

Each of the four simulated populations (Ve², generic modO, ABD_M and ABD_I) displayed a characteristic distribution of eye position sensitivities (Fig. 5d, gray). The Ve² population exhibited little eye position sensitivity. The generic modO population that contains the VPNI neurons exhibited a long-tailed distribution of eye position sensitivities with much higher mean sensitivity. The ABD population generated the largest eye position sensitivities, and these were bimodally distributed, with higher sensitivities for ABD_I than for ABD_M neurons. Qualitatively, the distribution of experimental eye position sensitivities for each population matched the model predictions quite well (Fig. 5d, green), although the bimodality seen in the model ABD populations was not clearly evident. We compared these simulations to three classes of neurons (VPNI, Ve² and ABD) in the experiments, which could not easily distinguish between ABD_M and ABD_I neurons.

To test whether potential connectivity would have been sufficient for our network modeling, we estimated N_ij by the number of potential synapses onto neuron i from neuron j and then computed the weight matrix W_ij by the normalization described above (Methods). When our weight matrix based on potential synapses was substituted into our network model, we obtained population distributions for eye position sensitivities that were considerably different from experimental measurements (Fig. 5e). This suggests that the functional properties of the circuit depend on the spatially precise connections revealed by the EM connectome.

We next tested whether the connectome could predict the low dimensionality of the system dynamics and any features of the distribution of single-neuron relaxation dynamics across the population. Regarding the low dimensionality, because the (real parts of the) eigenvalues of the weight matrix set the relative timescales of decay in the system and the matrix W_ij has two leading eigenvalues, this immediately predicts that the postsaccadic activity of the network will exhibit low-dimensional dynamics at long times (Fig. 5f; percent variance explained by the leading principal components of recorded and simulated network activity). To assess the dynamics of exponential decay of the network activity, we performed a double exponential fit to each neuron in the model and experimental recordings (Fig. 5c, bottom). We then calculated the amplitude of each neuron’s participation in the longest and second-longest exponential decays. Consistent with the experiments and previous findings¹⁶, the model predicted a negative correlation between the decay amplitudes (Fig. 5g).

To test the robustness of these results, we changed the model by adjusting the threshold for dividing neurons into center and periphery (Methods and Extended Data Fig. 7). Even if the number of neurons in the center was reduced by 25%, the model distributions of eye position sensitivities and exponential decay dynamics as well as the dimensionality of the network remained similar (Extended Data Fig. 7a,c). We also corrupted our reconstructed wiring diagram by simulating errors in the automated synapse detection and likewise found that the population distributions of eye position sensitivities and exponential decay dynamics and the dimensionality of the network remained very similar (Extended Data Fig. 7b,d).

Remarkably, the match between model predictions and experiments in Fig. 5d,f,g involved fitting only a minimal set of parameters. For the eye position sensitivity distributions of Fig. 5d, we set an overall scale factor such that the mean eye position sensitivity for the generic modO population in the model matched that of the experimental VPNI population; beyond this, as long as the intrinsic time constant was set short enough that the eye position sensitivity was dominated by the leading eigenvector of the weight matrix, the shapes of the distributions matched. Quantitatively matching the approximately two-dimensional persistent activity of the experimentally recorded population (Fig. 5f,g) required additionally setting the intrinsic time constant τ and the value of the overall weight scale parameter β in the model. Interestingly, we found that the best fits were obtained with an intrinsic time constant of ~1 s, a timescale estimated previously from inactivation experiments⁴³.

Finally, we examined which features of the wiring diagram are important for generating the above results. We did this by randomly shuffling the connections within a one-, two- and seven-block model of modO. The one-block model failed to capture the eye position sensitivity distributions, the network dynamics dimensionality and the exponential decay dynamics distributions, whereas both the two- and seven-block models reasonably captured the experimental findings (Extended Data Fig. 8a–c, top). The eye position sensitivities in the seven-block model were slightly less sensitive to shuffling than those in the two-block model (Extended Data Fig. 8a–c, middle). The observation that the two- and seven-block models behaved very similarly likely reflects that our experiments probed only the longest timescales of behavior (the two leading eigenvalues), and the two- and seven-block models have highly similar leading eigenvalues (Extended Data Fig. 8a–c, bottom, and Extended Data Fig. 9). Thus, although the seven-block model is required to recover the cyclic connectivity of the network and the separate vestibular-dominated block ${{\rm{O}}}_{{{\rm{Ve}}}^{2}}$ , it is not required to reveal the slowest timescale decay dynamics.

Discussion

This work shows how the neural coding properties throughout a hindbrain network relate to its synapse-resolution structure. Through connectomic reconstruction and analysis, we reveal that hindbrain regions previously thought to be largely unstructured instead contain a hierarchical modular organization that closely aligns with known motor functions. Directly inserting the connectome into a network model statistically predicts the single-neuron coding properties assessed with two-photon calcium imaging throughout the reconstructed volume, representing a breakthrough in deriving neuronal properties within a topographically heterogeneous network from the connectome. EM-resolution analysis was critical to revealing the above network structure and predictions, as there were large discrepancies when lower-resolution analysis based on potential synapses was used to assess connectivity. The resulting connectome and modular analysis tools are being provided as an open resource to the community.

Low-dimensional brain dynamics have been found to underlie a wide array of motor, navigational and cognitive functions. Twenty-five years ago, ‘line attractor’ and ‘ring attractor’ network models were proposed for low-dimensional neural dynamics in the oculomotor¹⁵ and head direction systems^51,52, respectively. Connectomic information from the Drosophila head direction system¹⁰ is currently being used to inform ring attractor network models. Our work similarly informs line attractor network models, for which the VPNI has served as a model system^53,54. The classic view of this system is as a line attractor with unstructured recurrent feedback leading to persistent activity¹⁵. This picture corresponds well to the one-block shuffled connectivity of our network analysis (Fig. 4a, left) but does not fit well with the true connectome. Rather, we found the connectome for the VPNI to contain two weakly coupled three-block cycles that lead to a corresponding pair of large real eigenvalues, providing an anatomical basis for the low-dimensional neural activity observed in this network^16,55. Future connectomes incorporating the different classes of anatomical inputs targeting the VPNI may additionally explain the previously observed context dependence of these low-dimensional dynamics¹⁶.

The reconstruction and analysis of modO provides further insights into the organization and function of the VPNI and its relation to the larger oculomotor system. The presence of two weakly coupled VPNI subnetworks, preferentially projecting either to ABD_M neurons (modO_M) or ABD_I neurons (modO_I), likely enables separate control of the two muscles (lateral rectus and medial rectus, respectively) controlling horizontal eye movements. Together, these subnetworks could provide for maintenance of conjugate and disconjugate eye movements^56,57,58. The identified cyclic structure implies a pair of complex conjugate eigenvalues in the eigenspectrum (Fig. 4a) that, dynamically, corresponds to overdamped oscillations (Extended Data Fig. 10). This may be required for control of the eye plant, which can be approximated as a highly damped mass-spring system whose impulse-response function exhibits overdamped oscillations⁵⁹ and may suggest the importance of a weak, but non-negligible, mass term in the plant that is commonly ignored in simplified descriptions^55,60. Finally, at a gross anatomical level, previous physiological studies of VPNI cells mainly focused on rhombomere 7/8 (R7/R8; refs. ^16,21,22,55). Our map of modO (Fig. 2c) suggests that, at least at this stage of development, the VPNI should also include R4–R6. This extension is consistent with a previous observation that VPNI function was only partially abolished by sizable inactivation of R7/R8 (ref. ⁵⁵) and with previous reports of eye position signals in R4–R6 neurons that are not from ABD neurons^50,57.

In most neural network models of brain function, the synaptic weight matrix has been regarded as the solution to an ‘inverse problem’. Given the observed effects (neural activity and behavior), the modeler attempts to identify the unobserved cause (weight matrix). Connectomics offers the possibility of treating network modeling as more of a ‘forward problem’. The forward approach has been feasible for small nervous systems in which the weight matrix can be directly observed and completely mapped by synaptic physiology⁶¹. In sensory systems, a forward approach has been used to demonstrate how connectomic information could be used to model whitening of odor representations in a vertebrate olfactory bulb⁵ and to constrain a model of orientation and direction selectivity in the Drosophila visual motion detection circuit⁶², although fine-tuning by backpropagation learning was necessary in the latter case. At a lower, ‘mesoscopic’ level of resolution, interarea projection maps in primates have been used to explain the temporal dynamics of cortical responses⁶³. Here, we have shown how, with a minimal number of assumptions, the forward modeling approach could predict the statistical distribution of single-neuron eye position sensitivities and exponential decay dynamics for several neural populations measured by calcium imaging of animals during ocular fixations (Fig. 5c). Some ‘inverse’ aspects to our modeling remained because we constrained ourselves to modeling eye movements in the absence of body movements. Our forward approach also incorporated a minimal set of key anatomical and physiological constraints such as known signs of connections^21,32, physiological observations about the approximate linearity of oculomotor responses^18,49 and functional independence of the bilateral halves of the circuit^20,56. Nevertheless, given the minimal nature of these physiological constraints, our success in matching physiological observations is perhaps surprising. This is especially so given our naive estimates of synaptic weights from the simple measure of number of synapses, simplified linear rate model neuron treatment of cell morphology, synaptic and intrinsic cellular biophysics and neglect of possible neuromodulatory effects⁶⁴. Comparisons between model predictions and physiological data at the level of single cells, rather than populations, might require more sophisticated modeling of cell- and synapse-specific biophysics.

As more connectomes become available in other settings, it will be important to consider which physiological constraints need to be incorporated to make appropriate use of these powerful datasets. Answers to this question depend on many factors, including the breadth of behaviors to be produced in a single model, the range of dynamics of the constituent components and the degree to which the model is to produce quantitative versus qualitative matches to data. Nevertheless, we hope that our work suggests how, when guided by knowledge of the various behaviors a circuit participates in and appropriate physiological constraints gleaned from recordings and perturbations of activity, it may be possible in even more complex circuits to identify physiological modules whose function can be well understood using the connectome-based analysis and modeling approach taken here.

Methods

Image acquisition and alignment

We acquired a dataset of the larval zebrafish hindbrain (mitfa^b692/+ (ABxWIK), 6 days after fertilization, sex unspecified, Zebrafish International Resource Center) that extends 250 μm rostrocaudally and includes R4 to R7/R8. The volume extends 120 μm laterally from the midline and 80 μm ventrally from the plane of the Mauthner cell axon. The serial-section Electron Microscopy dataset was an extension of the original dataset in Vishwanathan et al.²² and was extended by additional imaging of the same serial sections. The images were acquired using a Zeiss Field-Emission Scanning Electron Microscope using WaferMapper software⁶⁵. Only a few tens of neurons had been manually reconstructed in our original publication on the ssEM dataset in Vishwanathan et al.²². The dataset was stitched and aligned using the custom package Alembic (see Code availability). The tiles from each section were first montaged in two dimensions (two steps) and then registered and aligned in 3D as whole sections (four steps). Point correspondences for alignment were generated by block matching via normalized cross-correlations both between tiles and across sections. Errors in each step were found by a combination of programmed flags (such as lower than expected correspondences, small search radius, large distribution of norms or high residuals after mesh relaxation) and visual inspection. They were corrected by either changing the parameters or by manual filtering of points. In most cases, the template and the source were both passed through a bandpass filter. The final parameters that were used are listed in Table 1.

Table 1 Parameters used for image alignment

Full size table

Stitching of tiles (montaging) within a single section was split into a linear translation step (premontage) and a nonlinear elastic step (elastic montage). In the premontage step, individual tiles were assembled to have 10% overlap between neighboring tiles, as specified during imaging and by fixing a single tile (anchoring) in place. They were then translated row by row and column by column according to the signal correspondence found between the overlaps. In the elastic montage step, the locations of the tiles were initialized from the translated locations found previously, and blockmatches were computed every 100 pixels on a regular triangular mesh (see Table 1 for the parameters used). Once the correspondences were found, outliers were filtered by checking the spatial distribution of the cross-correlogram (σ filter), height of the peak of the correlogram (Pearson r value), dynamic range of the source patch contrast, kurtosis of the source patch, local consensus (average of immediate neighbors) and global consensus (inside the section). After the errors had been corrected by filtering bad matches, the tiles were modeled as individual spring meshes attached via the matches⁶⁶, and the linear system was solved using conjugate gradient descent. The mean residual errors were in the range of 0.5–1.0 pixels after relaxation.

The intersection alignment was split into a translation step (pre-prealignment), a regularized affine step (prealignment), a fast coarse elastic step (rough alignment) and a slow fine elastic step (fine alignment). In the pre-prealignment step, a central patch of the given montaged section was matched to the previous montaged section to obtain the rough translation between two montaged sections. In the prealignment step, the montaged images were offset by that translation, and then a small number of correspondences were found between the two montaged sections. The sections were then aligned using an affine transformation fit using least squares, regularized with a 10% (empirically derived) weighting of the identity transformation to reduce shear accumulating and propagating across multiple sections. Proceeding sequentially allowed the entire stack to get roughly in place. The mean residual errors after prealignment were in the range of 3.5 pixels after relaxation. In the subsequent rough and fine alignment steps, the blockmatches were computed and filtered between bandpassed neighboring sections in a regular triangular mesh, similar to the elastic montaging step, followed by a conjugate gradient descent on the linear system. Splitting the elastic step into rough and fine alignments allowed the search radius to be reduced during the fine alignment relative to what would have been necessary with a direct attempt, thereby lowering the likelihood of spurious matches as well as the computing time.

Convolutional net training

Dataset

Four expert brain image analysts (D. Visser, Princeton Neuroscience Institute; K.W.; M. Moore, Princeton Neuroscience Institute; and S.K.) manually segmented neuronal cell boundaries from six subvolumes of EM images with the manual annotation tool VAST⁶⁷, labeling 194.4 million voxels in total. These labeled subvolumes were used as the ground truth for training convolutional networks to detect neuronal boundaries. We used 187.7 million voxels for training and reserved 6.7 million voxels for validation.

Network architecture

To detect neuronal boundaries, we used a multiscale 3D convolutional network architecture similar to the boundary detector in Zung et al.⁶⁸. This architecture was similar to U-Net⁶⁹ but with more pathways between scales. We augmented the original architecture of Zung et al.⁶⁸ with two modifications. First, we added a ‘lateral’ convolution between every pair of horizontally adjacent layers (that is, between feature maps at the same scale). Second, we used batch normalization⁷⁰ at every layer (except for the output layer). These two architectural modifications were found to improve boundary detection accuracy and stabilize/speed-up training, respectively. For more details, we refer the reader to Zung et al.⁶⁸.

Training procedures

We implemented the training and inference of our boundary detectors with the Caffe deep learning framework⁷¹. We trained the networks on a single Titan X Pascal GPU. We optimized the binary cross-entropy loss with the Adam optimizer⁷² initialized with α = 0.001, β₁ = 0.9, β₂ = 0.999 and ε = 0.01. The step size α was halved when the validation loss plateaued after 135,000, 145,000 and 175,000 iterations. We used a single training example (minibatch of size 1) to compute gradients for each training iteration. The gradient for target affinities (the degree to which image pixels are grouped together) in each training example was reweighted dynamically to compensate for the high imbalance between target classes (that is, low and high affinities). Specifically, we weighted each affinity inversely proportional to the class frequency, which was computed independently within each of the three affinity maps (x, y and z) and dynamically in each training example. We augmented training data using (1) random flips and rotations by 90°, (2) brightness and contrast augmentation, (3) random warping by combining five types of linear transformation (continuous rotation, shear, twist, scale and perspective stretch) and (4) misalignment augmentation⁷³ with a maximum displacement of 20 pixels in the x and y dimensions. The training was terminated after 1 million iterations, which took about 2 weeks. We chose the model with the lowest validation loss at 550,000 iterations.

Convolutional net inference

The above trained network was used to produce an affinity map of the whole dataset using the ChunkFlow.jl package⁷⁴. Briefly, the computational tasks were defined in a JSON formatted string and submitted to a queue in Amazon Web Services Simple Queue Service (AWS SQS). We launched 13 computational workers locally with NVIDIA Titan X GPU. The workers fetched tasks from the AWS SQS queue and performed the computation. The workers first cut out a chunk of the image volume using BigArrays.jl and decomposed it into overlapping patches. The patches were fed into the convolutional network model to perform inference in PyTorch⁷⁵. The output affinity map patches were blended in a buffer chunk. The output chunk was cropped around the margin to reduce boundary effects. The final affinity map chunk, which is aligned with block size in cloud storage, was uploaded using BigArrays.jl. Both image and affinity map volumes were stored in Neuroglancer precomputed format (https://neurodata.io/help/precomputed/). The inference took about 17 days in total and produced about 26 TB of affinity map.

Chunk-wise segmentation

To perform segmentation of the entire volume, we divided the volume into ‘chunks’. Overlapping affinity map chunks were cut out using BigArrays.jl, and a size-dependent watershed algorithm⁷⁶ was applied to agglomerate neighboring voxels to make supervoxels. The agglomerated supervoxels are represented as a graph where the supervoxels are nodes and the mean affinity values between contacting supervoxels are the edge weights. A minimum spanning tree was constructed from the graph by recursively merging the highest weight edges. This oversegmented volume containing all supervoxels and the minimum spanning tree was ingested into Eyewire (https://eyewire.org) for crowdsourced proofreading.

Semiautomated reconstructions on Eyewire

Neurons were chosen for proofreading in Eyewire based on an initial set of ‘seed’ neurons that were identified as carrying eye position signals by co-registering the EM volume to calcium imaging performed on the same animal²². All pre- and postsynaptic partners of the initial seed of 22 neurons were reconstructed. Following this, we reconstructed partners of the neurons that were reconstructed in the initial round in a random manner. Eyewirers were provided the option of agglomerating (merging) supervoxels using a slider to change the threshold of agglomeration. To ensure accurate reconstructions, we did two things: (1) only players who met a certain threshold (see below), determined by their accuracy on a previously published retinal dataset^14,77, were allowed to reconstruct zebrafish neurons, and (2) the reconstructions were performed by two players in two rounds, in which the second player could modify the first player’s reconstruction (Supplementary Information). Finally, after two rounds of reconstruction, neuronal reconstructions were checked and if necessary corrected by expert in-house image analysts who each had more than 5,000 h of experience. The accuracy of the players in the crowd compared to experts (assuming experts are 100%) was >80% in the first round and ~95% after the second round of tracing. The validated reconstructions were subsequently skeletonized for analysis purposes. Player accuracy was calculated as an F1 score, where F1 = 2TP/(2TP + FP + FN), where TP represents true positives, FP represents false positives, and FN represents false negatives. All scores were calculated as a sum over voxels. TP was assigned when both the player and the expert agreed that the segmentation was correct, FN was assigned when the player missed segments that were added in by the expert, and FP was assigned when the player erroneously added segments that did not belong. Two F1 scores were calculated for each player, once for round 1 and once for round 2. No player played the same neuron in both rounds. Typically, at an agglomeration threshold of 0.3, the segmentation had an F1 score of 62%.

Only the most experienced players on Eyewire were allowed to participate. A small group of four highly experienced players received an invitation to test this new dataset. The players were given the title of ‘Mystic’, which was a new status created to enable gameplay and became the highest achievable rank in the game. Subsequent Mystic players had to apply for the status to unlock access to the zebrafish dataset within Eyewire. There was a high threshold of achievement required for a player to gain Mystic status. Each player was required to reach the previous highest status within the game as well as complete 200+ cubes a month and maintain 95% accuracy when resolving cell errors. Once a player was approved by the lab, they were granted access to the zebrafish dataset and given the option to have a live tutorial session with an Eyewire Admin. There was also a prerecorded tutorial video and written tutorial materials for players who could not attend the live session or who desired a review of the materials. Newly promoted players were also given a special badge, a new chat color and access to a new chat room exclusive to Mystics. These rewards helped motivate players by showing their elevated status within the game as well as giving them a space to discuss issues specific to the zebrafish dataset.

Cells were parceled out in batches to players. When a cell went live on Eyewire, it could be claimed by any Mystic. Each cell was reconstructed by only one player at a time. Once this player had finished their reconstruction, a second player would check over the cell for errors. After the second check, a final check was done by an Eyewire Admin. To mitigate confusion when claiming cells, a special graphical user interface was built into the game that allowed players to see the status of each cell currently online. A cell could be at one of the following five statuses: ‘Need Player A’, ‘Player A’, ‘Need Player B’, ‘Player B’ and ‘Need Admin’. These statuses indicated whether a cell needed to be checked or was in the process of being checked and whether it needed a first-, second- or Admin-level check. At each stage, the username of the player or Admin who had done a first, second or final check was also visible. It was made mandatory that the first and second checks were performed by two separate players.

Collaboration and feedback were important parts of the checking process. If a player was unsure about an area of the cell they were working on, they could leave a note with a screenshot and detail of the issue or create an alert that would notify an Admin. If a ‘Player B’ or an Admin noticed a mistake made earlier in the pipeline, they could inform the player of the issue via a ‘Review’ document or through an in-game notification (player tagging).

To differentiate the zebrafish cells from the regular dataset, each cell was labeled with a Mystic tag. This tag helped identify the cells as separate from the e2198 retinal Eyewire dataset and also populated them to a menu of active zebrafish cells within Eyewire.

Players were rewarded for their work in the zebrafish dataset with points. For every edit they made to a cell, they received 300 points. Points earned while playing zebrafish were added to a player’s overall points score for all gameplay done on Eyewire and appeared in the Eyewire leaderboard.

Skeletonization

The neuron segmentation IDs were ingested to an AWS SQS queue, and multiple distributed workers were launched in Google Cloud using kubernetes. Each worker fetched the segmentation chunks associated with a neuron ID. The segmentation voxels were extracted as a point cloud, and the distance to boundary field (DBF) was computed inside each chunk. Finally, a modified version of the skeletonization algorithm TEASAR was applied⁷⁸. Briefly, we constructed a weighted undirected graph from the point cloud, where the neighboring points are connected with an edge and the edge weight is computed from the DBF. We then took the point with the largest DBF as source and found the furthest point as the target. The shortest path from source to target in the graph was computed as the skeleton nodes. The surrounding points were labeled as visited, and the closest remaining unvisited point was taken as the new source. We repeated this process until all the points were visited. The skeleton node diameter was set as its DBF. The skeleton nodes were postprocessed by removing redundant nodes, removing ‘hairs’ based on diameter, removing branches inside somata, downsampling the nodes, merging single-child segments and smoothing the skeleton path. All skeletonization was performed at multum in parvo level 4.

Synapse detection

Synapses were automatically segmented in this dataset using neural networks to detect clefts and assign the correct partner as previously described³⁰. Briefly, a subset of the imaged data (219 μm³) was selected for annotation. The annotations were performed using the manual annotation tool VAST⁶⁷. Trained human annotators labeled the voxels that were part of the postsynaptic density and presynaptic docked vesicle pools. A convolutional neural network³⁰ was trained to match the postsynaptic density, using 107 μm³ as a training set and 36 μm³ as a validation set and leaving the remaining 76 μm³ as an initial test set. All of these sets were compared to the predictions of the model tuned to an F-score of 1.5 on the validation set to bias toward recall, where recall = TP/(TP + FN). Biasing the predictor toward recall reduces false negatives at the cost of more false positives, which are easier to correct. Apparent human errors were corrected, and training was restarted with a new model. We also later expanded the test set by proofreading similar automated results applied to new sections of the datasets (to increase representation of rare structures in the full image volume). The final model used an RS-U-Net architecture⁷³ implemented using PyTorch⁷⁵ and was trained using a manual learning rate schedule, decreasing the rate by a factor of 10 when the smoothed validation error converged. The final network reached 86% precision and 83% recall on the test set after 230,000 training iterations.

A convolutional network was also trained to assign synaptic partners to each predicted cleft as previously described³⁰. All 361 synapses in the ground truth were labeled with their synaptic partners, and the partner network used 204 synapses as a training set, 73 as a validation set and the remaining 84 as a test set. The final network was 95% accurate in assigning the correct partners of the test set after 380,000 training iterations.

The final cleft network was applied across the entire image volume and formed discrete predictions of synaptic clefts by running a distributed version of connected components. Each cleft was assigned synaptic partners by applying the partner network to each predicted cleft within nonoverlapping regions of the dataset (1,024 × 1,024 × 1,792 voxels each). In the case where a cleft spanned multiple regions, the assignment within the region that contained most of that cleft was accepted, and the others were discarded. Cleft regions whose centroid coordinates were within 1 μm and were assigned the same synaptic partners were merged together to merge artificially split components.

Finally, spurious synapse assignments (that is, postsynapses on axons and presynapses on dendrites) were cleaned by querying the identity of the ten nearest synapses to every synapse, where each synapse was associated with its closest skeleton node on both the pre- and postsynaptic sides. If the majority of the ten nearest neighbors were of the same identity (presynaptic or postsynaptic), then the synapse was assigned correctly. If the majority were of an opposing identity, these synapses were assigned incorrectly and were deleted. This process eliminated 1,975 falsely assigned synapses (~2% of the total).

Registration to reference atlas

Registration of the EM dataset to the Z-Brain reference atlas³¹ was performed in two stages. We created an intermediate EM stack from the low-resolution (270 nm per pixel) EM images of the entire larval brain tissue. This intermediate stack had the advantage of a similar field of view as the light microscopy reference volume while also being of the same imaging modality as the high-resolution EM stack. The low-resolution EM stack was registered to the reference brain (elavl3:H2b-RFP) by fitting an affine transform that maps the entire EM volume onto the light microscopy volume. To do this, we selected corresponding points such as neuronal clusters and fiber tracts using the tool BigWarp⁷⁹. These corresponding points were used to determine an affine transform using the MATLAB least squares solver (mldivide). Subsequently, the intermediate EM stack in the same reference frame as the Z-Brain atlas was used as the template to register the high-resolution EM stack. This was performed in a similar manner by selecting corresponding points and fitting an affine transform. This transform was used to map the reconstructed skeletons from high-resolution EM space to the reference atlas space.

Identification of ABD neurons

We identified ABD_M neurons (Fig. 1), which drive lateral movements of the ipsilateral eye, by their overlap with the Mnx transgenic line (Supplementary Fig. 3a, left). ABD_M axons exited R5 and R6 through the ABD (VIth) nerve (Fig. 1f, black box) as reported previously²².

Medial movements of the contralateral eye are driven by the medial rectus muscle, which is innervated by motor neurons in the contralateral oculomotor nucleus, which in turn are driven by internuclear neurons (ABD_I) in the ipsilateral ABD complex (Fig. 1f). ABD_I neurons were identified (Supplementary Fig. 3a, right) by their overlap with two nuclei in the Evx2 transgenic line that labels glutamatergic interneurons. The ABD_I neurons were just dorsal and caudal to the ABD_M neurons, and their axons crossed the midline⁸⁰.

Identification of Ve² neurons

We identified a class of Ve² neurons (Fig. 1h, dark brown) that received synapses from primary vestibular afferents. The latter were orphan axons in R4 identified as the vestibular branch of the vestibulocochlear nerve (VIIIth nerve) by comparison with the Isl2 line, which labels the major cranial nerves (Supplementary Fig. 3b, blue axons). We speculate that Ve² neurons are descending octavolateral neurons but withhold final annotation until functional confirmation is available³².

Identification of RS neurons

RS neurons were divided into large and small groups (Fig. 1h). Large RS neurons were the M, Mid2, MiM1, Mid3i and CaD neurons. Small RS neurons were RoV3, MiV1 and MiV2 cells. These were identified by their stereotypic locations⁸¹ and by comparison with the Z-Brain atlas (Supplementary Fig. 4).

Centrality-based division into the center and periphery

The division of the reconstructed wiring diagram into center and periphery is based on standard measures of ‘centrality’, which have been devised in network science⁸². We define the simplest measure, known as degree centrality, as the geometric mean of the in-degree and out-degree of a neuron (it is more common to choose one or the other). Another popular measure, known as eigenvector centrality or eigencentrality, is a neuron’s element in the eigenvector of the connection matrix N_ij (the number of synapses onto neuron i from neuron j) corresponding to the eigenvalue with maximal real part; this measure extends the simpler concept of degree centrality by weighing a neuron’s connections by its centrality (a neuron is more central to the network if it receives inputs from other high-centrality neurons). Mathematically, this defines the eigenvector problem ν_i = ΣN_ijν_j, where ν_i is the (input) centrality of neuron i. An analogous formula, but instead replacing N_ij by its transpose (thus, defining left rather than right eigenvectors), can be used to weight output connections by the (output) centrality of the neuron to which output is projected. The eigenvector elements can be chosen non-negative by the Perron–Frobenius theorem. It is standard to use either the left or right eigenvector, but we use both for our definition by computing the geometric mean of the left and right eigenvector elements. Degree centrality and eigencentrality are correlated but not perfectly (Supplementary Fig. 1c). For a visualization of the network based on eigencentrality, see Supplementary Fig. 1e.

We define the periphery as those cells with vanishing (<10⁻⁸) eigencentrality but nonzero degree centrality. All but 62 of these 2,344 neurons had strictly 0 eigencentrality. The remaining 540 recurrently connected neurons are defined as the ‘center’ of the graph. See below for the effects of varying the eigencentrality threshold for center–periphery division.

Graph clustering

We applied two graph clustering algorithms to divide the center into two supermodules and obtained similar results from each algorithm. The clustering from the Louvain algorithm is presented in the main text and that of the spectral algorithm in Extended Data Fig. 3. We further subdivided each supermodule into subclusters using an SBM, as described below.

Louvain clustering

Graph clustering was performed using the Louvain clustering algorithm⁸³ for identifying different ‘communities’ or ‘modules’ in an interconnected network by optimizing the ‘modularity’ of the network, where modularity measures the (weighted) density of connections within a module compared to between modules. Formally, the modularity measure maximized is Q_gen = ∑B_ijδ(c_i, c_j), where δ(c_i,c_j) equals 1 if neurons a and b are in the same module and 0 otherwise, where ${B}_{ij}=\frac{1}{\omega }({W}_{ij}-\gamma \frac{{s}_{i}{s}_{j}}{\omega })+{\mathrm{transpose}}[\frac{1}{\omega }({W}_{ij}-\gamma \frac{{s}_{i}{s}_{j}}{\omega })]$ . Here, s_i = ∑_cW_ic is the sum of weights into neuron i, s_j = ∑_bW_jb is the sum of weights out of neuron j, ω = ∑_cdW_cd is the total sum of weights in the network, and the resolution parameter γ determines how much the naively expected weight of connections $\gamma \frac{{s}_{i}{s}_{j}}{\omega }$ is subtracted from the connectivity matrix. Potential degeneracy in graph clustering was addressed by computing the consensus of the clustering similar to Sporns and Betzel⁸⁴. Briefly, an association matrix, counting the number of times a neuron is assigned to a given module, was constructed by running the Louvain algorithm 200 times. Next, a randomized association matrix was constructed by permuting the module assignment for each neuron. Reclustering the thresholded association matrix, where threshold was the maximum element of the randomized association matrix, provided consensus modules. We used the community_louvain.m function from the Brain Connectivity Toolbox package (https://sites.google.com/site/bctnet/home). In addition to the Louvain graph clustering algorithm, we also clustered the ‘center’ with an alternative graph clustering algorithm, spectral clustering, described below.

Spectral clustering

We used a generalized spectral clustering algorithm for weighted directed graphs to bisect the zebrafish ‘center’ subgraph, as proposed by Chung³⁹. Given a graph G(V, E) and its weighted adjacency matrix $A\in {{\mathbb{R}}}_{\ge 0}^{n\times n}$ , where A_ij indicates the number of synapses from neuron i to neuron j, one can construct a Markov chain on the graph with a transition matrix P_α, such that ${[{P}_{\alpha }]}_{ij}:= (1-\alpha )\cdot {A}_{ij}/{\sum }_{k}{A}_{ik}+\alpha /n$ . The coefficient α > 0 ensures that the constructed Markov chain is irreducible, and the Perron–Frobenius theorem guarantees that P_α has a unique positive left eigenvector π with eigenvalue 1, where π is also called the stationary distribution. The normalized symmetric Laplacian of the Markov chain is ${\mathcal{L}}=I-\frac{1}{2}\left({\Pi }^{1/2}{P}_{\alpha }{\Pi }^{-1/2}+{\Pi }^{-1/2}{P}_{\alpha }^{\top }{\Pi }^{1/2}\right)$ .

To approximately search for the optimal cut, we use the Cheeger inequality for a directed graph³⁹ that bridges the spectral gap of ${\mathcal{L}}$ and the Cheeger constant ϕ^*. As shown in Gleich⁸⁵, the eigenvector v corresponding to the second smallest eigenvalue of ${\mathcal{L}}$ , λ₂, results in optimal clusters. We obtained two clusters by a binary rounding scheme, that is, S₊ = {i ∈ V ∣ v_i ≥ 0} and S₋ = {i ∈ V ∣ v_i < 0}.

We modified the directed_laplacian_matrix function in the NetworkX package (https://networkx.github.io) to calculate the symmetric Laplacian for sparse connectivity matrices, with a default α = 0.05. The spectral gap for the eigenvector centrality subgraph is ${\lambda }_{2}^{{\mathrm{eigen}}}=0.137$ and for the partitioned oculomotor (modO) module is ${\lambda }_{2}^{{\mathrm{eigen}}_{{\mathrm{OM}}}}=0.256$ .

Degree-corrected SBM

Unlike the Louvain and spectral clustering algorithms that assume fewer intercluster connections than intracluster connections, SBMs assume only that there exist blocks such that the probability of two neurons being connected depends only on their block identities and degrees. SBMs are therefore capable of finding structures like cycles that Louvain clustering would not find because SBMs are built on intercluster connections. Starting with the two modules found by applying Louvain clustering to the ‘center’ subgraph, we applied an efficient statistical inference algorithm⁸⁶ to obtain SBMs for each module.

The traditional SBM⁸⁷ is composed of n vertices divided into B blocks with $\left\{{n}_{r}\right\}$ vertices in each block and with probability p_rs that an edge exists from block r to block s. Here, we use an equivalent definition that uses average edge counts from the observation e_rs = n_rn_sp_rs to replace the probability parameters, and we construct a degree-corrected stochastic block model⁸⁸ that further specifies the in- and out-degree sequences $\left\{{k}_{i}^{+},{k}_{i}^{-}\right\}$ of the graph as additional parameters.

To infer the best block membership $\left\{{b}_{i}\right\}$ of the vertices in the observed graph G, we maximize the likelihood ${\mathcal{P}}\left(G| \left\{{b}_{i}\right\}\right)=1/\Omega \left(\left\{{e}_{rs}\right\},\left\{{n}_{r}\right\},\left\{{k}_{i}^{+},{k}_{i}^{-}\right\}\right)$ , where $\Omega \left(\left\{{e}_{rs}\right\},\left\{{n}_{r}\right\},\left\{{k}_{i}^{+},{k}_{i}^{-}\right\}\right)$ is the total number of different graph realizations with the same degree distribution $\left\{{k}_{i}^{+},{k}_{i}^{-}\right\}$ and with $\left\{{e}_{rs}\right\}$ edges among and within blocks of sizes $\left\{{n}_{r}\right\}$ , corresponding to the block membership $\left\{{b}_{i}\right\}$ . Therefore, maximizing likelihood is equivalent to minimizing the microcanonical entropy⁸⁹ ${\mathcal{S}}\left(\left\{{e}_{rs}\right\},\left\{{n}_{r}\right\}\right)=\ln \Omega \left(\left\{{e}_{rs}\right\},\left\{{n}_{r}\right\}\right)$ , which can be calculated as ${\mathcal{S}}\simeq -M-{\sum }_{i}\left({k}_{i}^{+}!\right)-{\sum }_{i}\left({k}_{i}^{-}!\right)-{\sum }_{rs}{e}_{rs}\ln \left(\frac{{e}_{rs}}{{\sum }_{s}{e}_{rs}{\sum }_{r}{e}_{rs}}\right)$ , where M = ∑_rse_rs is the total number of edges. We used the minimize_blockmodel_dl function in the graph-tool package (https://graph-tool.skewed.de) to bisect the central subgraphs by setting B_min = B_max = 2 and degcorr = true.

Shuffled networks

To test the degree to which network properties were determined by specific connections as opposed to statistical connections between blocks, we shuffled synapses while preserving a given block organization. For each cell, we initialized counters to the values of its blockwise in-degree and out-degree. We then randomly reassigned synapses in the following manner to shuffle connections while preserving these degrees. We picked a random presynaptic cell i and a random postsynaptic cell j under the constraint that the counter for the out-degree of cell i to the block containing cell j was greater than 0 and the counter for the in-degree of cell j from the block containing cell i was greater than 0. Each of these counters was then decremented by 1, and the process was repeated until all synapses were reassigned. In this way, once a cell had been selected enough times to reach its original in-degree and out-degree, it could no longer be chosen when assigning the remaining synapses. The shuffling procedure continued until all synapses were reassigned and the blockwise in- and out-degree of each cell was exactly preserved; if this became impossible midshuffle, the procedure was restarted. We considered three different organizations of modO into one, two and seven blocks and calculated each cell’s relative eye position sensitivity and best-fit double exponential decay dynamics as well as the variance explained by the leading eigenvectors and the eigenvalue distribution of the network for 100 random shuffles under each block organization (Extended Data Figs. 8 and 9).

Potential synapse formalism

We define a potential synapse as a presynapse–postsynapse pair within a certain threshold distance (Fig. 2d). This definition is somewhat different from the light microscopic approach, which defines a potential synapse^37,90 as an approach of an axon and dendrite within some threshold distance. Also, we use neurons reconstructed from a single animal, whereas the light microscopic approach aggregates neurons from multiple animals or ‘clones’ a single neuron many times⁹¹.

Calcium imaging and eye position signals

The complete methods for recording calcium activity used to create the functional maps are reported in Ramirez and Aksay⁵⁰. Briefly, we used two-photon, raster-scanning microscopy to image calcium activity from single neurons throughout the hindbrain of 7- to 8-day-old transgenic larvae (sex indeterminate) expressing nuclear-localized GCaMP6f (Tg(HuC:GCaMP6f-H2B) strain cy73-431 from the laboratory of M. Ahrens, HHMI Janelia Research Campus). Horizontal eye movements were recorded simultaneously with calcium signals using a substage CMOS camera. We used CalmAn-MATLAB software to extract the neuronal locations from fluorescence movies⁹².

We analyzed STA activity to determine which neurons were related to eye movements⁵⁰. For each cell, we interpolated fluorescence activity that occurred within 5 s before or after saccades to a grid of equally spaced, 0.33-s time points and then averaged the interpolated activity across saccades to compute the STA. Separate STAs were taken for saccades toward the left and right. We performed a one-way analysis of variance (ANOVA) on each STA to determine which neurons had significant saccade-triggered changes in average activity (P < 0.01 using the Holm–Bonferroni method to correct for multiple comparisons). To determine which of these neurons had activity related to eye position and eye velocity, we first performed a principal components analysis on the STAs from neurons with significant saccade-triggered changes. We found that the first and second principal components had postsaccadic activity largely related to eye movement and eye velocity sensitivity, respectively (see Fig. 3a in Ramirez and Aksay⁵⁰). We characterized each STA using a scalar index, named ϕ in Ramirez and Aksay⁵⁰, created from that STA’s projections onto the first two principal components. We found that this index provides a good characterization of the average eye position- and eye velocity-related activity seen across the population (see Fig. 3c in Ramirez and Aksay⁵⁰ for a map of values and population average STAs). Eye position and eye velocity neurons were defined as neurons with an STA whose value of ϕ was within a specific range (–83 to 37 and 38 to 68, respectively). We removed any neurons with presaccadic activity that was significantly correlated with time until the upcoming saccade. The locations of each neuron were then registered to the Z-Brain atlas³¹ (see Ramirez and Aksay⁵⁰ for complete details).

Categorization of functional cells into cell type was guided by registration to the EM dataset. Candidate functional ABD neurons were identified by proximity (6 μm) to any ABD neuron from the registered EM dataset. Functional Ve² neurons were identified by proximity (6 μm) to any Ve² neuron from the registered EM dataset; because the above restrictions filtered out all of the candidate Ve² neurons, which had very weak sensitivity to spontaneous eye movements, in this case we simply used the constraint that these neurons be deemed active in the CaImAn-MATLAB step. Candidate VPNI cells were determined to be the remaining cells in R4 to R7/8 that satisfied the above restrictions and were more than 30 μm away from the center of the abducens.

Relationship between firing rate and eye position

To assess the functional characteristics of various oculomotor neurons (Fig. 5), we fit a traditional model of eye position sensitivity to neuronal firing rates extracted from our fluorescence measurements. Neuronal firing rates were estimated, up to an unknown scaling factor, by deconvolving the raw calcium fluorescence traces using the deconvolution algorithm with non-negativity constraint described in Pnevmatikakis et al.⁹³. Comparison of the resulting relative firing rates across neurons in different populations was justified as we observed similar sensor expression levels and baseline noise levels in these populations.

We modeled the dependence of firing rate on eye position using the equation $r={[\tilde{k}(E-{E}_{{\mathrm{th}}})]}_{+}$ , where $\tilde{k}$ is the relative eye position sensitivity, E_th is the threshold eye position at which r becomes positive, and the function [x]₊ = x if x > 0 and 0 otherwise⁴⁹. To facilitate comparison of results across animals, we normalized the units of eye position before fitting the model by subtracting the median eye position about the null position (measured as the average raw eye position) and then dividing by the 95th percentile of the resulting positions. Because our focus was on a cell’s position dependence, we also eliminated the eye velocity-dependent burst of spiking activity at the time of the saccade that ABD and VPNI neurons are known to display by removing samples that occur within 1.5 s before or 2 s after each saccade. Saccade times were found in an automated fashion by determining when eye velocity crossed a threshold value⁵⁰. Finally, because the eye position and fluorescence were recorded at different sampling rates, we linearly interpolated the values of neuronal activity at the eye position sample times.

To fit the value of $\tilde{k}$ , for each cell, we defined the eye movements toward the cell’s responsive direction as positive so that $\tilde{k}$ is positive by construction. We then determined the threshold E_th using an iterative estimation technique based on a Taylor series approximation of [x]+ described by Muggeo⁹⁴. Using the resulting estimate of E_th, we determined $\tilde{k}$ as the slope resulting from a linear regression (with offset) to r using E as a regressor. Because we do not know the cell’s responsive direction or preferred eye a priori, we ran the model four times (for each eye, once with movements to the left as positive and once with movements to the right as positive) and used the value of $\tilde{k}$ that resulted in the highest R² value.

As a goodness-of-fit measure, we required all neurons, except Ve² neurons, to have an R² value greater than 0.4. Additionally, non-Ve² neurons were required to have an STA with at least one significant time point (P < 0.01 by ANOVA using Holm–Bonferroni correction) as defined in Ramirez and Aksay⁵⁰ and to have a ΔF/F response that was loosely related to eye position (R² greater than 0.2 when we run the above model replacing r with ΔF/F). The relative eye position sensitivity, $\tilde{k}$ , for fluorescence data was then scaled to the average physiological VPNI response from goldfish⁴⁹.

All functional neurons associated with the circuit schematized in Fig. 5a were ipsiversive in position sensitivity. All ABD neurons increased activity with ipsiversive steps in eye position. The vast majority of putatutive VPNI neurons (98%) also had ipsiversive sensitivity; those few cells with contraversive sensitivity were deemed unlikely to be part of the VPNI based on previous work^49,56 and not considered further. Candidate Ve² neurons were mixed in their sensitivity to the direction of eye movement (52% ipsiversive); both types were included in the distributions of Fig. 5d. However, based on previous work³² demonstrating contrasting activity patterns for cells with ipsilateral versus contralateral projections, it was assumed in the below model that only those with ipsiversive sensitivity had ipsilateral projections as schematized in Fig. 5a.

Identification of excitatory versus inhibitory neurons

Prior work has identified two major classes of VPNI cells, those with ipsilateral recurrent excitatory projections and those with inhibitory contralateral projections^21,49. In R7/R8, the former were found to have somata primarily within a medial glutamatergic parasagittal stripe associated with the Alx transcription factor and the latter with a more lateral GABAergic parasagittal stripe associated with the Abx1b transcription factor^21,95,96. We therefore reasoned that membership in modO, which was predicated to have an ipsilateral recurrent connection with another member of modO, selected for glutamatergic VPNI cells and selected against GABAergic VPNI cells (which were retained in the ‘periphery’). Consistent with this, the VPNI neurons found in the center were located primarily in a medial stripe running from R4 to R7/R8 (Fig. 2c) that overlapped in the Z-Brain atlas with a region of alx expression (ZBB-alx-gal4). ModO neurons at the extreme lateral edge of R5/R6, identified as Ve² neurons (see above), also had ipsilateral projections, but these have been identified previously as being inhibitory through structure/function analysis³².

Network model based on synaptic wiring diagram

A unilateral model of the oculomotor integrator was built using the reconstructed synapses for ABD, Ve² and generic modO neurons that contain the majority of the VPNI population. Although the VPNI is a bilateral circuit, previous experiments^20,56 have shown that one-half of the VPNI is nevertheless capable of maintaining the ipsilateral range of eye positions after the contralateral half of the VPNI is silenced. This may reflect that most neurons in the contralateral half are below a threshold for transmitting synaptic currents to the opposite side when the eyes are at ipsilateral positions⁴³. Therefore, we built a recurrent network model of one-half of the VPNI circuit based on the modO neurons that we had reconstructed from one side of the zebrafish brainstem. We did not include modA neurons, assuming that input from modA cells fell below a threshold needed to drive modO neurons, similar to the manner in which bilateral interactions have been shown to be negligible for the maintenance of persistent activity^20,43. We added the ABD neurons and the feedforward connections from modO to ABD to the model because the ABD neurons are the ‘read-out’ of the oculomotor signals in the VPNI. Projections from the Ve² population were taken to be inhibitory, and all other connections were taken to be excitatory, as explained above. The remaining neurons in modO were included as part of the VPNI; we did not consider other vestibular populations because only Ve² neurons in modO received vestibular afferents from the VIIIth nerve. Saccadic burst neurons and recently discovered presaccadic ramp neurons located in R2/R3 (ref. ⁵⁰) were outside of our reconstructed volume and were therefore not included. Candidate burst neurons in dorsal regions of R5 and R6 have also been recently identified⁹⁷; to examine the possible impact of misclassifying putative burst neurons as VPNI cells, we ran the model with and without these candidates (identified via registration to the functional maps in Ramirez and Aksay⁵⁰).

Directed connection weights between each pair of neurons were set in proportion to the number of synapses from the presynaptic neuron onto the postsynaptic neuron divided by the total number of synapses onto the postsynaptic neuron, ${W}_{ij}=\pm \beta \frac{{N}_{ij}}{{\Sigma }_{k}{N}_{ik}}$ . Thus, we assume that each element W_ij corresponds to the fraction of total inputs to neuron i that are provided by neuron j. We then constructed a linear rate model for the network, governed by $\tau \frac{d{r}_{i}}{dt}=-{r}_{i}+{W}_{ij}{r}_{j}$ , where r_i is the firing rate of the ith neuron, W_ij are the connection weights, and τ is an intrinsic cellular time constant. τ for ABD and Ve² neurons was fixed at 100 ms. The scale factor β and the intrinsic time constant τ for the remaining putative VPNI neurons then provided two free parameters that could be used to set the timescale of persistence for the leading two eigenvectors of the network. For a given choice of scale factor and VPNI neuron time constant, we simulated the response of the network to a random pulse of excitatory input and compared the time course of the leading principal components of the simulated network to that of the calcium imaging data. These parameters were tuned so that the simulations roughly matched the decay times of the leading two principal components in the data, resulting in an intrinsic VPNI neuron time constant of 1 s and a scale factor set such that the leading eigenvalue of the network was equal to 0.9. These free parameters of the network model were not finely tuned, and the precise values used are not critical to the results obtained here.

Relative eye position sensitivities were determined numerically by initializing all neurons to zero firing rate and then simulating the response of the network to a series of three saccades occurring at intervals of 7 s. The input received by the network during a saccade was modeled as a 10-ms pulse along a vector composed of a fixed component that was an equal mixture of the leading three eigenvectors of the weight matrix and a random excitatory vector with an amplitude that was 15% of the amplitude of the fixed component. The random component was assigned at the beginning of the simulation, so the network received the same input for each saccade. The same method used to relate firing rates extracted from calcium imaging data to position sensitivity was used for the simulated firing rates, with the exception that the threshold estimation step was not required for the simulations because all firing rates in the simulations were above threshold. Saccade times in the simulations were taken to be at the onset of each input pulse. Finally, to obtain absolute position sensitivities that could be compared to experiments, we multiplied the relative sensitivities by a single global scale factor that was determined by matching the average VPNI population responses from the model to the average physiological VPNI responses from goldfish⁴⁹.

For both the simulated network and the calcium imaging data, we fit each neuron’s firing rate during the period between 2 s and 6 s following a saccade to a double exponential function, $r(t)={a}_{1}{e}^{-t/{\tau }_{1}}+{a}_{2}{e}^{-t/{\tau }_{2}}$ , where the time constants were fixed at 10 s and 4 s, respectively, leaving the amplitudes of each exponential as the only free parameters. These time constants were chosen because a simulated network with leading eigenvectors having these time constants had principal components of its neuronal firing rates that consisted of mixtures of exponentials whose summed time courses approximately matched those of the leading two principal components of the data. We found that the correlation between these two amplitudes in simulation better matched the correlation seen in the data when the saccadic input consisted of fixed and random components, as described above. The exact number of leading eigenvectors included in the fixed component is not important to achieve this tuning of the correlation between a₁ and a₂. We chose to use the leading three eigenvectors because the eigenvalues associated with them were clearly separated from the central disk of eigenvalues.

As a test of the robustness of our results to possible errors in connectome reconstruction, we generated a connectome that accounted for the estimated false-positive and false-negative rates of synapse detection by our connectome reconstruction procedure. We generated 100 models by randomly varying the identified synapses according to the estimated false-positive and false-negative rates and calculated the connection weights as described above. The eye position sensitivities with this synaptic detection jitter were reported as the average of these 100 models (Extended Data Fig. 7b). We also tested how robust our results were to the cutoff criterion for including neurons in the recurrently connected center by progressively increasing the minimum eigenvector centrality criterion for counting a neuron as belonging to the center as opposed to the periphery. We then plotted how the simulation model results changed as the number of center neurons was decreased and simultaneously reported the resulting number of generic (non-Ve²) modO neurons (Extended Data Fig. 7a). For each of these robustness tests, we also performed a singular value decomposition and double exponential fits on the simulated firing rates to compare to the functional imaging results (Extended Data Fig. 7c,d). To characterize the degradation in model performance when the actual connectome was replaced by connectomes generated by spatial proximity (potential connections), we reran the position sensitivity analysis using potential connectomes defined for connections within 5 μm (Fig. 5e).

Statistics and reproducibility

Here, we briefly summarize factors relevant for rigor and reproducibility that have been addressed in more detail elsewhere in the manuscript.

EM data were acquired from one reimaged animal (Fig. 1). No statistical method was used to predetermine sample size. Reconstructions were seeded by a group of 22 neurons that had previously been identified as VPNI cells through functional imaging²². Reconstructions were semiautomated and validated (and if necessary corrected) by in-house image analysis experts. Synapses were automatically segmented with a convolutional network that was 95% accurate when tested on ground truth data; an error correction step checking for preidentity/postidentity eliminated 2% of the total. During subsequent analysis of connectivity, whenever statistical tests were applied, distributions were assumed to be normal, but this was not formally tested. Randomization was used extensively in determination of modular organization. Investigators were not blinded during these analyses.

Fluorescence microscopy data were previously acquired from 20 animals⁵⁰. No statistical method was used to predetermine sample size. A one-way ANOVA (with multiple comparison correction; normality was assumed) on STA activity was used to identify activity of significance. Investigators were not blinded during these analyses.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability

Raw EM images, segmentation information, neuronal skeletons and identified synapse data are available at https://seung-lab.github.io/zebrafish/home/.

Code availability

WaferMapper is available at https://lichtman.rc.fas.harvard.edu/LGN/WaferMapper.html, Alembic is available at https://github.com/seung-lab/Alembic.git, BigArrays.jl is available at https://github.com/seung-lab/BigArrays.jl with Apache License Version 2.0., ChunkFlow.jl is available at https://github.com/seung-lab/ChunkFlow.jl with Apache License Version 2.0., Watershed is available at https://github.com/seung-lab/Watershed.jl with GNU General Public License v3.0., and Agglomeration is available at https://github.com/seung-lab/Agglomeration with MIT License. Skeletonization, morphology and functions can be found at https://github.com/seung-lab/RealNeuralNetworks.jl with Apache License 2.0. Postreconstruction analysis can be found at https://github.com/ashwinvish/Zfish_recon.git. Modularity analyses and network simulations are available at https://github.com/goldman-lab/Connectome-Model.git.

References

Briggman, K. L., Helmstaedter, M. & Denk, W. Wiring specificity in the direction-selectivity circuit of the retina. Nature 471, 183–188 (2011).
Article CAS PubMed Google Scholar
Anderson, J. R. et al. Exploring the retinal connectome. Mol. Vis. 17, 355–379 (2011).
PubMed PubMed Central Google Scholar
Liu, Z. et al. Central vestibular tuning arises from patterned convergence of otolith afferents. Neuron 108, 748–762 (2020).
Article CAS PubMed PubMed Central Google Scholar
Liu, Z. et al. Organization of the gravity-sensing system in zebrafish. Nat. Commun. 13, 5060 (2022).
Article CAS PubMed PubMed Central Google Scholar
Wanner, A. A. & Friedrich, R. W. Whitening of odor representations by the wiring diagram of the olfactory bulb. Nat. Neurosci. 23, 433–442 (2020).
Article CAS PubMed PubMed Central Google Scholar
Morgan, J. L., Berger, D. R., Wetzel, A. W. & Lichtman, J. W. The fuzzy logic of network connectivity in mouse visual thalamus. Cell 165, 192–206 (2016).
Article CAS PubMed PubMed Central Google Scholar
Bock, D. D. et al. Network anatomy and in vivo physiology of visual cortical neurons. Nature 471, 177–182 (2011).
Article CAS PubMed PubMed Central Google Scholar
Lee, W.-C. A. et al. Anatomy and function of an excitatory network in the visual cortex. Nature 532, 370–374 (2016).
Article CAS PubMed PubMed Central Google Scholar
Motta, A. et al. Dense connectomic reconstruction in layer 4 of the somatosensory cortex. Science 366, eaay3134 (2019).
Article CAS PubMed Google Scholar
Turner-Evans, D. B. et al. The neuroanatomical ultrastructure and function of a biological ring attractor. Neuron 108, 145–163 (2020).
Article Google Scholar
Lyu, C., Abbott, L. & Maimon, G. Building an allocentric travelling direction signal via vector computation. Nature 601, 92–97 (2022).
Article CAS PubMed Google Scholar
Westeinde, E. A. et al. Transforming a head direction signal into a goal-oriented steering command. Nature 26, 819–826 (2024).
Takemura, S.-y et al. A visual motion detection circuit suggested by Drosophila connectomics. Nature 500, 175–181 (2013).
Article CAS PubMed PubMed Central Google Scholar
Kim, J. S. et al. Space-time wiring specificity supports direction selectivity in the retina. Nature 509, 331–336 (2014).
Article CAS PubMed PubMed Central Google Scholar
Seung, H. S. How the brain keeps the eyes still. Proc. Natl. Acad. Sci. USA 93, 13339–13344 (1996).
Article CAS PubMed PubMed Central Google Scholar
Daie, K., Goldman, M. S. & Aksay, E. R. F. Spatial patterns of persistent neural activity vary with the behavioral context of short-term memory. Neuron 85, 847–860 (2015).
Article CAS PubMed PubMed Central Google Scholar
Yoon, K. et al. Specific evidence of low-dimensional continuous attractor dynamics in grid cells. Nat. Neurosci. 16, 1077–1084 (2013).
Article CAS PubMed PubMed Central Google Scholar
Aksay, E., Gamkrelidze, G., Seung, H. S., Baker, R. & Tank, D. W. In vivo intracellular recording and perturbation of persistent activity in a neural integrator. Nat. Neurosci. 4, 184–193 (2001).
Article CAS PubMed Google Scholar
Aksay, E., Baker, R., Seung, H. S. & Tank, D. W. Correlated discharge among cell pairs within the oculomotor horizontal velocity-to-position integrator. J. Neurosci. 23, 10852–10858 (2003).
Article CAS PubMed PubMed Central Google Scholar
Aksay, E. et al. Functional dissection of circuitry in a neural integrator. Nat. Neurosci. 10, 494–504 (2007).
Article CAS PubMed PubMed Central Google Scholar
Lee, M. M., Arrenberg, A. B. & Aksay, E. R. F. A structural and genotypic scaffold underlying temporal integration. J. Neurosci. 35, 7903–7920 (2015).
Article CAS PubMed PubMed Central Google Scholar
Vishwanathan, A. et al. Electron microscopic reconstruction of functionally identified cells in a neural integrator. Curr. Biol. 27, 2137–2147 (2017).
Article CAS PubMed PubMed Central Google Scholar
Allen, W. F. Formatio reticularis and reticulospinal tracts, their visceral functions and possible relationships to tonicity and clonic contractions. J. Wash. Acad. Sci. 22, 490–495 (1932).
Google Scholar
Ramón-Moliner, E. & Nauta, W. The isodendritic core of the brain stem. J. Comp. Neurol. 126, 311–335 (1966).
Article PubMed Google Scholar
Koyama, M., Kinkhabwala, A., Satou, C., Higashijima, S.-i & Fetcho, J. Mapping a sensory–motor network onto a structural and functional ground plan in the hindbrain. Proc. Natl Acad. Sci. USA 108, 1170–1175 (2011).
Article CAS PubMed PubMed Central Google Scholar
Seung, H. S. Reading the book of memory: sparse sampling versus dense mapping of connectomes. Neuron 62, 17–29 (2009).
Article CAS PubMed Google Scholar
Ding, H., Smith, R. G., Poleg-Polsky, A., Diamond, J. S. & Briggman, K. L. Species-specific wiring for direction selectivity in the mammalian retina. Nature 535, 105–110 (2016).
Article CAS PubMed PubMed Central Google Scholar
Schoonheim, P. J., Arrenberg, A. B., Del Bene, F. & Baier, H. Optogenetic localization and genetic perturbation of saccade-generating neurons in zebrafish. J. Neurosci. 30, 7111–7120 (2010).
Article CAS PubMed PubMed Central Google Scholar
Lee, K. et al. Convolutional nets for reconstructing neural circuits from brain images acquired by serial section electron microscopy. Curr. Opin. Neurobiol. 55, 188–198 (2019).
Article CAS PubMed PubMed Central Google Scholar
Turner, N. L. et al. Synaptic partner assignment using attentional voxel association networks. In 2020 IEEE 17th International Symposium on Biomedical Imaging (ISBI) 1–5 (IEEE, 2020).
Randlett, O. et al. Whole-brain activity mapping onto a zebrafish brain atlas. Nat. Methods 12, 1039–1046 (2015).
Article CAS PubMed PubMed Central Google Scholar
Pastor, A. M., Calvo, P. M., de la Cruz, R. R., Baker, R. & Straka, H. Discharge properties of morphologically identified vestibular neurons recorded during horizontal eye movements in the goldfish. J. Neurophysiol. 121, 1865–1878 (2019).
Article CAS PubMed Google Scholar
Gahtan, E., Sankrithi, N., Campos, J. B. & O’Malley, D. M. Evidence for a widespread brain stem escape network in larval zebrafish. J. Neurophysiol. 87, 608–614 (2002).
Article PubMed Google Scholar
Orger, M. B., Kampff, A. R., Severi, K. E., Bollmann, J. H. & Engert, F. Control of visually guided behavior by distinct populations of spinal projection neurons. Nat. Neurosci. 11, 327–333 (2008).
Article CAS PubMed PubMed Central Google Scholar
Huang, K.-H., Ahrens, M. B., Dunn, T. W. & Engert, F. Spinal projection neurons control turning behaviors in zebrafish. Curr. Biol. 23, 1566–1573 (2013).
Article CAS PubMed PubMed Central Google Scholar
Baldwin, C. Y. & Clark, K. B. Design Rules: The Power of Modularity, Vol. 1 (MIT Press, 2000).
Stepanyants, A. & Chklovskii, D. Neurogeometry and potential synaptic connectivity. Trends Neurosci. 28, 387–394 (2005).
Article CAS PubMed Google Scholar
Reichardt, J. & Bornholdt, S. Statistical mechanics of community detection. Phys. Rev. E 74, 016110 (2006).
Article Google Scholar
Chung, F. Laplacians and the cheeger inequality for directed graphs. Ann. Comb. 9, 1–19 (2005).
Article Google Scholar
Girko, V. L. Circular law. Theory Prob. Appl. 29, 694–706 (1985).
Article Google Scholar
Tao, T., Vu, V. & Krishnapur, M. Random matrices: universality of ESDs and the circular law. Ann. Probab. 38, 2023–2065 (2010).
Yang, R. et al. Cyclic structure with cellular precision in a vertebrate sensorimotor neural circuit. Curr. Biol. 33, 2340–2349 (2023).
Article PubMed PubMed Central Google Scholar
Fisher, D., Olasagasti, I., Tank, D. W., Aksay, E. R. F. & Goldman, M. S. A modeling framework for deriving the structural and functional architecture of a short-term memory microcircuit. Neuron 79, 987–1000 (2013).
Article CAS PubMed PubMed Central Google Scholar
Gonçalves, P. J., Arrenberg, A. B., Hablitzel, B., Baier, H. & Machens, C. K. Optogenetic perturbations reveal the dynamics of an oculomotor integrator. Front. Neural Circuits 8, 10 (2014).
Pujala, A. & Koyama, M. Chronology-based architecture of descending circuits that underlie the development of locomotor repertoire after birth. eLife 8, e42135 (2019).
Turrigiano, G. Homeostatic synaptic plasticity: local and global mechanisms for stabilizing neuronal function. Cold Spring Harb. Perspect. Biol. 4, a005736 (2012).
Article PubMed PubMed Central Google Scholar
Cannon, S. C., Robinson, D. A. & Shamma, S. A proposed neural network for the integrator of the oculomotor system. Biol. Cybern. 49, 127–136 (1983).
Article CAS PubMed Google Scholar
McFarland, J. L. & Fuchs, A. F. Discharge patterns in nucleus prepositus hypoglossi and adjacent medial vestibular nucleus during horizontal eye movement in behaving macaques. J. Neurophysiol. 68, 319–332 (1992).
Article CAS PubMed Google Scholar
Aksay, E., Baker, R., Seung, H. S. & Tank, D. W. Anatomy and discharge properties of pre-motor neurons in the goldfish medulla that have eye-position signals during fixations. J. Neurophysiol. 84, 1035–1049 (2000).
Article CAS PubMed Google Scholar
Ramirez, A. D. & Aksay, E. R. Ramp-to-threshold dynamics in a hindbrain population controls the timing of spontaneous saccades. Nat. Commun. 12, 4145 (2021).
Article CAS PubMed PubMed Central Google Scholar
Skaggs, W. E., Knierim, J. J., Kudrimoti, H. S. & McNaughton, B. L. A model of the neural basis of the rat’s sense of direction. Adv. Neural Inf. Process. Syst. 7, 173–180 (1995).
Zhang, K. Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory. J. Neurosci. 16, 2112–2126 (1996).
Article CAS PubMed PubMed Central Google Scholar
Major, G. & Tank, D. Persistent neural activity: prevalence and mechanisms. Curr. Opin. Neurobiol. 14, 675–684 (2004).
Article CAS PubMed Google Scholar
Joshua, M. & Lisberger, S. G. A tale of two species: neural integration in zebrafish and monkeys. Neuroscience 296, 80–91 (2015).
Article CAS PubMed Google Scholar
Miri, A. et al. Spatial gradients and multidimensional dynamics in a neural integrator circuit. Nat. Neurosci. 14, 1150–1159 (2011).
Article CAS PubMed PubMed Central Google Scholar
Debowy, O. & Baker, R. Encoding of eye position in the goldfish horizontal oculomotor neural integrator. J. Neurophysiol. 105, 896–909 (2011).
Article PubMed Google Scholar
Brysch, C., Leyden, C. & Arrenberg, A. B. Functional architecture underlying binocular coordination of eye position and velocity in the larval zebrafish hindbrain. BMC Biol. 17, 110 (2019).
Article PubMed PubMed Central Google Scholar
Feierstein, C. E. et al. Dimensionality reduction reveals separate translation and rotation populations in the zebrafish hindbrain. Curr. Biol. 33, 3911–3925 (2023).
Article CAS PubMed PubMed Central Google Scholar
Robinson, D. The mechanics of human saccadic eye movement. J. Physiol. 174, 245–264 (1964).
Article CAS PubMed PubMed Central Google Scholar
Sklavos, S., Porrill, J., Kaneko, C. R. & Dean, P. Evidence for wide range of time scales in oculomotor plant dynamics: implications for models of eye-movement control. Vision Res. 45, 1525–1542 (2005).
Article PubMed PubMed Central Google Scholar
Hartline, D. K. Pattern generation in the lobster (Panulirus) stomatogastric ganglion. II. Pyloric network simulation. Biol. Cybern. 33, 223–236 (1979).
Article CAS PubMed Google Scholar
Lappalainen, J. K. et al. Connectome-constrained networks predict neural activity across the fly visual system. Nature 634, 1132–1140 (2024).
Article CAS PubMed PubMed Central Google Scholar
Wang, X.-J., Pereira, U., Rosa, M. G. & Kennedy, H. Brain connectomes come of age. Curr. Opin. Neurobiol. 65, 152–161 (2020).
Article CAS PubMed PubMed Central Google Scholar
Bargmann, C. I. & Marder, E. From the connectome to brain function. Nat. Methods 10, 483–490 (2013).
Article CAS PubMed Google Scholar
Hayworth, K. J. et al. Imaging ATUM ultrathin section libraries with WaferMapper: a multi-scale approach to EM reconstruction of neural circuits. Front. Neural Circuits 8, 68 (2014).
Article PubMed PubMed Central Google Scholar
Saalfeld, S., Fetter, R., Cardona, A. & Tomancak, P. Elastic volume reconstruction from series of ultra-thin microscopy sections. Nat. Methods 9, 717–720 (2012).
Article CAS PubMed Google Scholar
Berger, D. R., Seung, H. S. & Lichtman, J. W. VAST (volume annotation and segmentation tool): efficient manual and semi-automatic labeling of large 3D image stacks. Front. Neural Circuits 12, 88 (2018).
Article PubMed PubMed Central Google Scholar
Zung, J., Tartavull, I., Lee, K. & Seung, H. S. An error detection and correction framework for connectomics. In Proc. 31st International Conference on Neural Information Processing Systems (eds. von Luxburg, U., Guyon, I., Bengio, S., Wallach, H. & Fergus, R.) 6821–6832 (Curran Associates, 2017).
Ronneberger, O., Fischer, P. & Brox, T. U-Net: convolutional networks for biomedical image segmentation. In 18th International Conference on Medical Image Computing and Computer-Assisted Intervention (eds. Navab, N., Hornegger, J., Wells, W. M. & Frangi, A. F.) 234–241 (Springer, 2015).
Ioffe, S. & Szegedy, C. Batch normalization: accelerating deep network training by reducing internal covariate shift. In Proc. 32nd International Conference on Machine Learning (eds. Bach, F. & Blei, D.) 448–456 (JMLR, 2015).
Jia, Y. et al. Caffe: convolutional architecture for fast feature embedding. In Proc. 22nd ACM International Conference on Multimedia (eds. Hua, K. A. et al.) 675–678 (Association for Computing Machinery, 2014).
Kingma, D. P. & Ba, J. Adam: a method for stochastic optimization. Preprint at https://doi.org/10.48550/arXiv.1412.6980 (2014).
Lee, K., Zung, J., Li, P., Jain, V. & Seung, H. S. Superhuman accuracy on the SNEMI3D connectomics challenge. Preprint at https://doi.org/10.48550/arXiv.1706.00120 (2017).
Wu, J., Silversmith, W. M., Lee, K. & Seung, H. S. Chunkflow: hybrid cloud processing of large 3D images by convolutional nets. Nat. Methods 18, 328–330 (2021).
Article CAS PubMed Google Scholar
Paszke, A. et al. Automatic differentiation in PyTorch. PyTorch https://pytorch.org/tutorials/beginner/basics/autogradqs_tutorial.html (2017).
Zlateski, A. & Seung, H. S. Image segmentation by size-dependent single linkage clustering of a watershed basin graph. Preprint at https://doi.org/10.48550/arXiv.1505.00249 (2015).
Bae, J. A. et al. Digital museum of retinal ganglion cells with dense anatomy and physiology. Cell 173, 1293–1306 (2018).
Article CAS PubMed PubMed Central Google Scholar
Sato, M., Bitter, I., Bender, M. A., Kaufman, A. E. & Nakajima, M. TEASAR: tree-structure extraction algorithm for accurate and robust skeletons. In Proc. 8th Pacific Conference on Computer Graphics and Applications (eds. Barsky, B. A., Shinagawa, Y. & Wang, W.) 281–449 (IEEE, 2000).
Bogovic, J. A., Hanslovsky, P., Wong, A. & Saalfeld, S. Robust registration of calcium images by learned contrast synthesis. In 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI), 1123–1126 (IEEE, 2016).
Cabrera, B., Torres, B., Pásaro, R., Pastor, A. M. & Delgado-Garcia, J. M. A morphological study of abducens nucleus motoneurons and internuclear neurons in the goldfish (Carassius auratus). Brain Res. Bull. 28, 137–144 (1992).
Article CAS PubMed Google Scholar
Metcalfe, W. K., Mendelson, B. & Kimmel, C. B. Segmental homologies among reticulospinal neurons in the hindbrain of the zebrafish larva. J. Comp. Neurol. 251, 147–159 (1986).
Article CAS PubMed Google Scholar
Newman, M. Networks (Oxford University Press, 2018).
Blondel, V. D., Guillaume, J.-L., Lambiotte, R. & Lefebvre, E. Fast unfolding of communities in large networks. J. Stat. Mech. 2008, P10008 (2008).
Article Google Scholar
Sporns, O. & Betzel, R. F. Modular brain networks. Annu. Rev. Psychol. 67, 613–640 (2016).
Article PubMed Google Scholar
Gleich, D. Hierarchical directed spectral graph partitioning. Information Networks https://www.stat.cmu.edu/~brian/780/week07/Zhang,%20Zhang/Gleich%202005%20-%20hierarchical%20directed%20spectral.pdf (2006).
Peixoto, T. P. Hierarchical block structures and high-resolution model selection in large networks. Phys. Rev. X 4, 011047 (2014).
Google Scholar
Holland, P. W., Laskey, K. B. & Leinhardt, S. Stochastic blockmodels: first steps. Soc. Networks 5, 109–137 (1983).
Article Google Scholar
Karrer, B. & Newman, M. E. J. Stochastic blockmodels and community structure in networks. Phys. Rev. E 83, 016107 (2011).
Article Google Scholar
Bianconi, G. Entropy of network ensembles. Phys. Rev. E 79, 036114 (2009).
Article Google Scholar
Stepanyants, A., Hof, P. R. & Chklovskii, D. B. Geometry and structural plasticity of synaptic connectivity. Neuron 34, 275–288 (2002).
Article CAS PubMed Google Scholar
Markram, H. et al. Reconstruction and simulation of neocortical microcircuitry. Cell 163, 456–492 (2015).
Article CAS PubMed Google Scholar
Giovannucci, A. et al. CaImAn an open source tool for scalable calcium imaging data analysis. eLife 8, e38173 (2019).
Pnevmatikakis, E. A. et al. Simultaneous denoising, deconvolution, and demixing of calcium imaging data. Neuron 89, 285–299 (2016).
Article CAS PubMed PubMed Central Google Scholar
Muggeo, V. M. R. Estimating regression models with unknown break-points. Stat. Med. 22, 3055–3071 (2003).
Article PubMed Google Scholar
Kimura, Y., Okamura, Y. & Higashijima, S.-I. alx, a zebrafish homolog of Chx10, marks ipsilateral descending excitatory interneurons that participate in the regulation of spinal locomotor circuits. J. Neurosci. 26, 5684–5697 (2006).
Article CAS PubMed PubMed Central Google Scholar
Kinkhabwala, A. et al. A structural and functional ground plan for neurons in the hindbrain of zebrafish. Proc. Natl Acad. Sci. USA 108, 1164–1169 (2011).
Article CAS PubMed PubMed Central Google Scholar
Wolf, S. et al. Sensorimotor computation underlying phototaxis in zebrafish. Nat. Commun. 8, 651 (2017).
Article PubMed PubMed Central Google Scholar

Download references

Acknowledgements

We thank M. Moore, K.W., R.W., S.K., S.M., B.S., S.R., A.P., S.W. and D. Visser for manual annotation and proofreading of neuron reconstructions. We thank W. Wong for help with image data transformation for Eyewire and A. Bae for help with skeletonization. We thank M. Tsodyks, D. Kleinfeld, C. Brody, M. Ahrens, M. Koyama, A. Arrenberg, C. Brysch, K. Daie, V. Mishra and C. Chun for their suggestions. We acknowledge support from R01 NS104926 and R01 EY027036 (E.R.F.A., M.S.G., A.V. and H.S.S.), R01 EY021581 and Simons Foundation Global Brain Initiative (E.R.F.A. and M.S.G.), NIH-NINDS Brain initiative award 5U19NS104648 (M.S.G.), K99 EY027017 (A.D.R.), NIH/NCI UH2 CA203710, ARO W911NF-12-1-0594, the Mathers Foundation, Intelligence Advanced Research Projects Activity via Department of Interior/Interior Business Center contract number D16PC0005 and assistance from Google, Amazon and Intel (H.S.S.). The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. The US Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon.

Author information

Jingpeng Wu
Present address: Center for Computational Neuroscience, Flatiron Institute, New York, NY, USA
Alexandro D. Ramirez
Present address: Department of Physiology and Pharmacology, SUNY Downstate Health Sciences University, New York, NY, USA
These authors contributed equally: Alex Sood, Jingpeng Wu, Alexandro D. Ramirez.
These authors jointly supervised this work: H. Sebastian Seung, Mark S. Goldman, Emre R. F. Aksay.

Authors and Affiliations

Princeton Neuroscience Institute, Princeton University, Princeton, NJ, USA
Ashwin Vishwanathan, Jingpeng Wu, Runzhe Yang, Nico Kemnitz, Dodam Ih, Nicholas Turner, Kisuk Lee, Ignacio Tartavull, William M. Silversmith, Chris S. Jordan, Celia David, Doug Bland, Amy Sterling, H. Sebastian Seung, Kyle Wille, Ben Silverman, Ryan Willie, Sarah Morejohn, Selden Koolman, Marissa Sorek, Devon Jones, Amy Sterling, Celia David, Sujata Reddy, Anthony Pelegrino & Sarah Williams
Center for Neuroscience, University of California, Davis, Davis, CA, USA
Alex Sood & Mark S. Goldman
Institute for Computational Biomedicine and the Department of Physiology and Biophysics, Weill Cornell Medical College, New York, NY, USA
Alexandro D. Ramirez & Emre R. F. Aksay
Computer Science Department, Princeton University, Princeton, NJ, USA
Runzhe Yang, Nicholas Turner & H. Sebastian Seung
Brain & Cognitive Sciences Department, Massachusetts Institute of Technology, Cambridge, MA, USA
Kisuk Lee
Department of Neurobiology, Physiology and Behavior, University of California, Davis, Davis, CA, USA
Mark S. Goldman
Department of Ophthalmology and Vision Science, University of California, Davis, Davis, CA, USA
Mark S. Goldman
Eyewire, Boston, MA, USA
Marissa Sorek, Devon Jones, Amy Sterling & Celia David

Consortia

the Eyewirers

Kyle Wille
, Ben Silverman
, Ryan Willie
, Sarah Morejohn
, Selden Koolman
, Marissa Sorek
, Doug Bland
, Devon Jones
, Amy Sterling
, Celia David
, Sujata Reddy
, Anthony Pelegrino
& Sarah Williams

Contributions

A.V. designed and conceptualized the study, collected and interpreted the data, acquired the EM data, performed registration and data analysis and wrote the paper. A. Sood designed and conceptualized the computational model, performed data analysis and interpretation and wrote the paper. J.W. developed the code for data generation and curation in Eyewire, meshing, skeletonization and convnet inference. A.D.R. collected, registered and analyzed the light microscopic data and interpreted the data. R.Y. performed the clustering algorithm comparisons. N.K. performed the Eyewire data assembly. D.I. performed EM image assembly. N.T. performed synapse detection and partner assignment. K.L. performed boundary detection. I.T. developed Eyewire algorithms. W.M.S. developed Eyewire algorithms and data manipulation software. C.S.J. developed Eyewire algorithms and performed Eyewire system administration. C.D. and D.B. performed Eyewire moderation and data curation. A. Sterling managed Eyewire. H.S.S. designed and conceptualized the study, performed data analysis and interpretation and wrote the paper. M.S.G. designed and conceptualized the computational model, performed data analysis and interpretation and wrote the paper. E.R.F.A. designed and conceptualized the study, interpreted the data and wrote the paper. Eyewirers performed neuron reconstruction online.

Corresponding authors

Correspondence to Ashwin Vishwanathan, Mark S. Goldman or Emre R. F. Aksay.

Ethics declarations

Competing interests

H.S.S., N.K., D.I., N.T. and K.L. have financial interests in Zetta AI, LLC. The remaining authors declare no competing interests.

Peer review

Peer review information

Nature Neuroscience thanks the anonymous reviewers for their contribution to the peer review of this work.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data

Extended Data Fig. 1 Composition of Axial module (modA).

a. Connectivity matrix (as in Fig. 2) with identification (arrows) of small RS neurons that were part of the center and inclusion into the periphery of the remaining small RS neurons along with the large RS neurons. b. Visualization of the large and small classes of RS neurons in the periphery (ro - rostral; c - caudal; m - medial; l - lateral). c. Visualization of small vSPNs that are part of the ‘center’ in modA and are indicated by arrows above the rows in a.

Extended Data Fig. 2 Module-specific synapse size and location distributions.

a. Histogram of neuronal pathlength in modA (n = 251) and modO (n = 289) for all cells (left), those with somata in the reconstructed volume (middle), and those with somata outside the reconstructed volume (right). Here and below, p is the significance value based on a two-sided Wilcoxon-rank sum test (RS-test). b. Histogram of synaptic size (detected PSD voxels [vx]) for connections within (left) and between (right) modules modA and modO. Table (here and below) summarizes the mean and standard deviations. c. Histogram of synapse locations, that is, distance from somata to synaptic site along the neurite, for neurons in modA and modO. Cohen’s D measures the effect size. d. Histogram of synapse sizes (detected PSD) of neurons within-module (left) and between-modules (right) for modO_I and modO_M. e. Histogram of synapse locations for neurons within the oculomotor modules modO_I and modO_M.

Extended Data Fig. 3 Clustering algorithm comparisons.

a. Connectivity when the center is organized by spectral clustering into two modules. b. Normalized number of synapses within and between relevant cell groups.

Extended Data Fig. 4 Clustering of center subgraph is insensitive to selection criteria and synapse detection errors.

a. Connectome for modO and its connections to the abducens with the eigencentrality threshold for inclusion in the center raised to exclude 25% of the cells that were previously included in the center. Cells are ordered by block identity after fitting a stochastic block model with 8 blocks. In this case, an eighth block was needed to sequester the vestibular cells into their own block. Colored rectangles highlight two cycles in the block structure, one connecting primarily to ABD_M (orange) and the other primarily to ABD_I(purple). b. Association matrix for block assignments when segmentation of the center is repeated 100 times after random addition and deletion of synapses according to the estimated precision and recall of synapse detection (see Methods). On average, a random variation of the connectome caused 0.5% of neurons to change block assignments.

Extended Data Fig. 5 Substructure of the center network.

Each dot indicates the presence of one or more synapses between two cells in the center network or from a cell in the center to a cell in the periphery. Cells are ordered by their block identity, determined by using Louvain clustering to split the center into an axial module (modA) and an oculomotor module (modO), followed by Stochastic Block Modeling to find substructure in each module. In the above, modA is split into 3 blocks, and modO is split into 7 blocks.

Extended Data Fig. 6 Example neuron firing rates in the simulated network.

Simulated firing rates for three neurons in response to a sequence of three simulated saccades at intervals of seven seconds (see Methods). Each trace is normalized to the peak firing rate reached by that neuron during the first fixation.

Extended Data Fig. 7 Behavior of the connectome-based model is insensitive to centrality threshold and errors in synapse detection.

a. Relative position sensitivites for different cell classes as the eigencentrality threshold for inclusion in the center is varied. Box bounds indicate 1st and 3rd quartiles; the whiskers extend from the box to the farthest data point lying within 1.5x the inter-quartile range (IQR) from the box. b. A comparison of relative eye position sensitivities between the actual connectome (gray) and 100 variations of the network produced by random addition and deletion of synapses according to the estimated error in synapse detection (green, see Methods). c. Comparison of estimated firing rates of identified integrator cells in larval zebrafish using calcium imaging (green) to firing rates simulated using a connectome-based model when the eigencentrality threshold is raised to remove 25% of the center (gray). (top) Cumulative variance explained for the leading principal components of the average activity between two and six seconds after a saccade. (bottom) Best-fit amplitudes determined as in Fig. 5g. Error bars on the cumulative variance explained (mean) show the standard deviation across the 100 random variations of the network. d. Same as c, but for the variations of the network used in b.

Extended Data Fig. 8 Block structure sufficient to maintain circuit function.

a. Position sensitivities (top, middle) and eigenvalues (bottom) for 100 networks in which all synapses within modO were shuffled while maintaining the in-degree and out-degree of each cell. In the top row, the position sensitivities of all shuffled networks (green) are compared to the actual network (gray) for each cell class (VPNI, Ve², ABD). In the middle row, shuffled networks (x-axis) are compared to the actual network (y-axis) on a cell-by-cell basis, with the colors indicating which block each cell belonged to and opacity indicating the relative frequency of the eye position sensitivities produced for each cell across all of the shuffles. b. Same as a except that modO was partitioned into two blocks using a stochastic block model (SBM) and each group of synapses within each block and between each pair of blocks was shuffled separately. c. Same as b, but for 7 blocks instead of two.

Extended Data Fig. 9 Sub-blocks within ModO are necessary to reproduce observed dynamics.

a. Comparison of estimated firing rates of identified integrator cells in larval zebrafish using calcium imaging (green) to firing rates simulated using a connectome-based model after shuffling the synapses within ModO (gray, see Methods). (left) Cumulative variance explained for the leading principal components of the average activity between two and six seconds after a saccade. Error bars about the mean show the standard deviation across 100 random shuffles of the network. (right) Best-fit amplitudes determined as in Fig. 5g. b, c. Same as a, but with ModO split into two sub-blocks (b) or seven sub-blocks (c) using a stochastic block model. Synapses were shuffled separately for each possible pairing of presynaptic block and postsynaptic block.

Extended Data Fig. 10 Complex modes of cycles in the VPNI exhibit overdamped oscillation.

Stimulating a complex mode of a network gives rise to a linear combination of two spatial patterns of activity, one corresponding to the real part of the eigenvector and one corresponding to the imaginary part. The level of activity in each pattern oscillates over time while the overall amplitude decays, so that each neuron follows a trajectory of the form Ae^−ηt sin (ωt + φ). If the decay is fast enough relative to the frequency of oscillation, as occurs for the complex modes of the internuclear and motor cycles, oscillations will not be evident in the single-neuron firing rates. The curves above show, for the complex mode of the internuclear cycle, a range of potential single-neuron firing rate trajectories corresponding to different phase offsets, from zero initial activity (purple) to maximal initial activity (blue).

Supplementary information

Supplementary Information

Supplementary Figs. 1–4 and Note.

Reporting Summary

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Cite this article

Vishwanathan, A., Sood, A., Wu, J. et al. Predicting modular functions and neural coding of behavior from a synaptic wiring diagram. Nat Neurosci 27, 2443–2454 (2024). https://doi.org/10.1038/s41593-024-01784-3

Download citation

Received: 12 December 2022
Accepted: 11 September 2024
Published: 22 November 2024
Issue Date: December 2024
DOI: https://doi.org/10.1038/s41593-024-01784-3

Subjects

This article is cited by

Wiring of a low-dimensional integrator network
- Bo Hu
- Rainer W. Friedrich
Nature Neuroscience (2024)