A Brain-wide Circuit Model of Heat-Evoked Swimming Behavior in Larval Zebrafish.

CONTACT FOR REAGENT AND RESOURCE SHARING

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Florian Engert ( florian@mcb.harvard.edu ).

EXPERIMENTAL MODEL AND SUBJECT DETAILS

All experiments were conducted on 6-7 days post fertilization zebrafish homozygous for the nacre mutation expressing the transgenes indicated below. The sex of the larva is not defined at this early stage. Fish were fed paramecia from day 5 onwards. All experiments followed the guidelines of the National Institutes of Health and were approved by the Standing Committee on the Use of Animals in Research of Harvard University. All experiments used for mapping heat responsive neurons and for model derivation used fish expressing nuclear Elavl3:H2B-GCaMP6s fish ( Freeman et al., 2014 ). Experiments combining heat and taps as well as trigeminal ablations were performed in fish expressing cytoplasmic Elavl3:Gcamp6s ( Kim et al., 2017 ). We note that while the decay time-constants for cytoplasmic and nuclear GCaMP are different, we found this effect to be negligible given the slow timescales of our temperature stimulus. To develop the scan stabilization protocol Elavl3:H2B-RFP fish expressing RFP in all neuronal nuclei were used ( Randlett et al., 2015 ). To identify neurotransmitter types of heat modulated neurons, progeny of crosses between Elavl3:H2B-Gcamp6s ( Freeman et al., 2014 ) and either vglut2a:mCherry ( Satou et al., 2013 ) or gad1B:dsRed ( Satou et al., 2013 ) expressing fish were used. METHOD DETAILS Since almost all analysis in this study was performed in an automated manner, no blinding or randomization was performed. All sample sized were fixed before the start of experiments and all animals were analyzed except those that freed themselves from the preparation during imaging or that died during the experimental procedure.

Imaging and Behavior

Larval zebrafish were embedded in 2.5% medium melt agarose (Fisher scientific, USA) and their tails were freed the night before the experiment. Experiments were conducted in a custom built 2-photon microscope and run using custom written software in C# (Microsoft, USA). Heat stimuli were delivered using a 1W 980nm fiber-coupled diode laser (Roithner, Austria) coupled into a collimator (Aistana Inc., USA) placed under the microscope objective 4 mm in front of and 1.2 mm above the head of the zebrafish larva pointing downwards at an angle of 16.5 degrees. The laser power was controlled by the computer via a laser diode driver (Thorlabs, USA). We note that the 980 nm laser itself did not excite GCaMP fluorescence due to the low photon density. The main mapping experiments consisted of the imaging of 30 individual planes, spaced 2.5 μm apart. In each plane 3 trials of the stimulus depicted in Figure 1B were presented. The heat and tap experiments consisted of imaging 4 individual planes, spaced 5 μm apart. In each plane 25 trials of the stimulus depicted in Figure S1G were presented. To avoid excessive heating of the preparation by scanning over the eyes, custom exclusion masks were created for each experiment in which the eyes were in the field of view, restricting the scan-lines such that the eyes were excluded from the field of view. Imaging was performed at 2.4 Hz - 3 Hz and all imaging data was interpolated to a 5 Hz timebase before further analysis. Tail-tracking data was acquired at 100 Hz. All behavioral features such as identifying and classifying bouts were performed at this timebase, however, for all comparisons of behavior and imaging, the behavioral data was downsampled to a 5 Hz timebase.

Stimulus characterization

To measure the temperature induced in larval zebrafish by the laser stimulus a thermistor with the same dimension and absorption characteristics as larval zebrafish was used (Warner instruments, USA; ( Haesemeyer et al., 2015 )). The thermistor was embedded in the same low-melting point agarose as the larval zebrafish and the laser was positioned accordingly ( Madelaine et al., 2017 ). For each stimulus used, temperature changes in the thermistor were subsequently recorded. While tissue absorption is not homogeneous for larval zebrafish, modeling of thermal flux revealed that over the short length scales at this age there won’t be an uneven distribution of temperature within the head of larval zebrafish due to the thermal conductivity of water within about 250 ms of stimulus onset. In other words, equilibration happens much faster than the timescales of our stimulus (not shown). Furthermore, the size of the laser spot was set such that it is larger than the head of the fish. Image stabilization To counter heat-induced deformations of the preparation, before each plane was scanned a +/- 5 μm sized pre-stack consisting of 21 slices spaced 0.4 μm apart was acquired. During scanning, each acquired plane was cross-correlated with each plane in the pre-stack and using a low-pass filter the position of the objective was adjusted online so as to minimize z-drift. Since the heat induced drift observed in Elavl3:H2B-RFP stacks followed very reproducible kinetics these were used to predict the movement induced by heating during our experiments to induce very slight movements in the predicted direction. This served to overcome the delay induced by the low-pass filter. A second alignment step was performed post-acquisition. Here each individual imaging plane was assigned to the most likely plane in the pre-stack that it corresponds to via correlation. The fluorescence timeseries of each segmented region was subsequently corrected by normalizing it with the resting fluorescence of this region in the pre-stack.

