Ongoing activity both arises from and contributes to fluctuations in neuronal excitability. To gain insight into the dominant network interactions that mediate these processes, we analysed ongoing activity in the optic tectum (OT). To measure this activity, we performed light-sheet calcium imaging of transgenic larvae expressing GCaMP6s (elavl3:H2B-GCaMP6s, 6 dpf, n=14), initially under constant environmental conditions (i.e. without presenting visual stimuli). We imaged a 375x410 × 75 μm volume encompassing OT and adjacent brain structures, at a rate of 5 volumes per second (Figure 1A–B). Following motion correction, each imaging plane was automatically segmented to individual neurons (Kawashima et al., 2016), and spikes were inferred by deconvolution of the extracted fluorescence signal (Friedrich et al., 2017). Imaged volumes were registered to a reference brain (Avants et al., 2009; Marquart et al., 2015) and somata in the stratum periventriculare (SPV), which comprise the majority of tectal neurons, were considered for further analysis (Figure 1B, 1.4‧10⁴±3‧10³ tectal SPV cells per fish, mean ± SD).

Figure 1 with 1 supplement see all

Download asset Open asset

Figure 1—source data 1

Data provided as a MATLAB structure.

https://cdn.elifesciences.org/articles/78381/elife-78381-fig1-data1-v1.zip

Download elife-78381-fig1-data1-v1.zip

The spatial pattern of activity correlations was compatible with recurrent excitatory interactions between tectal cells. To show this, we analyzed the Pearson correlation between the ongoing activity of each neuron (‘seed’), and the activity of all other neurons during the entire experimental session. Correlation coefficients smoothly and monotonically decreased as a function of the Euclidean distance between neurons (Figure 1C and Figure 1—figure supplement 1A–C). We confirmed this finding using two-photon imaging, suggesting it is not a consequence of light scattering in the light-sheet microscope but rather a feature of tectal activity (Figure 1—figure supplement 1C). We projected neuronal coordinates into two dimensions (anterior-posterior and medial-lateral tectal axes, see Materials and methods), producing a 2D correlation map for each seed neuron (Figure 1D). A 2D Gaussian was fitted to each map, excluding the seed neuron itself, yielding good fits for most cells (r²=0.42 ± 0.17, mean ± SD, Figure 1D and Figure 1—figure supplement 1A–C). The peaks of the fitted Gaussians showed close correspondence to the locations of seed neurons (Figure 1E and Figure 1—figure supplement 1A–C). Comparably good fits were obtained by fitting an exponential decay function (r²=0.39 ± 0.18, mean ± SD), and again the peaks of fitted exponentials corresponded to the locations of seed neurons (Figure 1—figure supplement 1A–C). These data are compatible with a recurrent network in which excitatory interactions smoothly decay with distance from each cell.

A salient aspect of tectal activity is the occurrence of ‘spontaneous’ bursts, during which spatially compact groups of neurons briefly fire in synchrony (Bianco and Engert, 2015; Romano et al., 2015; Avitan et al., 2017). To detect bursts, we evaluated the spatial distribution of cells active during one-second windows centred on population activity peaks (Figure 1F). For each window, we identified the cells participating in a burst using density-based DBSCAN clustering (Ester et al., 1996; Figure 1—figure supplement 1, Materials and methods). The temporal extent of each burst was then defined as the period of non-zero Poisson-filtered activity of participating neurons (Figure 1F, iv). Overall, 54,424 bursts were identified in 14 fish. On average, each activity peak was associated with 1.5±2.0 localised bursts (mean ± SD), bursts occurred at a rate of 46±11 per minute and appeared uniformly distributed in space and time (Figure 1G). Bursts were brisk, phasic events, with the majority of spiking occurring over a timecourse of around one second (Figure 1H).

The statistics of bursting suggested tectal activity is shaped by a second interaction, where activity results in prolonged reduction in the excitability of tectal cells. As has been previously reported, the sizes and durations of bursts were highly variable (95±183 neurons per burst, 2.5±3 s duration, mean ± SD), and distributed according to a power law (Figure 1I–J), suggestive of a critical phenomenon with ‘avalanche’ dynamics (Ponce-Alvarez et al., 2018). A universal feature in such phenomena is a self-limiting element responsible for both event termination and suppression of subsequent events. Therefore, we analysed how participating in a burst influences subsequent spiking activity (Figure 1K) and burst participation probability (Figure 1L) for individual neurons. Cells were less active and less likely to take part in another burst for tens of seconds. These observations are compatible with activity-dependent suppression that has a self-limiting effect on bursting as well as a lasting effect extending beyond burst termination.

In sum, ongoing tectal activity is compatible with two types of intercellular interaction that might shape network dynamics. First, activity correlation is suggestive of excitatory recurrent connectivity that smoothly decays with distance. Second, the dynamics of bursting suggest prolonged activity-dependent suppression of excitability.

Next, we used a computational modeling framework with defined interactions between spiking activity and excitability state to explore the possibility that tectal bursting emerges from the recurrent interaction motifs suggested by our analysis of ongoing activity. Theoretical network analysis has proposed that recurrent excitatory connectivity, balanced by long-range inhibition, produces stable, self-sustaining and spatially defined ‘blobs’ of activity (Amari, 1977). We reasoned that by making the inhibition in such a model depend on a slow time-integration of activity, transient bursts would emerge instead of stable blobs and would be followed by prolonged activity suppression.

To test this idea, we used a probabilistic spiking network simulation based on a linear-nonlinear-Poisson (LNP) formalism. This framework has been successfully used to relate activity in neuronal populations to past activity and sensory input (Pillow et al., 2005; Pillow et al., 2008; Truccolo et al., 2005), and its non-deterministic aspect fits the apparent stochastic nature of tectal bursting. In LNP models, the instantaneous firing rate of each neuron (λ) is determined by an exponentiated ‘linear drive’, computed as the sum of filtered inputs from all neurons in the network plus optional external inputs; spikes are then emitted by simulating a stochastic Poisson process according to the firing rate (Figure 2A, Pillow et al., 2008). In our model, the linear drive can be likened to the excitability state of model neurons. Our aim was to test if a network defined by a recurrent interaction motif, applied equally across all neurons, might be sufficient to reproduce the bursting activity observed in OT, rather than fit the connectivity of each neuron. Our model is thus defined by a set of only seven parameters, describing the spatial and temporal aspects of ‘excitatory’ (‘E’) and ‘inhibitory’ (‘I’) interactions in the network, resulting in increase and decrease in excitability, respectively. The connection weights between all pairs of neurons decay with Euclidian distance according to Gaussian functions with variances ${\sigma }^{\left(E\right)},{\sigma }^{\left(I\right)}$, scaled by fixed gains $g^{\left(E\right)},g^{\left(I\right)}$. Both interactions decay exponentially over time according to time constants ${\tau }^{\left(E\right)},{\tau }^{\left(I\right)}$. A fixed bias term, µ, is summed with all other inputs to determine the total linear drive (see Materials and methods). The locations of neurons in the simulation were taken from a single fish registered to the reference brain. Notably, given that the same connectivity rule is applied to every tectal cell, the network is fully described by the set of seven parameters and the locations of the neurons.

