Modularity and neural coding from a brainstem synaptic wiring diagram

Image acquisition and alignment

We acquired a dataset of the larval zebrafish hindbrain that extended 250 µm rostrocaudally and includes rhom-bomeres 4 through 7/8 (R4 to R7/8). The volume extends 120 µm laterally from the midline and 80 µm ventrally from the plane of the Mauthner cell axon. The ssEM dataset was an extension of the original dataset in ref ²⁹ and was extended by additional imaging of the same serial sections. Only a few tens of neurons had been manually reconstructed in our original publication on the ssEM dataset²⁹. The dataset was stitched and aligned using a custom package, Alembic (see Code availability). The tiles from each section were first montaged in 2D, and then registered and aligned in 3D as whole sections. Point correspondences were generated by block matching via normalized cross-correlations both between tiles and across sections. The final set of parameters that were used are listed in table.

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 band-pass filter. Stitching of tiles (montaging) within a single section was split into a linear translation step (premontage) and a non-linear 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 single 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 for parameters used).

Once the correspondences were found, outliers were filtered by checking the spatial distribution of the cross-correlogram (sigma filter), height of the peak of the correlogram (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 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 inter-section 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, which were solved for a least-squared-residual affine transform, regularized with 10% (empirically derived) of identity transformation to reduce shear from propagating across multiple sections. Proceeding sequentially allowed the entire stack to get roughly in place. The mean residual errors were in the range of 3.5 pixels after relaxation.

Convolutional Net Training

LM/EM registration

Registration of the EM dataset to the reference atlas (ZBrain, ⁴⁵) was carried out in two stages. We created an intermediate EM stack from the low resolution (270 nm/pixel) EM images of the entire larval brain tissue. This intermediate stack had the advantage of a similar field of view as compared to the LM 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 by fitting an affine transform that maps the entire EM volume onto the LM 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 ZBrain atlas, was used as the template to register the high-resolution EM stack onto it. This was performed in a similar manner by selecting corresponding points and fitting an affine transform. The resulting transform would transform points from the high-resolution EM space to the reference atlas space. This transform was used to map the reconstructed skeletons from high-resolution EM space to the reference atlas space.

Functional imaging and functional maps

The complete methods for recording calcium activity used to create the functional maps are reported in⁴⁴. Briefly, we used two-photon, raster-scanning microscopy to image calcium activity from single neurons throughout the hindbrain of 7-8 day old transgenic larvae expressing nuclear-localized GCaMP6f, Tg(HuC:GCaMP6f-H2B) strain cy73-431 from Misha Ahren’s lab. Horizontal eye movements using a sub-stage CMOS camera were recorded simultaneously with calcium signals. We used the CalmAn-Matlab software to extract the neuronal locations from fluorescence movies⁸⁵.

We analyzed saccade-triggered average (STA) activity to determine which neurons were related to eye movements (see⁴⁴ for complete details). For each cell, we interpolated fluorescence activity that occurred within five seconds before or after saccades to the left (right) to a grid of equally spaced, ⅓ second timepoints and then averaged the interpolated activity across saccades to compute the STA. The STA could occur around saccades to the left or to the right. We performed a one-way 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 (PCA) on the STAs from neurons with significant saccade-triggered changes. We found that the first and second principal components had post-saccadic activity largely related to position and velocity sensitivity respectively (see Figure 3A in ⁴⁴). We characterized each STA using a scalar index, called ϕ in⁴⁴, created from that STA’s projections onto the first two principal components and found that this index does a good job of characterizing the average position and velocity-related activity seen across the population (see Figure 3C in⁴⁴ for a map of ϕ values and population average STAs). Position and velocity neurons were defined as neurons with an STA whose ϕ value of was within a specific range (−83 to 37 and 38-68 respectively). We removed any neurons with pre-saccadic activity that was significantly correlated with time until the upcoming saccade. The locations of each neuron were then registered to the Z-Brain atlas using similar methods as listed in the previous section (see⁴⁴ for complete details).

Computing slopes of eye-position dependence

To assess the functional characteristics of various oculomotor cells (Figure 4) we fit a classical model of eye position sensitivity to neuronal firing rates extracted from our fluorescence measurements. We approximated the firing rate, r, of a cell using the deconvolution algorithm with non-negativity constraint described in⁸⁶. We modeled the dependence of firing rate on eye position using the equation, r = [k(E − E_th)]₊, where k and E_th are free parameters and the function [x]₊ = x if x > 0 and 0 otherwise⁸⁷. In order to best compare 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. Since our focus was on a cell’s position-dependence, we also removed the eye velocity dependent burst of spiking activity at the saccade time that Abducens and VPNI neurons are known to display by removing samples that occur within 1.5 seconds before or 2 seconds after each saccade. Saccade times were determined in an automated fashion by determining when eye velocity crossed a threshold value⁴⁴. Finally, since 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 k, for each cell, we defined the eye movements toward the cell’s responsive direction as positive so that k is positive by construction. We then determined the ratio of the threshold to eye position, E_th/K, using an iterative estimation technique based on a Taylor series approximation of [x]₊ described in ref ⁸⁸. Using the resulting estimate of E_th/K, we determined k as the slope resulting from a linear regression, with offset, to r using the regressor [E = E_th/K]₊. Since we did not know the cell’s responsive direction a priori, we ran the model twice -- once with movements to the left as positive and once with movements to the right as positive -- and used the value of k that resulted in the highest R² value. We fit the model using positions from the contralateral eye for the Abducens internuclear neurons and the ipsilateral eye for all other neurons. A neuron was identified as an Abducens (ABD), Abducens internuclear (ABD_I), or descending octavolateral (DO) neuron if its location on the ZBrain atlas was within 5 microns of an ABD_M, ABD_I, or DO cell determined from EM anatomy (since EM reconstruction was only performed for one hemisphere, we duplicated the locations of the EM-identified neurons and mirrored their location onto the opposite hemisphere). As a goodness-of-fit measure we required all neurons, except DO cells, to have an R² value greater than 0.4. Additionally, non-DO cells were required to have a saccade-triggered average with at least one significant time point (p<0.01 by an ANOVA test using Holm-Bonferroni correction) as defined in ref ⁴⁴ and to have a (which is defined as ), fluorescence response that was loosely related to eye position (R² greater than 0.2 when we run the above model replacing r with ). The eye position sensitivity for fluorescence data was then scaled to average physiological VPNI responses from goldfish ^48:89.

Modular analysis

We organized the reconstructed neurons into two groups: those in a recurrently connected “center”, where feedback interactions are expected to establish collective dynamics, and those in the “periphery” that project to or collect input from the center. To do so, we first defined a connection matrix N_ij by the number of synapses onto node i from node j. We computed the left and right eigenvectors of the connection matrix corresponding to the eigenvalue with maximal real part. The elements of these eigenvectors, which are all non-negative, are regarded as measures of network “centrality”. We define eigen-centrality as the geometric mean of the left and right eigenvectors corresponding to the maximal real part. This terminology is helpful because nodes with zero out-degree (orphan dendrites) have zero out-centrality, and nodes with zero in-degree (orphan axons) have zero in-centrality. We note for completeness that our definition of centrality is properly termed an eigen-centrality; similar organization could be generated by using the related concept of degree-centrality. Degree-centrality (geometric mean of in-degree and out-degree) and eigen-centrality are correlated, but not perfectly (Extended Data Fig. 1c).

Using the above method, we were able to define a ‘center’ that contained 540 recurrently connected neurons (eigen-centrality >10⁻⁸). We then applied three graph-clustering algorithms, the Louvain algorithm, Spectral clustering and the Stochastic block matching method to divide the center into modules. The Louvain based methods is presented in the main text, and the rest are described in the Extended data.

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 . Here s_i = ∑_c W_ic is the sum of weights out of node i, s_j = ∑_c W_jb is the sum of weights out of node j, ω = ∑_cd W_cd is the total sum of weights in the network, and the resolution parameter γ determines how much the naively expected weight of connections γ ^} is subtracted from the connectivity matrix. Potential degeneracy in graph clustering was addressed by computing the consensus of the clustering similar to ⁹⁰. Briefly, an association matrix, counting the number of times a node (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 node. Reclustering the thresholded association matrix, where threshold was the maximum element of the randomized association matrix, provided consensus modules. We used the commuinity_louvain.m function from the Brain Connectivity Toolbox package (BCT, https://sites.google.com/site/bctnet/Home). In addition to the Louvain graph-clustering algorithm, we also clustered the ‘center’ with two alternate graph-clustering algorithms; spectral clustering and stochastic block matching, described below.

Spectral Clustering

We employed a generalized spectral clustering algorithm for weighted directed graphs to bisect the zebrafish ‘center’ subgraph as proposed by³⁸. Given a graph G(V, E) and its weighted adjacency matrix , A_ij where 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_α]_ij := (1 − α) · A_ij/ ∑_k A_ik + α/n. The coefficient α > 0 ensures that the constructed Markov chain is irreducible, and Perron-Frobenius theorem guarantees 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 .

To approximately search for the optimal cut, we utilize the Cheeger inequality for a directed graph³⁸ that bridges the spectral gap of ℒ and the Cheeger constant ϕ^∗. Ref⁹¹ shows that the eigenvector v corresponding to the second smallest eigenvalue of ℒ, λ₂ results in optimal clusters. We obtained two clusters by a binary rounding scheme, i.e., 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 gaps for the eigenvector-centrality subgraph is and for the partitioned the oculomotor (mod_O) module is .

Degree-corrected Stochastic Block Matching (SBM)

Unlike the Louvain and spectral clustering algorithms that assume fewer intra-cluster connections than inter-cluster connections, we applied an efficient statistical inference algorithm ³⁹ to obtain the stochastic block models (SBM) that best describes the ‘center’ subgraph.

The traditional SBM⁹² is composed of n vertices, divided into B blocks with {n_r} vertices in each block, and with the probability, p_rs, that an edge exists from block r to block s. Here we use another equivalent definition, to use average edge counts from the observation e_rs = n_rn_sp_rs to replace the probability parameters. The degree-corrected stochastic block model⁹³ further specifies the in- and out-degree sequences of the graph as additional parameters.

To infer the best block membership {b_i} of the vertices in the observed graph G, we maximize likelihood Р(G|{b_i}) = 1/Ω({e_rs}, {n_r}), where Ω({e_rs}, {n_r}) is the total number of different graph realizations with the same degree distribution , and {e_rs} edges among and within blocks of sizes {n_r}, corresponding to the block membership {b_i}. Therefore, maximizing likelihood is equivalent to minimizing the microcanonical entropy⁹⁴ S ({e_rs}, {n_r}) = ln Ω({e_rs}, {n_r}), which can be calculated as , where M =∑_rs e_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.

Network level modeling

A full-scale, one-sided model was built using the reconstructed synapses for the ABD, DO, and VPNI populations. A full-scale, one-sided model was built using the reconstructed synapses for the ABD, DO, and VPNI populations. Although the VPNI is a bilateral circuit, previous experiments^87;89 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 their thresholds 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. For simplicity, we did not include the modA neurons, assuming that the weak connectivity between the axial and oculomotor modules makes their contributions negligible, similar to the manner in which bilateral interactions have been shown to be negligible for the maintenance of persistent activity. 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.

The elements of the model weight matrix are . An overall scale factor β was applied to tune the principal eigenvalue to one, which is a necessary condition for generating stable temporal integration in the network model. The automated synapses detection identified a small number of synapses (8 synapses from 3 ABD neurons) projecting from the abducens onto integrator and vestibular neurons. Manual inspection of the synapses indicated that they were false positives, so they were neglected when building the model. Projections from the DO population were taken to be inhibitory while all other connections were taken to be excitatory (supplementary information). The remaining neurons in the oculomotor module were included as part of the VPNI--we did not consider other vestibular populations since only DO neurons in modO received vestibular afferents from the Viiith nerve; saccadic burst neurons are likely located in R2/3⁴⁴, outside of our reconstructed volume; and recently discovered pre-saccadic ramp neurons⁴⁴, are quite sparse in comparison to neurons with position signal.

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 synapses onto the post-synaptic neuron. The weights were then scaled to achieve perfect integration in a linear rate model governed by where r_i is the firing rate of the i^th neuron, W_ij are the connection weights, and τ, the intrinsic time constant, is 100ms. Position sensitivities were determined for the neurons in each model by simulating the response of the network to a pulse of input along the integrating direction. The position sensitivity for each neuron was then determined by averaging across models and scaled to average physiological VPNI responses from goldfish^48;89.

As a control, we generated a connectome that accounted for the false positive and negative rate of synapse detection by the neural nets; synaptic detection jitter. We generated 1000 models by randomly varying the identified synapses according to the false positive and false negative rates and calculated the connection weights as described above. The eye position sensitivities with synaptic detection jitter were reported as the average of these 1000 models (Extended Data Fig. 10c).

Abducens neuron identification

We identified abducens motor neurons (Fig. 1c, Extended Data Fig. 3) by their overlaps with abducens motor neurons (ABD_M) labeled by the mnx transgenic line (Extended Data Fig. 3a). The ABD_M axons exited R5 and R6 through the abducens (VIth) nerve (Fig. 1c, black box) as reported previously [1].

Contraversive horizontal movements of the eye are driven by the medial rectus muscle, which are innervated by motor neurons in the oculomotor nucleus, which in turn are driven by internuclear neurons (ABD_I) in the contralateral abducens complex (Fig. 1c). ABD_I neurons were identified (Extended Data Fig. 3a) by their overlaps with two nuclei in the evx2 transgenic line [2] that labels glutamatergic interneurons. The ABD_I neurons were just dorsal and caudal to the ABD_M neurons, and their axons crossed the midline [3].

Identification of DO neurons

We identified a class of secondary vestibular neurons; Descending Octavolateral (DO) neurons (Fig. 1d, brown). We observed that DO cells received synapses from primary vestibular afferents. The latter were orphan axons in R4 identified as the vestibular branch of the vestibulocochlear nerve (VIIIn) by comparison with the isl-2 line, which labels the major cranial nerves (Extended Data Fig. 3b, blue axons).

Identification of vSPNs

Ventromedially located Spinal Projection Neurons (vSPNs) were identified by their stereotypic locations [4] and by comparing the registered ssEM volume to spinal backfilled neurons part of the Zbrain resource (Extended Data Fig. 4). The vSPNs were divided into two groups, large and small vSPNs (Fig. 1d). Large vSPNs were M, Mid2, MiM1, Mid3i and CaD neurons. Small vSPNs were RoV3, MiV1, MiV2.

Assigning sign to VPNI neurons

We first attempted to determine the signs of the neural connections by comparison with molecular maps in the Z-Brain atlas. Integrator neurons in the dorsomedial stripe running from R4 to R7/8 overlap with a region of alx expression (Extended Data Fig. 9c, ZBB-alx-gal4) where most neurons are glutamatergic [5]. More lateral and caudal neurons overlap with markers for glycine and glutamate. Finally, neurons in the oculomotor module do not overlap with the GABA expression (Extended Data Fig. 9c, ZBB-gad1b). For these reasons were assigned neurons in the oculomotor module to be excitatory.

Gamification of the zebrafish dataset and crowdsourced reconstructions

In addition to in-house reconstructions by experts, zebrafish cells were made available to citizen scientist gamers through the existing retinal crowdsourcing game “Eyewire.”

To test out this new dataset only the most experienced players on Eyewire were allowed to participate. A small group of 4 highly experienced players received an invite 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 pre-recorded 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 GUI was built into the game that allowed players to see the status of each cell currently online. A cell could be at one of 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 done 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 to 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.

The following players in Eyewire were the admins who validated crowd sourced neuronal reconstructions:

Hoodwinked, BenSilverman, Hightower, sarah.morejohn, SeldenK, sorek.m, zorek.m, twisterZ, hjones.jr, devonjones, amy, EinsteintheRapper, zkem, celiad, celiaz,sunreddy, peleaj43, sarahaw

The following players helped reconstruct neurons on Eyewire, collectively called - Mystic players:

r3, Atani, Nseraf, susi, eldendaf, Frosty, a5hm0r, hiigaran, kinryuu, Manni_Mammut, aesanta1, LotteryDiscountz, galarun, annkri, dragonturtle, LynneC, Cliodhna, jax123, KrzysztofKruk, Kfay, rinda, crazyman4865, JousterL, randompersonjci, Caffeine, Baraka, ggreminder, TR77, hewhoamareismyself, nagilooh, Oppen_heimer, Gruenewitwe, cognaso, twotwos, hawaiisunfun, danielag, lemongrab, zope, MysticM, kondor, frankenmsty, zfishman.