Image segmentation

Cells in each plane were segmented anatomically. To correct for motion artifacts individual planes in each timeseries were re-aligned based on image cross-correlations ( Miri et al., 2011 ). For nuclear GCaMP experiments, individual nuclei were subsequently segmented using Cell Profiler ( Carpenter et al., 2006 ). To resolve merged nuclei in areas of low contrast, objects larger than a typical nuclear size were divided into subregions based on pixel timeseries correlations. For experiments using cytoplasmic GCaMP, custom written software was used to identify individual cells. Each image was processed using a minimum and a maximum filter tuned to nuclear size. Points where the differences between the two filtered results crossed a threshold (at least 20 % of the average region brightness) were considered potential cell centroids. After including seed pixels around the centroids individual regions were grown in a greedy manner, incorporating pixels that were better correlated to their own current average activity than to neighboring averages. Resulting correlations masks were intersected with anatomical masks based on cell size to obtain the final segmentation. Registration and annotation 3D image registration based on CMTK ( Rohlfing and Maurer, 2003 ) was used to create a nuclear GCaMP-6s reference stack to which all experimental stacks were registered as described elsewhere ( Portugues et al., 2014 ; Randlett et al., 2015 ). The reference brain was used to annotate anatomical regions of interest based on Z-Brain annotations ( Randlett et al., 2015 ).

Swim bout identification and classification

Swim bouts were identified as described previously ( Portugues et al., 2015 ), namely based on the windowed standard deviation of the tail cumulative angle trace crossing a threshold. The mode-centered cumulative angle trace a of each bout was subsequently used to assign the following score: bias = ∑ a a ∑ a | a | Right flicks were defined as bouts with a bias < -0.8, left flicks as bouts with a bias > 0.8 and all other bouts were defined as “swims”. These cutoffs were chosen since the histogram of all bias scores has a minimum at those points. Clustering of heat and motor related activity To identify motor related activity, motor regressors were created by convolving a bout start trace with an exponential calcium kernel with a decay half time of 3s. This decay time was derived from motor triggered averages across hindbrain neurons. Behavioral subtype regressors were created by only considering bouts of a given type. The behavioral regressors that differentiate stimulus and rest periods were created by only considering motor events during those respective phases. Every cell with a correlation of at least 0.6 to at least one motor regressor was considered a “motor cell”. Since not all behaviors were observed in all imaging planes, cells were only assigned to a more specific category (such as swims versus all bouts) if the correlation to the more specific regressor was significantly higher (p < 0.01, bootstrap hypothesis test) than to any more general regressor. To identify heat responses across the whole brain spectral clustering was performed on a subset of cells since computing a pairwise correlation matrix of all cells and then performing clustering was not feasible. To arrive at a “canonical subset” for clustering the following filtering steps were performed. First for each cell the activity standard deviation during stimulus periods versus the standard deviation during rest periods was computed. All cells for which the standard deviation during stimulation was not greater than the standard deviation during rest were excluded, reasoning that those cells are unlikely to be stimulus modulated. Furthermore, all cells that had a correlation > 0.4 to a motor regressor were removed. This reduced the original number of cells from 699,840 down to 53,728 cells. For all these cells the pairwise correlations were computed. Any cells that did not have at least 20 partners across at least 3 fish that each explained 40 % of the cells variance were discarded. The reasoning behind this step was to enrich the considered set for canonical stimulus responses at the expense of discarding any rare response types. After this step 19,389 cells or 2.8 % of the original amount remained. We note however, that asking for merely one additional partner is responsible for 80 % of the observed reduction. After arriving at the filtered set of cells the trial-average response of each remaining cell was computed and spectral clustering with pair-wise correlations as similarity matrix asking for 6 total clusters was performed. Subsequently the cluster averages were used as regressors, probing the entire set of 699,840 cells. Each cell that had a correlation ≥ 0.6 was included in the final cell clusters. Across the resulting clusters about half the cells came from the original “canonical set” while the other half was made up of previously discarded cells. To identify heat modulated cells in specific brain regions we used the transformations computed during image registration of our experimental stacks together with our manual segmentation to assign each nuclear centroid to a brain region. From each brain region all cells that were correlated > 0.4 to any motor regressor were discarded. Since the number of cells in each region was much smaller than for the whole brain no further filtering was necessary and spectral clustering was performed on all the non-motor correlated cells. Overall the aim was to obtain as many different response types as possible. Setting the number of retrieved clusters to 6, at least one unstructured cluster was obtained in each region. The averages of remaining clusters were subsequently again used as regressors to identify responsive cells in the given region via a correlation ≥ 0.6 to the cluster average. Highly correlated clusters (r > 0.9) were subsequently merged. This simplified display of the data but had no influence on the further analysis (data not shown) as the highly correlated sister-clusters did not add explanatory power. To calculate ΔF/F0 values for reporting cell fluorescence we used the average across the first baseline period as the resting fluorescence F0.