Figure 2

Download asset Open asset

Figure 2—source data 1

Data provided as a MATLAB structure.

https://cdn.elifesciences.org/articles/78381/elife-78381-fig2-data1-v1.zip

Download elife-78381-fig2-data1-v1.zip

In accordance with our hypothesis, using this generative model to simulate network activity resulted in spatially localised bursting (Video 1). Bursting was a robust emergent property that was observed across a broad range of parameter values (see below). Bursts could be identified by analysing the simulation spiking output using the same method as for experimental data (described above) and were uniformly distributed in space and time with a range of sizes and durations (Figure 2B–C).

Video 1

Download asset

In previous studies, localised bursting in the zebrafish OT has been attributed to sub-networks with enhanced connectivity, termed neuronal ‘assemblies’ (Romano et al., 2015; Avitan et al., 2017; Marachlian et al., 2018; Mölter et al., 2018; Diana et al., 2019). A PCA-promax algorithm was previously used to identify such assemblies based on the tendency of neurons to be co-active (Romano et al., 2015). We applied this method to our simulated tectal network and despite the fact that, by design, there were no subgroups of cells with preferential connectivity, PCA-promax yielded multiple spatially localised ‘assemblies’ (Figure 2D), as was also the case for recorded data (Figure 2E). This analysis demonstrates that repeated, synchronised bursting of tectal neurons can occur in a network where every cell has the same distance-based recurrent connectivity and therefore that this feature of activity does not, by itself, imply the existence of preferentially connected subnetworks.

Our simulated tectal network enabled us to examine the linear drive of model neurons — in effect their excitability state. This revealed (Figure 2F, Video 2) that bursts tend to occur in regions with higher than average excitability (bright green) and leave behind a region of low excitability (black), supressing activity in the same region for a long period of time.

Video 2

Download asset

In sum, our network model reveals that a simple recurrent interaction motif is sufficient to explain tectal bursting.

We next evaluated the capacity of the model to produce bursting with biologically relevant statistics. We found that bursting occurred robustly across a range of parameter values, albeit with burst statistics – such as average burst size and duration – varying across parameter space.

Systematically varying the properties of the inhibitory interaction (space constant ${\sigma }^{\left(I\right)}$, time constant ${\tau }^{\left(I\right)}$ and gain $g^{\left(I\right)}$) produced a wide range of mean burst frequencies (0.1–245 bursts per minute) and sizes (18–356 neurons per burst), while some combinations of parameters completely abolished bursting, at least for the 30 min duration tested (Figure 3A). For example, eliminating the inhibitory interaction ($g^{\left(I\right)}=0$) abolished bursting (Figure 3A, green dashed line) since all neurons remain tonically active without a suppressive interaction to balance recurrent excitation. Very large values of ${\sigma }^{\left(I\right)}$, approach the condition where inhibition is independent of distance: This resulted in very infrequent but large bursts (Figure 3A, red dashed line), as the entire network simultaneously became non-excitable and then synchronously recovered. Increasing the temporal decay constant (${\tau }^{\left(I\right)}$) produced rare, small bursts (Figure 3A, orange dashed line) by prolonging the recovery time. We also confirmed the capacity of the model to produce bursting with an exponential rather than Gaussian relationship between distance and connection strength (Figure 3—figure supplement 1A), and the robustness of the bursting phenomenon to stochastic non-uniformities in connections by pruning random subsets of excitatory and inhibitory connections (Figure 3—figure supplement 1B–C). These analyses indicate that bursting is a robust phenomenon but that bursting statistics are modulated by the spatiotemporal parameters defining recurrent interactions.

Figure 3 with 1 supplement see all

Download asset Open asset

Figure 3—source data 1

Data provided as a MATLAB structure.

https://cdn.elifesciences.org/articles/78381/elife-78381-fig3-data1-v1.zip

Download elife-78381-fig3-data1-v1.zip

Next, we sought to optimise the model parameters such that the stochastic, emergent property of bursting would best match tectal bursting statistics. We started with the set of parameters used in the previous section (‘baseline model’), which generated bursting but deviated from biologically observed statistics with regard to frequency, size and duration of burst events (Figure 2C-D). We then used evolutionary multi-objective optimization (EMOO, Kalyanmoy, 2001, Figure 3—figure supplement 1D and Materials and methods) to tune the parameters according to three target objectives, namely the distribution of burst sizes (Figure 1I) and durations (Figure 1J), and the timecourse of burst-triggered activity suppression (Figure 1K). With each successive generation, EMOO produces a population of models that better approximates the Pareto optimal set of solutions (Pareto front), that is the set of models where no individual model performs better than any other in all three objectives (Figure 3—figure supplement 1E). Following thirty EMOO generations, we chose the solution at the elbow of the estimated Pareto front (the model closest to the origin in the z-scored loss space, Figure 3—figure supplement 1E, green), and used it as a starting point for a local search to fine-tune the parameters.

The resulting optimized model was successful in reproducing the bursting statistics of the zebrafish tectum (Figure 3—figure supplement 1E, blue). This was true for both the features used as optimization objectives (Figure 3B–D), as well as two that were not used for parameter tuning, namely burst-triggered burst participation probability (Figure 3E) and the burst temporal profile (Figure 3F).

The parameter values of the optimized model support the idea that bursting in OT arises from distinct spatiotemporal patterns of excitatory versus suppressive network interactions. Excitatory interactions are strong, short range (${\sigma }^{\left(E\right)}=4.5\mu m$) and decay rapidly (${\tau }^{\left(E\right)}=0.05s$), whereas suppressive interactions extend over longer distances (${\sigma }^{\left(I\right)}=40\mu m$) and decay slowly (${\tau }^{\left(I\right)}=24.1s$) and thus integrate activity for long periods of time (Figure 3G–H). As a result, a single spike in one neuron strongly excites its immediate neighbours for a brief moment (<1 s), while slightly reducing the excitability of a large part of the tectal hemisphere for tens of seconds (Figure 3I).

Because our optimized network model could reproduce spontaneous activity with biologically realistic statistics, we reasoned that the recurrent interactions it encapsulates may provide useful insight into the relationship between recent neuronal activity and the internal state (excitability) of the tectal network. Therefore, we next explored the extent to which model-inferred network state could account for variability in stimulus-evoked activity.

Fluctuations in the state of excitability of neuronal populations caused by recent network activity is likely to contribute to variability in neural representations of sensory stimuli. To explore this, we applied our network model to recorded tectal activity data in order to estimate the incidental state of the biological network immediately prior to individual visual stimulus presentations and thereby predict the influence of ongoing activity upon stimulus-evoked responses.

We presented zebrafish with small moving spots, which evoke robust tectal activity (Figure 4A, Niell and Smith, 2005; Bianco and Engert, 2015). Spots swept horizontally across the visual field at two elevations, one slightly below and one above the horizon and to mitigate interactions between consecutive stimuli we initially used a long inter-stimulus interval (ISI, 5 min). The responses of tectal neurons displayed retinotopic mapping and direction selectivity (Figure 4B, Stuermer, 1988; Niell and Smith, 2005; Hunter et al., 2013).

Figure 4 with 1 supplement see all

Download asset Open asset

Figure 4—source data 1

Data provided as a MATLAB structure.

https://cdn.elifesciences.org/articles/78381/elife-78381-fig4-data1-v1.zip

Download elife-78381-fig4-data1-v1.zip

Using this data, we examined the relationship between response variability and model-estimated network state. First, we used our model to estimate the incidental state of OT neurons immediately prior to each stimulus presentation, based on the recent history of ongoing activity in the network. To do this, spikes inferred from calcium imaging data during the minute prior to stimulus onset (specifically, –60 to –5 s before stimulus onset) were spatiotemporally filtered and summed for each cell using the model parameters (Figure 4C, Materials and methods). In this way, we estimated the `linear drive` of biological neurons (treated here as equivalent to their excitability), in the same way as for the previous simulations, except that we used recorded, rather than simulated, activity.

Model-predicted excitability explained a substantial fraction of the trial-to-trial variability of visually evoked neural activity. Specifically, we found that for individual neurons, a threshold-linear function described the visual responses of the cell across single trials in terms of the pre-stimulus linear drive, estimated using the model, for those same trials (Figure 4D–E). Model-inferred linear drive explained response variability significantly better than a ‘null’ estimate of network state, using randomly sampled timepoints, or states estimated using the baseline model parameters (Figure 4E, Materials and methods).

These results suggest that the recurrent interactions implemented within our network model can predict incidental network state and thereby account for a fraction of the trial-to-trial variability of visually evoked tectal activity. Because the excitatory interaction is short-lived, this effect of recent ongoing activity is largely determined by spatiotemporal properties of the more prolonged inhibitory interaction. For example, the improved performance between the baseline and optimized model was largely due to the decrease in ${\tau }^{\left(I\right)}$ (from 180 to 24 s), with a small contribution from the increase in ${\sigma }^{\left(I\right)}$ (from 20 µm to 40 µm, Figure 4—figure supplement 1). In sum, recent network activity in the vicinity of a neuron tends to supress its excitability, resulting in a relatively reduced response to visual sensory input.

Next, we explored the reciprocal interaction, namely if sensory stimuli exert a lasting effect on the state of the network that consequently perturbs ongoing activity.

First, we used our model to predict the effects of visual input on network state. We expanded the model to incorporate external (sensory) input such that model neurons were activated according to the visual receptive fields of their corresponding tectal neurons. We estimated receptive fields by fitting Gaussian functions to the measured activity of each cell as a function of stimulus angle (Figure 5A, Materials and methods). Receptive fields collectively covered visual space, with an average width (2.5σ) of 38°±24° (mean ± SD, Figure 5B), comparable to previous reports (Niell and Smith, 2005). Because visual space is retinotopically mapped within OT, neurons responsive to high and low elevation stimuli occupied distinct anatomical locations (Figure 5C, Materials and methods).

Figure 5

Download asset Open asset

Figure 5—source data 1

Data provided as a MATLAB structure.

https://cdn.elifesciences.org/articles/78381/elife-78381-fig5-data1-v1.zip

Download elife-78381-fig5-data1-v1.zip

Running the network simulation with sensory input predicted a substantial, long-term reduction in tectal activity for tens of seconds following visual stimulation (Figure 5D). This prediction is concordant with the long inhibition time constant in the model (24.1 s) and, as a result of the spatially weighted recurrent connectivity, was biased to recently stimulated regions of OT. Thus, cells responsive to the recently presented stimulus were predicted to show a greater suppression of ongoing activity in the period following stimulus offset as compared to neurons responsive to the other elevation (Figure 5D). This spatial-selective suppression of activity was explained by a lasting state of low linear drive extending beyond the duration of the sensory stimulus (Figure 5F, Video 3).

Video 3

Download asset

To test this prediction, we repeated the same analysis with calcium imaging data from five fish. As predicted by the model simulation, neurons responsive to the recently presented stimulus showed a marked reduction of ongoing activity for tens of seconds, whereas neurons responsive to the other elevation showed a lesser suppression of post-stimulus activity (Figure 5E).

In sum, our network model can account for bidirectional interactions between neuronal activity and network state that lead to interactions between ongoing and stimulus-evoked tectal activity. Ongoing activity shapes the incidental state of neuronal excitability, which in turn modulates visual responses, while these responses exert a long-lasting suppressive effect upon network excitability and ongoing activity.

These observations lead to the prediction that visual stimuli will have a lasting effect on the state of the tectal network that may influence subsequent sensory responses. Thus, we next used network modeling and calcium imaging to explore how the recent history of sensory stimulation impacts visually evoked activity and behavior.