Identification of heat and tap responsive cells

To identify heat responsive, tap responsive or mixed cells in the heat and tap experiments the timeseries of individual cells were categorized in the following manner. First using behavioral regressors all cells with a motor correlation > 0.4 were sorted out. Subsequently the period of sinewave heat stimulation (as this created the largest heat responses overall) as well as the period of the tap was used to calculate an activation score a averaged across the 25 trials T: a = 1 25 ∑ T | F pre ¯ − F stιm ¯ | σ pre Cells for which a > 2 were considered responsive for the given stimulus, cells where a < 0.75 were considered unresponsive and cells with intermediate scores were discarded. These criteria defined inclusion in the heat-responsive, tap-responsive or mixed categories. We note that these scores are fairly strict and overall identified fewer cells across Rh5/6 than the regressor based identification used in Figure 6 .

Nearest neighbor distance metrics

To compute nearest neighbor distance metrics for all cells the transformations obtained during image registration were used to map the centroids of all segmented nuclei into a common reference frame. For each cell of interest, the average distance to its two nearest neighbors of the same or comparison type was subsequently computed. Since nearest neighbor distances are influenced by the number of cells in each type as well, when comparing distances, the number of cells was down-sampled to the amount present in the smaller sample used for the comparison.

Information metrics

To compute the mutual information scores, the required entropies and joint entropies were computed using a jackknife estimate in order to reduce bias on small sample sizes ( Zahl, 1977 ). To compute an approximation of the mutual information between all region activity and the motor output principal component analysis was used for dimensionality reduction. This was necessary as memory requirements for computing the mutual information estimates grow exponentially with the number of traces considered. All mean cluster activities from all segmented regions were combined and principal components were extracted. The first five principal components explained > 95% of the total variance. These were subsequently used to compute the mutual information between these components and the motor output.

Determination of neurotransmitter types

To determine the neurotransmitter type of heat responsive cells in Rhombomeres 5/6 nuclear GCaMP6s was imaged using the heat stimulus in the presence of either a red vGlut2a (excitatory neurons) or red Gad1b (inhibitory neurons) label. Cells were subsequently segmented as described above and cell type cluster averages were used to identify cells of each type present in Rhombomere 5/6. Cells were subsequently cross-referenced with the transmitter label and expression ratios (positive / negative) were determined. Trigeminal ablations For trigeminal ablations the hindbrain of fish expressing cytoplasmic GCaMP6s was first imaged using the heat stimulus to identify heat responsive neurons. Subsequently individual cells in one trigeminal ganglion were ablated using 25 ms long pulses of high laser power (300 mW at sample at 850 nm), starting ventrally in the trigeminal and moving dorsally. In total around 40 pulses were delivered in each experiment resulting in partial destruction of the trigeminal ganglion. Damage was subsequently assessed by anatomical imaging. 2 fish were discarded right after the anatomical assessment, because no damage to the trigeminal was visible. A further set of 4 fish died before completion of the experiment and were therefore not analyzed. All 5 fish that survived the functional post-experiment were still healthy the following day and were included in the analysis. No fish were discarded after analyzing the functional data.

Circuit model