Based on the timescales and spatial ranges apparent in tectal dynamics, we designed an experimental paradigm incorporating a mixture of temporal and spatial relationships between visual stimuli. During 30 min ‘stimulation blocks’, small moving spots were presented at two different elevations (out of a possible three). One stimulus was presented at high frequency (‘common’, 30 s inter-stimulus interval, ISI) and the second at low frequency (‘deviant’, 5 min ISI). Each deviant stimulus was followed by a 90 s break. Stimulation blocks were separated by 30 min ‘rest blocks’, during which no stimuli were presented (Figure 6A).

Figure 6 with 1 supplement see all

Download asset Open asset

Figure 6—source data 1

Data provided as a MATLAB structure.

https://cdn.elifesciences.org/articles/78381/elife-78381-fig6-data1-v1.zip

Download elife-78381-fig6-data1-v1.zip

We first evaluated the model prediction for this experimental protocol by simulating it five times. To take into account adaptation of retinal inputs, we imaged the tectal neuropil in isl2b:GAL4;UAS:SyGC6s fish expressing GCaMP6s at the synaptic terminals of retinal ganglion cells (RGCs), and modulated the input gain for the simulated high-frequency stimulus accordingly (Figure 6—figure supplement 1 and Materials and methods). The model predicted an adaptation of tectal visual responses to the common stimulus, but not to the deviant stimulus (Figure 6A). For example, during the 2nd stimulation block, high and low elevation stimuli were presented as common and deviant stimuli, respectively (Figure 6A, blue and red). Adaptation reduced the mean response of model cells tuned to the high elevation stimulus to about 70% of their initial amplitude followed by a partial recovery during the 90 s breaks (Figure 6A, stars). This is due to a combination of gradual input (RGC) adaptation and additional activity-dependent suppression intrinsic to the tectum (Figure 6—figure supplement 1C–D). Conversely, the response to the low elevation stimulus remained approximately constant.

To test these predictions against experimental observations, we used the same stimulation protocol with five fish (Figure 6B). In agreement with the model prediction, high-frequency stimulus presentation caused a spatially specific response adaptation to the common stimulus. The adaptation pattern matches the predicted combined effect of RGC adaptation and intrinsic tectal dynamics (Figure 6A), rather than either one alone (Figure 6—figure supplement 1C–D).

To examine how this response modulation varies across cells, we compared the average response of individual cells during the 2nd stimulation block to their baseline response (specifically, their average response to low-frequency stimulation, large circles in Figure 6B). The response to the common (high elevation) stimulus was attenuated in most cells (Figure 6C–D, cyan), and attenuation magnitude was linearly correlated with the amplitude of the baseline response (Figure 6E). By contrast, response modulation to the deviant (low elevation) stimulus was less consistent (Figure 6C–D, yellow), although some attenuation was observed in cells that were strongly activated by the common stimulus (Figure 6E). These results show that the common stimulus drives a subset of tectal cells to a supressed state, which modulates their subsequent responses to either stimulus.

To conclude, OT displays experience-dependent response modulation as predicted by our network model. This effect allows for integration of past stimuli and suppression of subsequent responses to similar stimuli.

Lastly, we examined how activity-dependent changes in network state relate to visually guided prey-catching behavior, which is mediated by OT (Gahtan et al., 2005; Muto et al., 2013; Bianco and Engert, 2015; Helmbrecht et al., 2018). Because hunting routines of larval zebrafish invariably commence with eye convergence, we used the occurrence of convergent saccades (CS) to identify prey-catching behavioral responses (Figure 7A inset, Bianco et al., 2011; Trivedi and Bollmann, 2013; Bianco and Engert, 2015).

Figure 7

Download asset Open asset

Figure 7—source data 1

Data provided as a MATLAB structure.

https://cdn.elifesciences.org/articles/78381/elife-78381-fig7-data1-v1.zip

Download elife-78381-fig7-data1-v1.zip

We first analyzed ‘spontaneous’ CSs, which occasionally occur in the absence of visual stimuli (Bianco and Engert, 2015), to assess the link between ongoing tectal activity and behavioral output. By examining activity in 1 s time windows centred on CSs (‘CS-associated activity’, Figure 7A, orange shaded area), we found that spontaneous CSs are associated with an increase in tectal activity, particularly in the anterior tectum (Figure 7A–B, anterior tectum p=0.002, n=11 fish, Wilcoxon signed rank test, anterior vs posterior activity p=0.02, n=11 fish, Wilcoxon signed rank test). Moreover, CS-associated activity was lateralized in accordance with the direction of behavior. We showed this by comparing activity between the tectal hemispheres for spontaneous CSs that directed gaze to the left or right (>10° post-saccade gaze angle). In anterior (but not posterior) OT, CS-associated activity was significantly biased towards the hemisphere contralateral to post-saccadic gaze (Figure 7C–D, p=0.003 and 0.7 for anterior and posterior OT, n=408 CSs, two-way ANOVA). These findings are in agreement with previous reports linking anterior OT activity with prey-catching behavior directed to the contralateral side (Fajardo et al., 2013; Bianco and Engert, 2015).

These observations are compatible with intrinsic tectal dynamics contributing to the generation of spontaneous prey-catching responses. Alternatively, CS-associated activity could be a consequence of the eye movement (visual or proprioceptive sensory feedback), or represent a motor efference copy. Therefore, to investigate if tectal activity is likely to cause spontaneous behavior, we used our network model to estimate the state of tectal neurons (their linear drive) immediately prior to spontaneous CSs based on the recent activity history of the network (same method as above, Figure 4). This analysis showed that linear drive in anterior OT was significantly elevated prior to CSs (Figure 7E–F, pre-CS vs baseline p=0.004, n=11 fish, one-tailed t-test), and this difference was larger compared to posterior OT (p=0.04, n=11 fish, one-tailed paired t-test). At a single-cell level, pre-CS linear drive was positively correlated with CS-associated activity, both for individual CSs (Figure 7H) as well as for data pooled over all spontaneous CSs (Figure 7I). These results are compatible with ongoing activity in the tectal network modulating the excitability state of OT neurons and contributing to the generation of spontaneous prey-catching behavior.

Because we have shown that the recent history of visual stimuli has lasting effects upon tectal activity, we next asked how past experience of visual inputs modulates spontaneous and visually evoked prey-catching behavior. To do this, we used the mixed-ISI stimulation protocol (Figure 6), while tracking behavior.

When we analysed spontaneous behavior, we found that the rate of spontaneous CSs in the 2 min that followed any visual stimulus was reduced by approximately one-third compared to rest blocks (Figure 7J, n=13 fish, p=0.002, Wilcoxon signed rank test). This suppression of behavior is concordant with the lasting reduction in excitability state and ongoing tectal activity described above. Moreover, the probability of visually evoked prey-catching was modulated by stimulus history in a qualitatively similar way to visually evoked tectal activity. Specifically, response rates declined for common, but not deviant stimuli (Figure 7K, compare Figure 6B).

Taken together, these observations suggest that the state of the tectal network, reflecting the integration of previous sensory inputs as well as ongoing activity, impacts both spontaneous and visually evoked prey-catching behavior.

Motivated by the idea that ongoing activity is both a cause and consequence of changes in the excitability state of neural populations, we inferred a minimal recurrent circuit motif that was sufficient to generate stochastic activity sequences that closely resembled ongoing tectal activity. We achieved this by combining classical neural field models (Amari, 1977), where two opposing interactions are deterministic functions of intercellular distance, with the stochasticity of spiking Linear-Nonlinear-Poisson (LNP)-type models (Pillow et al., 2005; Pillow et al., 2008). Despite being defined by only seven parameters, this minimal motif supported the emergence of stochastic tectal bursting with biologically realistic statistics. The model also explained some of the variability in visually evoked activity, the effect of that activity on network state, visual response adaptation and spatially specific behavioral habituation, suggesting this recurrent interaction motif is a dominant feature of the optic tectum.

However, applying the same connectivity rule, with identical parameters, across every tectal neuron is certainly an oversimplification and our model should not be regarded as a detailed approximation of tectal connectivity. Indeed, although it qualitatively reproduced several biological observations, model predictions deviated quantitatively and more expressive models that accommodate cell-specific connectivity parameters are expected to provide more accurate quantitative predictions. Future extensions to the model, incorporating known heterogeneities in neurotransmission, morphology and intrinsic properties and with many more free parameters, should expand the accuracy and breadth of the model assuming such additions are adequately constrained by functional and/or anatomical data. In this regard, recent studies in zebrafish have characterised various aspects of molecular-genetic, morphological, and physiological diversity in OT (Helmbrecht et al., 2018; Antinucci et al., 2019; Shainer et al., 2022) and high-quality ultrastructural data is now available, allowing definitive reconstruction of synaptic connectivity (Svara et al., 2022). The model can also be extended to incorporate interconnections between OT and many other brain regions (e.g. nucleus isthmi, Henriques et al., 2019), guided by data from whole-brain activity recordings (van der Plas et al., 2023) and brain wide anatomical atlases (Marquart et al., 2015; Randlett et al., 2015; Hildebrand et al., 2017; Shainer et al., 2022).

Recent studies have suggested that bursting activity observed in OT represents recruitment of sub-networks with enhanced connectivity, or neuronal ‘assemblies’ (Romano et al., 2015; Avitan et al., 2017; Marachlian et al., 2018; Mölter et al., 2018; Diana et al., 2019). At the core of this model is the idea that to explain OT bursting there must be nonuniformities in connectivity in which cells forming an assembly are preferentially interconnected and operate as a local attractor network. However, here we show that biologically realistic localised bursting can emerge in a network in which, by design, there are no assemblies; rather the same connectivity motif (in which recurrent interactions smoothly decay with distance) is implemented identically for every cell. Although, as discussed above, our model is certainly and oversimplification and does not rule out the existence of assemblies, it highlights that localised bursting, by itself, is insufficient evidence for such heterogeneities in connectivity.

The recurrent interaction motif in our network model predicts that the activity of a given OT neuron strongly but transiently increases the excitability of its immediate neighbours, while causing sustained suppression of a broader population of neurons in its vicinity (Figure 3I). What biological mechanisms might underlie such interactions?

The short time scale we estimated for excitatory interactions (${\tau }^{\left(E\right)}=0.05s$) is compatible with direct, recurrent synaptic connectivity in OT. Such connections have been demonstrated by intracellular recording of delayed excitatory post-synaptic potentials in tadpole tectum following optic tract stimulation (Pratt et al., 2008). Delayed glutamatergic potentials were similarly recorded in slices from rat superior colliculus and are thought to mediate pre-saccadic bursting (Saito and Isa, 2003).

By contrast, the inhibitory interactions in our model caused activity-dependent suppression lasting tens of seconds (${\tau }^{\left(I\right)}=24.1s$), too long to be accounted for by direct ionotropic synaptic inhibition. This interaction might be mediated by pre-synaptic depression, which can persist for tens of seconds in the avian OT (Luksch et al., 2004). Alternatively, slow suppression might be mediated by activity-dependent changes in intracellular properties such as ionic concentrations or membrane conductances, potentially lasting for several minutes (Desai et al., 1999; Karmarkar and Buonomano, 2006; Gasselin et al., 2015; Zylbertal et al., 2017b).

Long-term suppression might also derive from activity in inhibitory populations external to OT. There is currently no evidence to suggest that tectal interneurons are involved in such suppression, since their visually evoked activity appears to be limited to transient responses (Preuss et al., 2014; Dunn et al., 2016). However, brain regions interconnected with OT, such as GABAergic populations in the vicinity of the isthmic hindbrain (Gebhardt et al., 2019; Henriques et al., 2019), might provide long lasting feedback inhibition.

Our results indicate that the state of the tectal network is predictive of prey-catching response probability, in agreement with the role of OT in generating premotor signals for this visuomotor behavior (Herrero et al., 1998; Gahtan et al., 2005; Bianco and Engert, 2015; Antinucci et al., 2019). What adaptive purpose might the network interactions we describe play during natural behavior?

The optic tectum (and its mammalian homolog the superior colliculus), is thought to be involved in detection and selection of salient inputs and generation of rapid behavioral responses (Boehnke and Munoz, 2008; Dutta and Gutfreund, 2014). The experience-dependent activity suppression we observe is compatible with this role, as it filters out expected inputs which are less likely to be associated with immediate danger or prey-catching opportunity. Similarly, zebrafish OT exhibits response adaptation to looming stimuli that is thought to play a role in the habituation of escape behavior (Marquez-Legorreta et al., 2019). Comparable stimulus-selective adaptation and behavioral habituation has been described for auditory stimuli in barn owl OT (Reches and Gutfreund, 2008; Netser et al., 2011), and for overhead looming stimuli in the mouse superior colliculus (Lee et al., 2020). Notably, adaptation in the barn owl is selective for multiple auditory stimulus features such as interaural time and level differences, as well as frequency and amplitude. In the simple model we propose, the strength of interactions between neurons depends solely on Euclidean distance and therefore the spatial proximity of their receptive fields. An interesting avenue for future study will be to examine whether recurrent interactions weighted by similarity of feature tuning can account for such feature-selective adaptation.