The circuit model was constructed as a feed-forward model from stimulus to behavioral output. The activation A in each stage (Trigeminal ganglion - Rh5/6 region of the hindbrain - Motor cells - Behavior) was described as a linear combination of the activations in the previous stage, optionally convolved with a linear filter f (Trigeminal and Rh5/6) and passed through an output nonlinearity g (Trigeminal and Rh5/6): A i = g ( ( A i − 1 β T ) ∗ f ) The linear filter f was necessary to explain the transformations between different response dynamics. The linear parts of each model stage ( β and f ) were fit via Markov-chain Monte Carlo using PyMC3 ( Salvatier et al., 2016 ). Since all our likelihoods were differentiable allowing gradient calculation we used the no-uturn sampler in PyMC3 for sample generation ( Hoffman and Gelman, 2014 ). Cubic output nonlinearities were subsequently fit using least squares optimization in Python. The time-resolution Δt of all models was set to 0.2 s the same as the time resolution of our interpolated imaging traces. We note that our model does not operate on raw neuronal or behavior traces, but the inputs and outputs of all stages were z-scored. This is a more conservative approach due to the nonlinear relationship between calcium signals and neuronal activity. Stimulus encoding in the trigeminal ganglion To capture the encoding of the heat stimulus in the observed trigeminal calcium activity the filter was parametrized as an exponential on/off filter akin to a calcium kernel. f ( Δ t ; τ off , τ on ) = e − Δ t τ off − ( 1 − e − Δ t τ on ) | − 20 s < Δ t ≤ 0 To simplify computation, convolution in the model was expressed via piecewise multiplications such that the predicted output activity A(t) of the trigeminal ON and OFF cells at each timepoint t was described in terms of the temperature stimulus T(t) as A ( t ; β , f ) = β ∑ Δ t T ( t − Δ t ) f ( Δ t ) ∣ − 20 s < Δ t ≤ 0 The prior parameter distributions used for the MCMC process were as follows: β ~ N ( 0 , 2 ) τ on ~ N ( 0 , 5 ) τ off ~ N ( 0 , 200 ) σ ~ | N ( 0 , 1 ) | The likelihood of the model was subsequently defined in terms of a normal distribution centered according to: ℒ ( A ( β , f ) ) ∣ A ~ N ( A ( t ; β , f ) , σ ) Transformation of trigeminal activity in Rh5/6 To express the activity in this hindbrain region in terms of the activity of trigeminal cell-types we used a filter-parametrization that would allow for integrating as well as differentiating filters: f ( Δ t ; s , τ 1 , τ 2 ) = s Δ t e − Δ t τ 0 + ( 1 − Δ t ) e − Δ t τ 1 ∣ − 4 s ≤ Δ t ≤ 0 To simplify computation, convolution in the model was again expressed via piecewise multiplications such that the predicted output activity A(t) of a cell type in Rh5/6 was described in terms of the input activities ON(t) and OFF(t) as: A ( t ; β ON , β OFF , f ) = β ON ∑ Δ t ON ( t − Δ t ) f ( Δ t ) + β OFF ∑ Δ t OFF ( t − Δ t ) f ( Δ t ) ∣ − 4 s ≤ Δ t ≤ 0 Since trigeminal neurons are largely glutamatergic, we enforced all trigeminal inputs to have activating effects in our models. However, to express the Fast-ON, Fast-OFF and Delayed-OFF types in terms of their inputs required inhibition. Those Rh 5/6 types were therefore expressed in terms of a strictly activating trigeminal input and a potentially inhibitory input from either the Slow-ON or Slow-OFF type. The prior distributions of β reflected this constraint: Slow ­ ON , Slow ­ OFF β ON ~ | N ( 0 , 2 ) | β OFF ~ | N ( 0 , 2 ) | Fast ­ ON , Fast ­ OFF β ON ~ | N ( 0 , 2 ) | β OFF ~ N ( 0 , 2 ) Delayed OFF β ON ~ N ( 0 , 2 ) β OFF ~ | N ( 0 , 2 ) | The priors of the remaining parameters were shared between all models: τ 1 ~ | N ( 0 , 10 ) | τ 2 ~ | N ( 0 , 10 ) | s ~ | N ( 0 , 5 ) | σ ~ | N ( 0 , 1 ) | The likelihood of the model was subsequently defined in terms of a normal distribution centered according to: ℒ ( A ( β ON , β OFF , f ) ) ∣ A ~ N ( A ( t ; β ON , β OFF , f ) , σ ) Rate coding steps The last two stages of the model were implemented as simple linear regression steps relating the weighted sum of input activity to output activity. A out ( β ) = A in β T The prior distributions on the parameters β and the error standard deviation s were defined as: β i ~ N ( 0 , 2 ) σ ~ | N ( 0 , 1 ) | And the model likelihood was expressed according to: ℒ ( A ( β ) ) ∣ A ~ N ( A out ( β ) , σ ) Simulation of filter derivation To test the influence of noise in our data on filtering, a simple simulation was performed. The temperature stimulus T(t) was used as the input to a sharp high pass filter ( Figure S5 I , orange line) to create the discrete difference trace dT(t) of the temperature stimulus: δT ( t ) = T ( t ) − T ( t − 1 ) Both the original temperature trace and the difference trace were subsequently corrupted by i.i.d. gaussian noise: I ( t ) = T ( t ) + ε ( t ) ∣ ε ( t ) ~ N ( 0 , 0.2 ) O ( t ) = δT ( t ) + ε ( t ) ∣ ε ( t ) ~ N ( 0 , 0.2 ) These traces were subsequently used as inputs and outputs to derive a model using the same filter parametrization and strategy as above for cell types in Rhombomeres 5 and 6. As can be seen in Figure S5I , the derived filter (black trace) extends over considerably longer timescales than the true differencing filter (orange trace).

QUANTIFICATION AND STATISTICAL ANALYSIS

All data analysis was performed in python. All confidence intervals and standard errors displayed in figures are derived using bootstrapping across cells and appropriate n-values are indicated in the figure legend. As indicated in the figure legends bootstrap hypothesis testing was used to determine significance testing except in the case of trigeminal ablations where a nonparametric ranksum test was used across 3 fish to assess significance. Counts either refer to cells or individual fish as indicated in the figure legends where applicable. Data shuffles Since the stimulus presented on each imaging plane is repeated three times, cell activity that reflects the stimulus should follow the same repeat structure. The data shuffles used for the analysis of heat modulated cells in Supplementary Figures 1 D-E and 4 E-F were therefore generated as follows. Each cell’s activity trace was split into the three repeats and each repeat was independently circularly permuted. This should break the repeat structure but leave in-repeat timescale structures intact. The activity generated this way was subsequently subjected to the same whole-brain or regional clustering approaches detailed above. Since regression-based analysis was used to identify motor-correlated cells a different approach was used to create the controls in supplemental Figure 3 A-C . In this case the whole activity trace of a cell was circularly permuted with respect to the motor regressors again keeping the overall structure of the activity data intact. The shuffled activity was subsequently probed with the motor regressors as described above.

DATA AND SOFTWARE AVAILABILITY

Raw experimental data (~40 GB) and custom written analysis as well as acquisition software will be made available by the authors upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

All experiments were conducted on 6-7 days post fertilization zebrafish homozygous for the nacre mutation expressing the transgenes indicated below. The sex of the larva is not defined at this early stage. Fish were fed paramecia from day 5 onwards. All experiments followed the guidelines of the National Institutes of Health and were approved by the Standing Committee on the Use of Animals in Research of Harvard University. All experiments used for mapping heat responsive neurons and for model derivation used fish expressing nuclear Elavl3:H2B-GCaMP6s fish ( Freeman et al., 2014 ). Experiments combining heat and taps as well as trigeminal ablations were performed in fish expressing cytoplasmic Elavl3:Gcamp6s ( Kim et al., 2017 ). We note that while the decay time-constants for cytoplasmic and nuclear GCaMP are different, we found this effect to be negligible given the slow timescales of our temperature stimulus. To develop the scan stabilization protocol Elavl3:H2B-RFP fish expressing RFP in all neuronal nuclei were used ( Randlett et al., 2015 ). To identify neurotransmitter types of heat modulated neurons, progeny of crosses between Elavl3:H2B-Gcamp6s ( Freeman et al., 2014 ) and either vglut2a:mCherry ( Satou et al., 2013 ) or gad1B:dsRed ( Satou et al., 2013 ) expressing fish were used.

METHOD DETAILS Since almost all analysis in this study was performed in an automated manner, no blinding or randomization was performed. All sample sized were fixed before the start of experiments and all animals were analyzed except those that freed themselves from the preparation during imaging or that died during the experimental procedure.

Imaging and Behavior

Larval zebrafish were embedded in 2.5% medium melt agarose (Fisher scientific, USA) and their tails were freed the night before the experiment. Experiments were conducted in a custom built 2-photon microscope and run using custom written software in C# (Microsoft, USA). Heat stimuli were delivered using a 1W 980nm fiber-coupled diode laser (Roithner, Austria) coupled into a collimator (Aistana Inc., USA) placed under the microscope objective 4 mm in front of and 1.2 mm above the head of the zebrafish larva pointing downwards at an angle of 16.5 degrees. The laser power was controlled by the computer via a laser diode driver (Thorlabs, USA). We note that the 980 nm laser itself did not excite GCaMP fluorescence due to the low photon density. The main mapping experiments consisted of the imaging of 30 individual planes, spaced 2.5 μm apart. In each plane 3 trials of the stimulus depicted in Figure 1B were presented. The heat and tap experiments consisted of imaging 4 individual planes, spaced 5 μm apart. In each plane 25 trials of the stimulus depicted in Figure S1G were presented. To avoid excessive heating of the preparation by scanning over the eyes, custom exclusion masks were created for each experiment in which the eyes were in the field of view, restricting the scan-lines such that the eyes were excluded from the field of view. Imaging was performed at 2.4 Hz - 3 Hz and all imaging data was interpolated to a 5 Hz timebase before further analysis. Tail-tracking data was acquired at 100 Hz. All behavioral features such as identifying and classifying bouts were performed at this timebase, however, for all comparisons of behavior and imaging, the behavioral data was downsampled to a 5 Hz timebase.

Stimulus characterization

To measure the temperature induced in larval zebrafish by the laser stimulus a thermistor with the same dimension and absorption characteristics as larval zebrafish was used (Warner instruments, USA; ( Haesemeyer et al., 2015 )). The thermistor was embedded in the same low-melting point agarose as the larval zebrafish and the laser was positioned accordingly ( Madelaine et al., 2017 ). For each stimulus used, temperature changes in the thermistor were subsequently recorded. While tissue absorption is not homogeneous for larval zebrafish, modeling of thermal flux revealed that over the short length scales at this age there won’t be an uneven distribution of temperature within the head of larval zebrafish due to the thermal conductivity of water within about 250 ms of stimulus onset. In other words, equilibration happens much faster than the timescales of our stimulus (not shown). Furthermore, the size of the laser spot was set such that it is larger than the head of the fish. Image stabilization To counter heat-induced deformations of the preparation, before each plane was scanned a +/- 5 μm sized pre-stack consisting of 21 slices spaced 0.4 μm apart was acquired. During scanning, each acquired plane was cross-correlated with each plane in the pre-stack and using a low-pass filter the position of the objective was adjusted online so as to minimize z-drift. Since the heat induced drift observed in Elavl3:H2B-RFP stacks followed very reproducible kinetics these were used to predict the movement induced by heating during our experiments to induce very slight movements in the predicted direction. This served to overcome the delay induced by the low-pass filter. A second alignment step was performed post-acquisition. Here each individual imaging plane was assigned to the most likely plane in the pre-stack that it corresponds to via correlation. The fluorescence timeseries of each segmented region was subsequently corrected by normalizing it with the resting fluorescence of this region in the pre-stack.