Zebrafish (Danio rerio) larvae were reared on a 14/10 hr light/dark cycle at 28.5 °C. For all experiments, we used zebrafish larvae homozygous for the mitfa^w2 skin-pigmentation mutation (Lister et al., 1999). For Ca2 +imaging experiments, we used larvae homozygous for Tg(elavI3:H2B-GCaMP6s)^jf5Tg (Vladimirov et al., 2014, ZFIN ID: ZDB-ALT-141023–2). For imaging RGC axonal projections, larvae were double transgenic for Tg(isl2b:Gal4)^zc60Tg (Ben Fredj et al., 2010, ZFIN ID: ZDB-ALT-101130–1) and Tg(UAS:SyGCaMP6s)^a155Tg (Dunn et al., 2016, ZFIN ID: ZDB-ALT-160406–1). All larvae were fed Paramecia from 4 dpf onward. Animal handling and experimental procedures were approved by the UCL Animal Welfare Ethical Review Body and the UK Home Office under the Animal (Scientific Procedures) Act 1986, under Home Office project licence 70/8162 awarded to Isaac Bianco.

For calcium imaging we used a custom-built digitally scanned light-sheet microscope. The excitation path included a 488 nm laser source (OBIS 488–50 LX, Coherent, Santa Clara, California), a pair of galvanometer scan mirrors (6210 H, Cambridge Technology, Bedford, Massachusetts) and objective (Plan 4 X, 4 x/0.1 NA, Olympus, Tokyo, Japan). A water-immersion detection objective (XLUMPLFLN, 20 x/1.0 NA, Olympus), a tube lens (f=200 mm), two relay lenses (f=100 mm) in a 4-f configuration, and sCMOS camera (Orca Flash 4.0, Hamamatsu, Hamamatsu, Japan) were used in the orthogonal detection path. For remote focusing (Fahrbach et al., 2013), an electrically tunable lens (ETL, EL-16–40-TC-VIS-20D, Optotune, Dietikon, Switzerland) was installed between the relay lenses, conjugate to the back focal plane of the objective. Volumes (375x410 x 75 μm) comprising 19 imaging planes spaced 4 μm apart, were acquired at 5 volumes/s. Each plane received laser excitation for 1ms (duty cycle 9%) resulting in average laser power at sample of 12.4 μW. To keep the observed population of neurons in each plane in focus throughout long imaging sessions, we implemented an automatic correction for slow drift in the Z direction. At the beginning of the experiment, we acquired two reference stacks, centred on two of the imaging planes, by incrementally biasing the Z scanning mirror and ETL in steps of 1.5% of their scan amplitude. During the course of the experiment, Z drift was estimated every 30 seconds by comparing recent images to these reference stacks, finding the reference images with the maximal XY cross-correlation, and averaging the two drift estimates. Every 5 min, the Z scan mirror and ETL signals were biased to offset any detected drift, according to the average of the ten most recent Z drift estimates.

For functional imaging, larval zebrafish were mounted in a custom 3D printed chamber (SLS Nylon 12, 3DPRINTUK, London, United Kingdom) in 3% low-melting point agarose (Sigma-Aldrich, St. Louis, Missouri) at 5 dpf and allowed to recover overnight before functional imaging at 6 dpf. Visual stimuli were back-projected (ML750ST, Optoma, New Taipei City, Taiwan) onto a curved screen forming the wall of the imaging chamber in front of the animal, at a viewing distance of ~10 mm. A colored filter (Follies Pink No. 344, Roscolux, Stamford, Connecticut) was placed in front of the projector to block green light from the collection optics. Visual stimuli were designed in Matlab (MathWorks, Natik, Massachusetts) using Psychophysics toolbox (Brainard, 1997). Stimuli comprised 10° dark spots on a bright magenta background, moving at 20°/s either left→right or right→left across ~110° of frontal visual space. Two or three elevation angles were used, calibrated for each fish during preliminary imaging by finding elevations separated by at least 15° that produce robust observable tectal activation (typically a very low elevation stimulus ~25° below the horizon, a low elevation stimulus ~10° below the horizon, and a high elevation stimulus ~5° above the horizon).

Eye movements were tracked during imaging experiments at 50 Hz under 850 nm illumination using a sub-stage GS3-U3-41C6NIR-C camera (Point Grey, Richmond, Canada). The angle of each eye was inferred online using a convolutional neural network (three 5x5 convolutional layers with 1, 1 and 4 channels, each followed by stride-2 max pooling layer, and a single fully connected layer), pre-trained by annotating images from multiple fish covering a wide range of eye positions. Eye movements were categorized as a convergent saccade if both eyes made nasally directed saccades within 150 ms of one another. Microscope control, stimulus presentation and behavior tracking were implemented using LabVIEW (National Instruments, Austin, Texas) and Matlab.

Functional 2-photon calcium imaging was carried out as described in Antinucci et al., 2019 to record tectal activity in a single imaging plane in the absence of visual stimulus presentation for 30 min. Prior to imaging initiation, the larva was allowed to adapt to the imaging rig for 15 min.

All calcium imaging data analysis was performed using Matlab scripts. Volume motion correction was performed by 3D translation-based registration using the Matlab function ‘imregtform’, with gradient descent optimizer and mutual information as the image similarity metric. A registration template was generated as the time-average of the first 10 volumes and then iteratively updated following each block of 10 newly registered volumes from the first 500 frames. This template was then used to register all remaining volumes. For elavl3:H2B-GCaMP6s experiments, 2D regions of interest (ROIs) corresponding to cell nuclei were computed from each template imaging plane using the cell detection code provided by Kawashima et al., 2016. For RGC axonal arbor imaging, two ROIs encompassing the tectal neuropil were manually defined for each imaging plane. The time-varying raw fluorescence signal F_raw(t) for each ROI was extracted by computing the mean value of all pixels within the ROI mask at each time-point. A slowly varying baseline fluorescence F₀(t) was estimated by taking the 10th percentile of a sliding 20 volume window, and was used to calculate the proportional change in fluorescence:

These values were subsequently zero-centred by subtracting the mean for each ROI, and ROIs with a slow drift in their baseline fluorescence (for which the standard deviation of mean-normalised F₀(t) was greater than 0.45) were discarded.

To estimate spike trains, the zero-centred $\frac{∆F}{F_0}$ time series for each cell was deconvolved using a thresholded OASIS (Friedrich et al., 2017) with a first-order autoregressive model (AR(1)) and an automatically-estimated transient decay time constant for each ROI.

To standardize the 3D coordinates of detected cell nuclei, template volumes were registered onto the Tg(elavl3:H2B-RFP) reference brain in the ZBB brain atlas (Marquart et al., 2017) using the ANTs toolbox version 2.1.0 (Avants et al., 2011) with affine and warp transformations. As an example, to register the 3D image volume in ‘fish1_01.nrrd’ to the reference brain ‘ref.nrrd’, the following parameters were used:

antsRegistration -d 3 -float 1 -o [fish1_, fish1_Warped.nii.gz] -n BSpline -r [ref.nrrd, fish1_01.nrrd, 1] -t Rigid[0.1] -m GC[ref.nrrd, fish1_01.nrrd, 1, 32, Regular, 0.25] -c [200×200 × 200×0,1e-8, 10] -f 12×8 × 4×2 s 4×3 × 2×1 t Affine[0.1] -m GC[ref.nrrd, fish1_01.nrrd, 1, 32, Regular, 0.25] -c [200×200 × 200×0,1e-8,10] -f 12×8 × 4×2 s 4×3 × 2×1 t SyN[0.1,6,0] -m CC[ref.nrrd, fish1_01.nrrd, 1, 2] -c [200×200 × 200x200×10,1e-7,10] -f 12×8 × 4x2×1 s 4×3 × 2x1×0

Following registration, tectal ROIs were labeled using a manually created 3D mask (Source data 1).

We used the same procedure to detect localised bursts of tectal activity in both recorded and simulated data. First, we found local maxima in tectal population activity after averaging the inferred spike trains of all tectal cells and smoothing with a sliding window (3 frames, 600ms). We next identified the population of cells that were active during a 1 s window centred on each peak and used the DBSCAN density-based clustering algorithm (Ester et al., 1996) to define spatial clusters based on Euclidean distances between the centroids of active cells. We excluded peaks during episodes of strong bilateral activation (commonly associated with vigorous swims or struggles) where more than 10% of the total tectal neurons were active and less than 70% of the active cells were located in the same hemisphere (16% ± 5% population peaks excluded per fish). We set the clustering distance threshold ϵ to 15 μm (~3 cell diameters) and the minimal number of cells in a neighbourhood, MinPts, to 12. MinPts was selected as the minimal number for which no clusters were detected after repeatedly performing circular permutations the time bases of all neurons, Figure 1—figure supplement 1D–E. After detecting the clusters associated with an activity peak, each was analysed independently to determine the initiation and termination times of the bursting activity of its constituent cells. A sliding window, six imaging frames in width (1.2 s) was evaluated at single frame increments. At each position, spike counts were discarded if they exceeded the 60^th percentile of a Poisson cumulative distribution function with λ equal to the mean number of inferred spikes across all cells in the cluster. Finally, the mean number of remaining spikes was computed and used to produce a vector describing the average Poisson-filtered firing within the cluster. Its first and last non-zero values were defined as the burst initiation and termination times.

To assess the correlation structure of ongoing tectal activity, we calculated the Pearson correlation of the spiking activity of each neuron (‘seed’) and the activity of all other tectal neurons. Spiking vectors were smoothed by filtering with a Gaussian window with σ=1.4 frames (280 ms). To generate correlation maps, we first projected the spatial coordinates of tectal neurons onto two dimensions: anterior-posterior and medial-lateral tectal axes. This was done by manually tracing a curve along the tectal anterior-posterior axis and calculating the distance of each neuron to this curve (its medial-lateral location) and the location along the curve where it passes closest to each neuron (its anterior-posterior location). Next, the two-dimensional correlation map of each seed (excluding the seed itself, but including all other neurons) was fitted to a general bivariate normal distribution density function:

where $A$ is the gain, b is the bias, ${\mu }_x$ and ${\sigma }_x$ are the A-P axis mean and standard deviation, ${\mu }_y$ and ${\sigma }_y$ are the M-L axis mean and standard deviation, and $\rho$ is the correlation coefficient.

Alternatively, the correlation map was fitted to an exponential decay function:

A random location within 100 μm of the seed cell was used as an initial guess for the peak ${\mu }_x,{\mu }_y$ . The following bounds were used when fitting:

The network model was based on the LNP stochastic spiking network formalism (Pillow et al., 2005; Pillow et al., 2008; Truccolo et al., 2005). Interactions between neurons (including feedback from a neuron onto itself) were modeled by temporal coupling filters and external (visual) input was considered instantaneous (see below).

The instantaneous firing rate of neuron j at time t, ${\lambda }_j\left(t\right)$, was derived by exponentiating its linear drive $\phi \left(t\right)$, or sum of inputs:

where ${\displaystyle x}$ is the external stimulus, $I_{ij}^{\left(E\right)}$ and $I_{ij}^{\left(I\right)}$ are excitatory and inhibitory coupling filters, $y_i$ are spike train histories of each neuron at time t, and μ is a global baseline log-firing rate (DC input).

The excitatory and inhibitory coupling filters were each determined by three global parameters: gain (g), spatial influence standard deviation (σ) and temporal exponential decay time constant (τ).

We also tested a model variant with exponential spatial decay function for connection weights (Figure 3—figure supplement 1A):

where X denotes excitatory (E) or inhibitory (I) interactions and $d_{ij}$ is the Euclidian distance between neuron i and neuron j in 3D standard coordinates following registration. If the two neurons reside in opposite tectal hemispheres, their coupling strengths were both scaled by 0.01.

In another model variant we eliminated 20% randomly chosen connections (both excitatory and inhibitory). To compensate for the reduced excitatory and inhibitory gain in the circuit we tuned the global gain parameters ($g^{\left(X\right)}$) accordingly.