Image segmentation

Cells in each plane were segmented anatomically. To correct for motion artifacts individual planes in each timeseries were re-aligned based on image cross-correlations ( Miri et al., 2011 ). For nuclear GCaMP experiments, individual nuclei were subsequently segmented using Cell Profiler ( Carpenter et al., 2006 ). To resolve merged nuclei in areas of low contrast, objects larger than a typical nuclear size were divided into subregions based on pixel timeseries correlations. For experiments using cytoplasmic GCaMP, custom written software was used to identify individual cells. Each image was processed using a minimum and a maximum filter tuned to nuclear size. Points where the differences between the two filtered results crossed a threshold (at least 20 % of the average region brightness) were considered potential cell centroids. After including seed pixels around the centroids individual regions were grown in a greedy manner, incorporating pixels that were better correlated to their own current average activity than to neighboring averages. Resulting correlations masks were intersected with anatomical masks based on cell size to obtain the final segmentation. Registration and annotation 3D image registration based on CMTK ( Rohlfing and Maurer, 2003 ) was used to create a nuclear GCaMP-6s reference stack to which all experimental stacks were registered as described elsewhere ( Portugues et al., 2014 ; Randlett et al., 2015 ). The reference brain was used to annotate anatomical regions of interest based on Z-Brain annotations ( Randlett et al., 2015 ).

Swim bout identification and classification

Swim bouts were identified as described previously ( Portugues et al., 2015 ), namely based on the windowed standard deviation of the tail cumulative angle trace crossing a threshold. The mode-centered cumulative angle trace a of each bout was subsequently used to assign the following score: bias = ∑ a a ∑ a | a | Right flicks were defined as bouts with a bias < -0.8, left flicks as bouts with a bias > 0.8 and all other bouts were defined as “swims”. These cutoffs were chosen since the histogram of all bias scores has a minimum at those points. Clustering of heat and motor related activity To identify motor related activity, motor regressors were created by convolving a bout start trace with an exponential calcium kernel with a decay half time of 3s. This decay time was derived from motor triggered averages across hindbrain neurons. Behavioral subtype regressors were created by only considering bouts of a given type. The behavioral regressors that differentiate stimulus and rest periods were created by only considering motor events during those respective phases. Every cell with a correlation of at least 0.6 to at least one motor regressor was considered a “motor cell”. Since not all behaviors were observed in all imaging planes, cells were only assigned to a more specific category (such as swims versus all bouts) if the correlation to the more specific regressor was significantly higher (p < 0.01, bootstrap hypothesis test) than to any more general regressor. To identify heat responses across the whole brain spectral clustering was performed on a subset of cells since computing a pairwise correlation matrix of all cells and then performing clustering was not feasible. To arrive at a “canonical subset” for clustering the following filtering steps were performed. First for each cell the activity standard deviation during stimulus periods versus the standard deviation during rest periods was computed. All cells for which the standard deviation during stimulation was not greater than the standard deviation during rest were excluded, reasoning that those cells are unlikely to be stimulus modulated. Furthermore, all cells that had a correlation > 0.4 to a motor regressor were removed. This reduced the original number of cells from 699,840 down to 53,728 cells. For all these cells the pairwise correlations were computed. Any cells that did not have at least 20 partners across at least 3 fish that each explained 40 % of the cells variance were discarded. The reasoning behind this step was to enrich the considered set for canonical stimulus responses at the expense of discarding any rare response types. After this step 19,389 cells or 2.8 % of the original amount remained. We note however, that asking for merely one additional partner is responsible for 80 % of the observed reduction. After arriving at the filtered set of cells the trial-average response of each remaining cell was computed and spectral clustering with pair-wise correlations as similarity matrix asking for 6 total clusters was performed. Subsequently the cluster averages were used as regressors, probing the entire set of 699,840 cells. Each cell that had a correlation ≥ 0.6 was included in the final cell clusters. Across the resulting clusters about half the cells came from the original “canonical set” while the other half was made up of previously discarded cells. To identify heat modulated cells in specific brain regions we used the transformations computed during image registration of our experimental stacks together with our manual segmentation to assign each nuclear centroid to a brain region. From each brain region all cells that were correlated > 0.4 to any motor regressor were discarded. Since the number of cells in each region was much smaller than for the whole brain no further filtering was necessary and spectral clustering was performed on all the non-motor correlated cells. Overall the aim was to obtain as many different response types as possible. Setting the number of retrieved clusters to 6, at least one unstructured cluster was obtained in each region. The averages of remaining clusters were subsequently again used as regressors to identify responsive cells in the given region via a correlation ≥ 0.6 to the cluster average. Highly correlated clusters (r > 0.9) were subsequently merged. This simplified display of the data but had no influence on the further analysis (data not shown) as the highly correlated sister-clusters did not add explanatory power. To calculate ΔF/F0 values for reporting cell fluorescence we used the average across the first baseline period as the resting fluorescence F0.