The model simulation was implemented in Matlab (code available at https://github.com/azylbertal/TectalLNP copy archived at Zylbertal, 2023), based on code provided by Pillow et al., 2008. It uses a 50 ms timestep (resulting in 4:1 ratio with imaging rate) and comprises 14733 simulated neurons with locations inferred from the tecta of one fish. To significantly reduce memory and computational costs, a much longer (5 s) time step was used to simulate the slow inhibitory inputs to all cells. The inhibitory input was then interpolated before being summed with the excitatory input to produce the total linear drive. Spikes were drawn at random based on an inhomogeneous Poisson process with expected rates ${\lambda }_j\left(t\right)$ , which were updated according to the equations above whenever at least one neuron spiked. To compare the simulation results with imaging results, simulated spiking output was binned to 200 ms windows, equivalent to imaging volumes.

We tuned the seven model parameters ($g^{\left(E\right)}{,g}^{\left(I\right)},{\sigma }^{\left(E\right)}{,\sigma }^{\left(I\right)},{\tau }^{\left(E\right)},{\tau }^{\left(I\right)},\mu$) to reproduce observed bursting characteristics using the Matlab implementation of controlled elitist evolutionary multi-objective optimization with a genotype-based crowding distance measure (Kalyanmoy, 2001). Specifically, we used three optimization objectives: (a) the log-rates of bursts binned according to the logarithm of the number of participating neurons; (b) the log-rates of bursts binned according to the logarithm of the burst duration; (c) mean burst-triggered spiking activity (specifically, the average activity of each neuron triggered on the peaks of bursts it participated in and subsequently averaged across all neurons, smoothed by spline interpolation and scaled from zero to one according to its 5th and 95th percentiles).

In each generation of the algorithm, a population of (initially random) 200 models was evaluated in parallel using UCL’s Myriad computing cluster. Each simulation was run without external inputs for 306,000 steps (the equivalent of two hours and 15 min of imaging), and the initial 18000 steps (15 min) were discarded to avoid initialization transients. Localised bursts were detected based on the simulated spiking as they were detected for experimental data (see above). Since the inter-hemispheric interactions in the model are very weak (1% of the intra-hemispheric interactions), we were able to reduce computational cost during optimization by evaluating the models on one hemisphere. The loss value for each objective is the sum of squared errors between the simulation and experimental results. Parents for the next generation are selected by tournament: subsets of models are chosen by random (with replacement), and the best one from each set, that is the one with the highest ranking from the previous generation, is selected as a parent. 80% of the children are produced by intermediate crossover, where their parameters are chosen randomly from a uniform distribution bounded by the parameters of pairs of parents. The remaining 20% are produced by adaptive mutation, where they are randomly mutated with a direction that is maintained upon an increase in mean population fitness from the previous generation, or changed otherwise. The extended population (original population and children) is then evaluated and ranked based on Pareto front order, and within each front by genotype crowding distance (individuals with less close neighbors are ranked higher). Finally, the extended population is trimmed to retain the best 200 individuals for the next generation.

The evolution was stopped when the population scores stopped improving (~30 generations), producing a set of models on the approximated Pareto front – the solutions for which no individual model performs better than any other in all objectives. We chose the solution at the ‘elbow’ of the resulting front (the one closest to the origin in the z-scored loss space) and used it as a starting point for a local pattern search for final parameter refinement (Figure 3—figure supplement 1D–E).

We estimated the linear drive, or ‘excitability’ of each OT neuron prior to stimulus presentation (Figure 4) or spontaneous convergent saccade (Figure 7) by applying the spatiotemporal filters described by the LNP model to spiking history inferred from imaging data, regardless of whether other spontaneous events occurred during this time period (Figure 4C). For each event $i$, we calculated 60 s long vectors ${\phi }_{ij}\left(t\right)$ representing the sum of inputs to each cell $j$, based on the activity of all other OT cells during the –60 to –5 s period prior to the event, weighted by time and distance per the model parameters (Equation 1 above). Since the model temporal filters assign larger weights to recent spikes, we ignored the spikes detected during the final 5 s to account for the inherent uncertainty in inferring spikes from calcium imaging data. To account for cell-specific biases, such as spatial edge effects, we corrected ${\phi }_{ij}\left(t\right)$ by subtracting cell-specific ‘baseline’ mean linear drive vectors $\overline{{\phi }_j}\left(t\right)$, estimated by repeating the calculation for 60 epochs of ongoing activity. The reported linear drive point-estimate was calculated by averaging the final 10 s of the baseline-corrected vectors. For each neuron, we calculated the percentage of variance in activity across presentations of a given stimulus type explained by the estimated linear drive according to a threshold-linear fit (Figure 4D–E):

The simulated external input to each model neuron, $x_j\left(t\right)$ , was estimated from the mean recorded activity of the corresponding tectal neuron as a function of stimulus azimuth (${\lambda }_j\left(\alpha \right)$). To estimate receptive fields, the following function was fitted to this data:

producing values for gain ($A_j$), peak azimuth (${\mu }_j$) and width (${\sigma }_j$) for a specific neuron, stimulus elevation and motion direction (Figure 5A–B). To reduce the impact of noise, neurons with very low peak response or very narrow (${\sigma }_j$ <3°) receptive fields were excluded from receiving external input in the model.

Finally, to infer external input drive from measured spiking output, we scaled the gain and receptive field widths (to account for apparent expansion caused by recurrent connectivity) to obtain:

Our results are robust to substantial variations in the choice of these scaling factors.

We identified tectal neurons as responsive to one of the stimulus elevations (Figures 5–6) if their mean response firing rate exceeded the 90th percentile over the responses of all neurons and was at least five times higher than their baseline firing rate.

p-values were obtained from two-tailed t-tests unless otherwise noted. For the two-way ANOVA used to model lateral activity bias (Figure 7D), post-saccadic gaze direction (rightward or leftward) and fish identity were used as independent factors, without interaction terms. The reported p-value is associated with the gaze direction.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. For the purpose of Open Access, the authors have applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.

The authors thank members of the Bianco lab, James Fitzgerald, Eyal Wigderson and Avraham M Libster for helpful discussions and critical feedback on the manuscript and UCL Fish Facility staff for fish care and husbandry. AZ was supported by a BBSRC Discovery Fellowship (BB/S010564/1). IHB was supported by a Wellcome Trust and Royal Society Sir Henry Dale Fellowship (101195/Z/13/Z) and a Wellcome Trust Senior Research Fellowship (220273/Z/20/Z).

Animal handling and experimental procedures were approved by the UCL Animal Welfare Ethical Review Body and the UK Home Office under the Animal (Scientific Procedures) Act 1986, under Home Office project licence 70/8162 awarded to Isaac Bianco.

© 2023, Zylbertal and Bianco

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.