Identification of heat and tap responsive cells

To identify heat responsive, tap responsive or mixed cells in the heat and tap experiments the timeseries of individual cells were categorized in the following manner. First using behavioral regressors all cells with a motor correlation > 0.4 were sorted out. Subsequently the period of sinewave heat stimulation (as this created the largest heat responses overall) as well as the period of the tap was used to calculate an activation score a averaged across the 25 trials T: a = 1 25 ∑ T | F pre ¯ − F stιm ¯ | σ pre Cells for which a > 2 were considered responsive for the given stimulus, cells where a < 0.75 were considered unresponsive and cells with intermediate scores were discarded. These criteria defined inclusion in the heat-responsive, tap-responsive or mixed categories. We note that these scores are fairly strict and overall identified fewer cells across Rh5/6 than the regressor based identification used in Figure 6 .

Nearest neighbor distance metrics

To compute nearest neighbor distance metrics for all cells the transformations obtained during image registration were used to map the centroids of all segmented nuclei into a common reference frame. For each cell of interest, the average distance to its two nearest neighbors of the same or comparison type was subsequently computed. Since nearest neighbor distances are influenced by the number of cells in each type as well, when comparing distances, the number of cells was down-sampled to the amount present in the smaller sample used for the comparison.

Information metrics

To compute the mutual information scores, the required entropies and joint entropies were computed using a jackknife estimate in order to reduce bias on small sample sizes ( Zahl, 1977 ). To compute an approximation of the mutual information between all region activity and the motor output principal component analysis was used for dimensionality reduction. This was necessary as memory requirements for computing the mutual information estimates grow exponentially with the number of traces considered. All mean cluster activities from all segmented regions were combined and principal components were extracted. The first five principal components explained > 95% of the total variance. These were subsequently used to compute the mutual information between these components and the motor output.

Determination of neurotransmitter types

To determine the neurotransmitter type of heat responsive cells in Rhombomeres 5/6 nuclear GCaMP6s was imaged using the heat stimulus in the presence of either a red vGlut2a (excitatory neurons) or red Gad1b (inhibitory neurons) label. Cells were subsequently segmented as described above and cell type cluster averages were used to identify cells of each type present in Rhombomere 5/6. Cells were subsequently cross-referenced with the transmitter label and expression ratios (positive / negative) were determined. Trigeminal ablations For trigeminal ablations the hindbrain of fish expressing cytoplasmic GCaMP6s was first imaged using the heat stimulus to identify heat responsive neurons. Subsequently individual cells in one trigeminal ganglion were ablated using 25 ms long pulses of high laser power (300 mW at sample at 850 nm), starting ventrally in the trigeminal and moving dorsally. In total around 40 pulses were delivered in each experiment resulting in partial destruction of the trigeminal ganglion. Damage was subsequently assessed by anatomical imaging. 2 fish were discarded right after the anatomical assessment, because no damage to the trigeminal was visible. A further set of 4 fish died before completion of the experiment and were therefore not analyzed. All 5 fish that survived the functional post-experiment were still healthy the following day and were included in the analysis. No fish were discarded after analyzing the functional data.

Circuit model

The circuit model was constructed as a feed-forward model from stimulus to behavioral output. The activation A in each stage (Trigeminal ganglion - Rh5/6 region of the hindbrain - Motor cells - Behavior) was described as a linear combination of the activations in the previous stage, optionally convolved with a linear filter f (Trigeminal and Rh5/6) and passed through an output nonlinearity g (Trigeminal and Rh5/6): A i = g ( ( A i − 1 β T ) ∗ f ) The linear filter f was necessary to explain the transformations between different response dynamics. The linear parts of each model stage ( β and f ) were fit via Markov-chain Monte Carlo using PyMC3 ( Salvatier et al., 2016 ). Since all our likelihoods were differentiable allowing gradient calculation we used the no-uturn sampler in PyMC3 for sample generation ( Hoffman and Gelman, 2014 ). Cubic output nonlinearities were subsequently fit using least squares optimization in Python. The time-resolution Δt of all models was set to 0.2 s the same as the time resolution of our interpolated imaging traces. We note that our model does not operate on raw neuronal or behavior traces, but the inputs and outputs of all stages were z-scored. This is a more conservative approach due to the nonlinear relationship between calcium signals and neuronal activity. Stimulus encoding in the trigeminal ganglion To capture the encoding of the heat stimulus in the observed trigeminal calcium activity the filter was parametrized as an exponential on/off filter akin to a calcium kernel. f ( Δ t ; τ off , τ on ) = e − Δ t τ off − ( 1 − e − Δ t τ on ) | − 20 s < Δ t ≤ 0 To simplify computation, convolution in the model was expressed via piecewise multiplications such that the predicted output activity A(t) of the trigeminal ON and OFF cells at each timepoint t was described in terms of the temperature stimulus T(t) as A ( t ; β , f ) = β ∑ Δ t T ( t − Δ t ) f ( Δ t ) ∣ − 20 s < Δ t ≤ 0 The prior parameter distributions used for the MCMC process were as follows: β ~ N ( 0 , 2 ) τ on ~ N ( 0 , 5 ) τ off ~ N ( 0 , 200 ) σ ~ | N ( 0 , 1 ) | The likelihood of the model was subsequently defined in terms of a normal distribution centered according to: ℒ ( A ( β , f ) ) ∣ A ~ N ( A ( t ; β , f ) , σ ) Transformation of trigeminal activity in Rh5/6 To express the activity in this hindbrain region in terms of the activity of trigeminal cell-types we used a filter-parametrization that would allow for integrating as well as differentiating filters: f ( Δ t ; s , τ 1 , τ 2 ) = s Δ t e − Δ t τ 0 + ( 1 − Δ t ) e − Δ t τ 1 ∣ − 4 s ≤ Δ t ≤ 0 To simplify computation, convolution in the model was again expressed via piecewise multiplications such that the predicted output activity A(t) of a cell type in Rh5/6 was described in terms of the input activities ON(t) and OFF(t) as: A ( t ; β ON , β OFF , f ) = β ON ∑ Δ t ON ( t − Δ t ) f ( Δ t ) + β OFF ∑ Δ t OFF ( t − Δ t ) f ( Δ t ) ∣ − 4 s ≤ Δ t ≤ 0 Since trigeminal neurons are largely glutamatergic, we enforced all trigeminal inputs to have activating effects in our models. However, to express the Fast-ON, Fast-OFF and Delayed-OFF types in terms of their inputs required inhibition. Those Rh 5/6 types were therefore expressed in terms of a strictly activating trigeminal input and a potentially inhibitory input from either the Slow-ON or Slow-OFF type. The prior distributions of β reflected this constraint: Slow ­ ON , Slow ­ OFF β ON ~ | N ( 0 , 2 ) | β OFF ~ | N ( 0 , 2 ) | Fast ­ ON , Fast ­ OFF β ON ~ | N ( 0 , 2 ) | β OFF ~ N ( 0 , 2 ) Delayed OFF β ON ~ N ( 0 , 2 ) β OFF ~ | N ( 0 , 2 ) | The priors of the remaining parameters were shared between all models: τ 1 ~ | N ( 0 , 10 ) | τ 2 ~ | N ( 0 , 10 ) | s ~ | N ( 0 , 5 ) | σ ~ | N ( 0 , 1 ) | The likelihood of the model was subsequently defined in terms of a normal distribution centered according to: ℒ ( A ( β ON , β OFF , f ) ) ∣ A ~ N ( A ( t ; β ON , β OFF , f ) , σ ) Rate coding steps The last two stages of the model were implemented as simple linear regression steps relating the weighted sum of input activity to output activity. A out ( β ) = A in β T The prior distributions on the parameters β and the error standard deviation s were defined as: β i ~ N ( 0 , 2 ) σ ~ | N ( 0 , 1 ) | And the model likelihood was expressed according to: ℒ ( A ( β ) ) ∣ A ~ N ( A out ( β ) , σ ) Simulation of filter derivation To test the influence of noise in our data on filtering, a simple simulation was performed. The temperature stimulus T(t) was used as the input to a sharp high pass filter ( Figure S5 I , orange line) to create the discrete difference trace dT(t) of the temperature stimulus: δT ( t ) = T ( t ) − T ( t − 1 ) Both the original temperature trace and the difference trace were subsequently corrupted by i.i.d. gaussian noise: I ( t ) = T ( t ) + ε ( t ) ∣ ε ( t ) ~ N ( 0 , 0.2 ) O ( t ) = δT ( t ) + ε ( t ) ∣ ε ( t ) ~ N ( 0 , 0.2 ) These traces were subsequently used as inputs and outputs to derive a model using the same filter parametrization and strategy as above for cell types in Rhombomeres 5 and 6. As can be seen in Figure S5I , the derived filter (black trace) extends over considerably longer timescales than the true differencing filter (orange trace).

Supplementary Material 1