Functional brain imaging in larval zebrafish for characterising the effects of seizurogenic compounds acting via a range of pharmacological mechanisms

Functional brain imaging using genetically encoded Ca2+ sensors in larval zebrafish is being developed for studying seizures and epilepsy as a more ethical alternative to rodent models. Despite this, few data have been generated on pharmacological mechanisms of action other than GABAA antagonism. Assessing larval responsiveness across multiple mechanisms is vital to test the translational power of this approach, as well as assessing its validity for detecting unwanted drug‐induced seizures and testing antiepileptic drug efficacy.

consistency broadly mapping onto what is known about neurotransmitter-specific circuitry in the larval zebrafish brain. Many seizure-associated compounds also exhibited altered whole brain functional connectivity compared with controls.
Conclusions and Implications: This work represents a significant step forward in understanding the translational power of 4dpf larval zebrafish for use in neuropharmacological studies and for studying the events driving transition from small-scale pharmacological activation of local circuits, to the large network-wide abnormal synchronous activity associated with seizures.
3Rs, CNS safety pharmacology, drug discovery/target validation, functional neuroimaging, neuropharmacology, seizures, zebrafish 1 | INTRODUCTION Functional brain imaging using genetically encoded Ca 2+ sensors in optically transparent larval zebrafish offers an extremely powerful approach for analysing the in vivo response of the vertebrate brain to drug treatment. Despite its undoubted power, this approach has yet to be widely used in neuropharmacological studies. A research area that is beginning to exploit this methodology is the study of seizures and the chronic seizurogenic condition, epilepsy (Burrows et al., 2020). Traditionally, seizures have been classified as the culmination of abnormal neuronal excitation that progresses to hypersynchronous excitation across large-scale neural circuits. Most animal models of seizures and epilepsy rely upon unbalancing excitatory and inhibitory synaptic influences using pharmacological (or electrical) stimulation (Staley, 2015). Whilst this process is often exploited for anti-epileptic drug development (Kandratavicius et al., 2014), this can also occur through the inadvertent action of chemicals, including drugs (Ruffmann et al., 2006). The exact mechanisms that lead to seizures are not fully understood (Scharfman, 2007) and a major challenge is that, to understand the mechanistic basis of seizure generation and identify novel drug targets, studies of brain networks linking local and global scales are crucial (Rossi et al., 2018;Staley, 2015). Importantly, and in contrast with other current animal models, functional brain imaging in the larval zebrafish allows this and also provides a more direct indicator of neural activity than, for example, behavioural assessments of convulsions (Winter et al., 2008).
Consequently, a number of recent studies have focussed on functional imaging in the zebrafish larval brain after exposure to the model seizure-precipitant and GABA A antagonist pentylenetetrazole (PTZ), as well as a small number of other seizurogenic tool compounds (Diaz Verdugo et al., 2019;Liu & Baraban, 2019;Rosch et al., 2018;Turrini et al., 2017;Winter et al., 2017) or have demonstrated the application of this methodology in assessing the efficacy of antiepileptic drugs in genetic or PTZ-induced epilepsy models (Ghannad-Rezaie et al., 2019;Lin et al., 2018). The fact that functional brain imaging studies using larval zebrafish have focussed on one or two mechanisms of initiation, however, means that we still do not have a clear idea of whether neural network hyperexcitation induced by other drug types can be detected using this approach and what those neural imaging signatures look like (Burrows et al., 2020). Such information is critical both for assessing the validity of this approach for the detection of unwanted neuroactivity during drug development and its application for the discovery of new antiepileptic drugs.
Furthermore, understanding the neuropharmacological responsiveness of the larval zebrafish provides important validation data on the applicability of this model for more fundamental mechanistic studies of seizures and epilepsy, and on the translational power of the larval What is already known • Functional brain imaging in larval zebrafish offers a powerful approach for profiling neuroactive drugs.
• Translational power of this approach across multiple drug types has yet to be fully tested.

What this study adds
• We show that the larval zebrafish brain is responsive to a range of seizurogenic drugs.
• Distinct phenotypes are associated with seizurogenic compound exposure and with certain initiating mechanisms of action.

What is the clinical significance
• This approach appears appropriate for neuropharmacological profiling and for detecting the seizurogenic potential of drugs.
• This approach also has potential as an alternative in vivo model for studying epileptogenic mechanisms. zebrafish for use in imaging-focussed neuroscience studies more widely.
To address this knowledge gap, using four-dimensional functional brain imaging, we systematically tested the responsiveness of 4 days post fertilisation (dpf) larval zebrafish to 57 compounds with a known link to seizure generation in mammals and spanning more than 12 drug classes, alongside eight compounds with no such link. We opted to use 4dpf zebrafish as these are not considered protected vertebrates under European animal legislation and therefore ethically preferable, particularly for use in higher throughput screens where animal usage can be high. In addition, the use of 4dpf larvae allows unfettered functional imaging in the vertebrate brain in the absence of the confounding effects of a general anaesthetic. With this experimental approach, we sought to address these key questions:-(1) Are 4dpf larval zebrafish responsive to a range of molecular mechanisms linked to network hyperexcitation in mammals? (2) Are there phenotypic features that allow for differentiation between nonseizurogenic and proseizurogenic compounds and/or particular pharmacologies? (3) Are there brain regions that appear particularly important for the initiation and/or propagation of seizure-like network activity regardless of the initiating mechanism? (4) Are the resultant spatiotemporal patterns of neural activity consistent with our knowledge of neurotransmitterspecific brain circuitry in the larval zebrafish?

| METHODS
The experimental approach is detailed in the following sections and the process is summarised in Figure 1. Throughout, the data shown in graphs and tables are the means ± SEM of data points from individual animals (n). Statistical analysis was undertaken using larvae as the individual replicate values and only where the sample size was at least n = 5. Due to data loss, one dataset contained an n = 4 (picrotoxin at the lowest treatment level of 15 μM) and in this case, this treatment was not subject to statistical assessment. 'Outliers' were not excluded from the analyses or presentation of data. This manuscript complies with British Journal of Pharmacology recommendations and requirements on experimental design and analysis (Curtis et al., 2018).

Misha Ahrens (Janelia Research Campus, Howard Hughes Medical
Institute, Ashburn, Virginia, USA), were held in aquaria at the University of Exeter, at 28 ± 1 C under optimal conditions for spawning (12-h light: 12-h dark cycle, with 20-min dusk-dawn transition periods). Culture water consisted of mains tap water which was filtered by reverse osmosis (Environmental Water Systems UK Ltd.) and then reconstituted with Analar-grade mineral salts to a standard synthetic freshwater composition (final ion concentrations: 117 mgÁL −1 CaCl 2 .2H 2 0, 25.0 mgÁL −1 NaHCO 3 , 50 mgÁL −1 MgSO4.7H20, 2.3 mgÁL −1 KCl, 1.25 mgÁL −1 Tropic Marine Sea Salt, giving a conductivity of 300 mS). The culture water was aerated and heated to 28 ± 1 C before being supplied to individual zebrafish tanks on a semire-circulatory system. Water was routinely monitored for temperature, pH, conductivity, ammonia, nitrite and nitrate, all of which were maintained within appropriate limits for zebrafish. Spawning occurred at the time of lights on and was encouraged by provision of a spawning substrate. Fertilized eggs were collected shortly after spawning, bleached for 1 min in 1% w/v chloramine T (Sigma-Aldrich, Poole, UK) in culture water, rinsed with fresh culture water and transferred to Petri dishes which were filled with culture water and held

| Test compounds
The range of drug types assessed spanned agents targeting multiple mechanisms implicated in the initiation of neural network hyperexcitation (Easter et al., 2009) and comprised many compounds which have not, to our knowledge, been tested before in zebrafish.
The compounds tested encompassed those implicated in the induction of seizures in humans or nonclinical models via a range of molecular mechanisms and those with no known association with seizures ( Figure 2). Compounds were selected following an extensive review of available published data. This review allowed us both to identify compounds with a link to seizures in preclinical or clinical studies and to select appropriate pharmacological probes known to act via mechanisms implicated in the generation of neural network hyperexcitation.
The compiled list was not considered exhaustive but rather designed to encompass several examples acting through related mechanisms of action. Based on the strength of the supporting literature, we categorised the selected compounds as follows: category A tool compounds typically used in animal models of seizures or epilepsy; category B compounds for which there is strong evidence for an association with seizures from clinical or nonclinical animal data; category C compounds for which there is some evidence supporting a link with (non) clinical seizures; category D compounds for which there are limited reports of an association or where seizures are induced only under specific exposure scenarios such as during overdose (e.g morphine) or withdrawal (e.g. ethanol); category E compounds for which there is a clear mechanistic link but for which literature data on seizure induction are absent (typically pharmacological probes) and category F compounds for which no evidence or mechanistic association with seizure induction was evident (negative controls). The final classification was agreed a priori by all study participants based upon the information available following the reviewed literature.

| Selection of exposure concentration ranges and durations
The concentration ranges and exposure durations used in the functional imaging in 4dpf larvae were selected using information gained from a separate assessment of compound effects on locomotion in 7dpf larval, using a previously described standardised assay (see Winter et al., 2008). Bioanalysis of internal (whole body) compound concentrations were carried out in the same 7dpf animals using liquid chromatography with tandem mass. The results of this work on 7dpf larvae are summarised in Figure S1 and the methods used briefly outlined in the Supporting Information. The internal compound concentrations in 7dpf larvae only provide an estimate of actual internal compound concentrations in the 4dpf larvae, but collectively, these data helped both to select the initial exposure conditions used for the 4dpf functional imaging assessment (e.g., high uptake supported a short exposure period and lower concentration range and vice versa), and to interpret the resultant functional imaging data (e.g., where the absence of a biological effect might be explained by poor compound uptake).

| Drug treatment and light sheet microscopic imaging
We aimed to generate equally sized datasets of n = 8 animals per treatment in our exposures for which the imaging logistics required the treatments to be run over two separate days on two batches of embryos (i.e. from two separate spawnings). Each treatment was equally represented over the 2 days of experimentation and the experimental process undertaken identically. Although treatments were not blinded, each larva was randomly selected and assigned to test compound at one of three concentrations or to a fish culture system water (vehicle) control for the exposure durations detailed in Figure 2. Following exposure, each animal was imaged as described below. The order in which animals were exposed and imaged was rotated through all treatments on each day to remove any bias associated with the time of day used for assessment. The need to hold larvae motionless during imaging to avoid movement artefacts, whilst avoiding introducing the widespread confounding effects of using a general anaesthetic, necessitated the use of a neuromuscular blocker.
Tubocurarine was selected as a non-depolarizing neuromuscular blocker (nicotinic acetylcholine antagonist) with a stable and long duration of action (Appiah- Ankam & Hunter, 2004) and was applied at a concentration of 4 mM to induce rapid paralysis appropriate for the acute exposure regimen of our study. Although tubocurarine has been associated with histamine release and resultant hypotension clinically, importantly it exhibits no muscarinic receptor (mAChR) antagonist activity and is water soluble. Furthermore, it has been used to paralyse zebrafish in previous functional imaging studies (Favre-Bulle et al., 2018;Miyazawa et al., 2018;Tao et al., 2011;Vanwalleghem et al., 2020). Adding to this, we observed no evidence of the attenuated neuromuscular blockade during acetylcholinesterase (AChE) inhibitor application over the time scale we use for imaging after neuromuscular blocker application. Each larva was transferred to an aqueous mixture of compound and tubocurarine (4 mM) until muscle tone was lost (typically $1 min), then moved to tubocurarine and test compound in 1.4% low melting point agarose and finally drawn into a clear borosilicate glass capillary plugged with 2% low melting point agarose. No organic solvents were employed to minimise potential confounding effects and vehicle (water) control larvae were treated identically (including the presence of neuromuscular blocker) but without test compound. Larvae were positioned on the light sheet microscope ( Figure 1a) and dorsal-ventral optical sectioning undertaken in the horizontal plane (10 equidistant Z-slices across 220 μm taken in 1.875 s, repeated 200 times). Following imaging, larval health was confirmed through observation of a normal heart rate and blood flow.
In any functional imaging study there is trade-off between spatial coverage and temporal resolution. Our imaging conditions were optimised to obtain the fastest cycle time whilst maintaining whole brain volume coverage at an appropriate temporal resolution. Some previous studies have chosen to obtain higher temporal resolutions across single focal plains, however in order to fully capture the spatial component of neural network hyperexcitation and answer specific questions about regional importance in seizure generation, we opted for whole brain volume imaging at a mesoscale spatial resolution. The importance of this scale of interrogation in the study of seizures has recently been highlighted by Wykes et al. (2019) who pointed out that seizure foci do not 'operate in isolation', exhibiting connectivity with 'short and long range projections' and that mesoscale study is required in order to 'record propagation pathways', areas of the brain that are recruited during seizure development and to detect abnormalities between 'functionally connected but distant areas of the brain'. Similarly, Rosch et al. (2018) acknowledged that epileptic seizures are an 'emergent property' at the level of neuronal populations and mesoscale modelling may provide insights that assessments at the microscale do not. Accordingly, we used anatomically anchored image registration to provide mesoscale segmentation of the whole brain in 3D and which provided the ability to compare equivalent regions across different larvae and treatments using anatomical location as the reference point for comparison. Further to this, in initial studies we undertook in support of this work we identified that a 6-min assessment period was sufficient to capture seizure-like activity after drug treatment, whilst providing datasets of a manageable size (Winter et al., 2017).

| Light sheet microscopic imaging settings
The light sheet microscope was custom built using the Open Access platform for SPIM (Huisken, 2012;Pitrone et al., 2013) as previously detailed in Winter et al. (2017). The system comprised a 488 nm Argon laser (Melles Griot, Didam, Netherlands) for illumination, a 20×/0.5 NA objective lens (Olympus, Southend-on-Sea, UK) with an intermediary magnification of 1×, followed by a 525/50 nm emission filter and 495 nm bandpass filter (Chroma, Olching, Germany) to collect emitted light and a 5.5-MP Zyla sCMOS camera (30 frames per second, 640 × 540 pixels, 4 × 4 binning, 40-ms exposure, Andor, Belfast, Northern Ireland) to capture images. The camera and rotational axis (Picard Instruments, Albion, USA) were controlled though μManager and the OpenSPIM plugin (Pitrone et al., 2013). The Point Spread Function of our system was determined by performing Z-stacks of 50 × 0.5 μm fluorescent beads, embedded in 1% low melting point agarose within a glass capillary. Images were acquired at z-intervals of 1.5 μm using the same system settings employed during data acquisition (detailed above) and the PSF calculated using the MOSAICsuite point spread function tool (http://imagej.net/ MOSAICsuite). The spatial resolution of our system was measured using a segment from a multi-grid stage micrometer (Edmund Optics, UK) with 100-μm divisions. The grid was imaged in the sample chamber using the same optics and medium to determine the pixel length of each 100-μm division. Each 100-μm division at 1 × 1 binning was 304 pixels and at 4 × 4 binning was 76 pixels.

| Image processing and peak parameter measurements
Data extraction from resultant images was undertaken automatically using a custom Python image processing pipeline developed at the University of Exeter and as previously detailed in Winter et al.
(2017) (available at: https://gitlab.com/exeter-zebrafish-research/ zebrafish-neuroanalyser3d) ( Figure 1b). Briefly, each image was initially down-sampled (3x3x1) and then 3D-aligned against a labelled reference brain (Randlett et al., 2015) using registration against 14 user-selected key-points from a representative 3D image stack. This brain atlas was established for use in 5-7dpf larvae, so we undertook a full 3D affine transformation (optimised using the Umeyama algorithm) that enabled alignment of our image series' to the reference brain and allowed extraction of voxel intensity data from defined anatomical regions of interest selected to encompass only major (larger) brain structures. This brain atlas has been widely used for anatomical registration including for a number of previous functional imaging studies (e.g. Migault et al., 2018;Rosch et al., 2018;Vanwalleghem et al., 2020). During processing, images were automatically shift-corrected to account for any sample drift.
For this, the shape index function as defined by Koenderink and van Doom (1992) was used to simultaneously correct for x, y and zdrift using functions from the Scikit-image library (van der Walt et al., 2014). Next, quantitative data were extracted for each region of interest as the mean of incorporated voxels for each time point, to generate a temporal profile of activity for that region.
These profiles were baseline subtracted using a sliding minimum filter comprising 50 time-points or the signal length if less than 50 time-points and peaks in activity were subsequently identified and analysed (Figure 1c (i)). Fluorescence-intensity profiles then underwent fine-scale smoothing using Gaussian filtering (σ = 1.5) to remove high-frequency noise and activity peaks spanning 2 or more consecutive imaging cycles were identified by subtracting the baseline and applying a threshold of 2σ. From the identified activity peaks in the profiles for each region of interest, time-averaged peak height and Area Under the Curve (AUC-calculated as the sum of the intensity values within each defined peak region) were calculated for each animal (referred to as 'peak parameters') as these two metrics were found to be the most useful for defining activity phenotypes (Figure 1c (ii)). We have previously described the performance characteristics of our imaging process (Winter et al., 2017). Briefly, the error of the 3D affine registration processes was calculated as 27.56 ± 0.43 μm, the theoretical point spread function was 0.751 μm, the measured lateral resolution was 0.76 pixels per micrometre and the temporal resolution was 1.875 s per full brain volume.

| Data analysis-Describing the resting state network
To describe the resting state network of 4dpf larvae analysed using this methodology, analysis of the peak parameters across all vehicle (water) control groups was undertaken. This analysis allowed us to describe the relative levels of neural activity in different brain regions under 'resting' conditions to establish which regions were comparatively (in)active during imaging. When comparing between regions of interests under such circumstances, it is important to control for potential regional variations in elavl3:GCaMP6s expression and, therefore, data were first normalised by subtracting each voxel's baseline, estimated using a sliding minimum filter (comprising 50 time-points or the signal length if less than 50 time-points) and then dividing by this baseline value. Functional connectivity in vehicle control larval brains was also assessed using the method detailed below.
2.8 | Data analysis-Analysis of region-specific GCamp6s fluorescence across the whole larval brain between regions of interests in treated animals, thus avoiding potential bias associated with regional differences in GCaMP expression in these cases.

| Data analysis-Multivariate exploration of larval brain functional connectivity
Functional connectivity in the context of neuroimaging pertains to the relationship between the neuronal activity patterns of anatomically separated brain regions, which in turn reflects the level of communication between these regions (van den Heuvel & Hulshoff Pol, 2010).
Practically, these relationships can be represented by measures of correlation that reveal information on mesoscale or macroscale networks responding to stimulation, in this case mediated by the action of seizurogenic compounds. Here, multivariate analysis of functional connectivity phenotypes was hypothesised to reveal information about the activation or suppression of specific neural circuits following compound treatment, by reflecting the spatial distribution of targets across the brain.
First, a connectivity matrix was generated for each animal by calculating the Pearson's correlation coefficient (r) between the time series' of all possible pairs of regions of interests (Figure 1d (i)). Initial assessment of concentration response curves showed that in the majority of cases the maximal response of larvae (versus control data) to compound treatment was apparent at the highest or second highest treatment concentrations. Consequently, it was decided to use these two treatment levels in the exploratory multivariate analysis. The lowest treatment concentration data were not included as in most cases activity levels were indistinguishable from controls and as such would likely reduce the ability to distinguish the phenotypic effect of drug treatment.
The data for DMCM, domoic acid, Cis-ACPD, TEA, bethanechol, theobromine and paraxanthine were also excluded from these analyses as we found no evidence of compound uptake (see Section 3).
Multivariate analysis required a consistent number of regions of interests across all experiments and therefore, regions of interests that were inconsistently registered (i.e. were absent in >2% of larvae analysed) during image processing as they occurred at the limits of Z-plane penetration were excluded prior to analysis, leaving 39 analysed regions of interests per animal. All values from each connectivity matrix were then reshaped into a connectivity vector (generating 741 regions of interest connections per larva) and the mean of the vectors for all animals calculated to provide one vector per compound (Figure 1d (ii)). To allow visualisation of similarities or differences between resultant profiles and to generate groups on which post hoc statistical analysis could be undertaken, a simple unsupervised agglomerative hierarchical cluster analysis ('unweighted pair group method with arithmetic mean' or UPGMA) was undertaken. This was considered the best approach where no prior information on the relative importance of any particular data feature was known. For this, the Euclidean distances between vectors of all of the compounds were calculated to provide a value for similarity between each compound within the multivariate feature space,and UPGMA cluster analysis undertaken. The same analysis was also undertaken on all control animal data to provide confidence that the observed clustering of 'related' compounds was not driven by potential differences in the baseline levels of GCaMP fluorescence between experimental runs ( Figure S2).
2.10 | Data analysis-Post hoc assessment of phenotype-defining data features Next, to assess whether the groupings identified by cluster analysis were based upon phenotypic similarities amongst grouped compounds and to test the hypothesis that certain data features differ between primary pharmacological mechanisms of action, post hoc statistical comparisons between compound-treated larvae and the resting state network were undertaken (Figure 1d (ii)). For this, Mann-Whitney U tests with a highly conservative Bonferroni correction (α = .000055) (due to the high number of comparison undertaken on each matrix) were used to identify statistically significant differences between compound-treated and control functional connectivity phenotypes, by comparing the r values obtained for all pairs of regions of interests. A difference between these two values would suggest an increase or decrease in the strength of that connection between the two groups, and help to determine whether any brain regions or functional connections were particularly prevalent within a specific grouping, and provide more confidence in the significance of groupings identified from the cluster analyses.

| Data analysis-Overall functional connectivity as an indicator of altered network synchronicity
Traditionally, a key defining feature of seizures is the occurrence of hypersynchronicity between disparate brain regions. It has, however, been acknowledged that the interictal-ictal cycle is likely to incorporate aspects of both hypersynchronicity and hyposynchronicity depending on the point in the cycle at which assessment is made and the spatiotemporal resolution of the detection methods used (Jiruska et al., 2013). With this in mind and given there is currently no universally accepted classification of epileptiform brain activity within Ca 2+ imaging datasets (Burrows et al., 2020), here we compared compound-treated and resting state network average functional connectivity (both average r values and numbers of connections with an r value >0.7) using Mann-Whitney U tests (p < .05) (Figure 1d (iii)).
These data were taken as an indicator of the degree to which synchronisation between regions was occurring compared with that observed in the resting state network. It was hypothesised that an increase in the number of functional connections and/or the strength of those connections may allow differentiation between nonseizurogenic and pro-seizurogenic compounds acting via a range of primary initiating mechanisms. Analysis of average functional connectivity was undertaken across all brain regions combined and also separately on the average of those connections that involve only regions within the telencephalon, diencephalon, mesencephalon, or the rhombencephalon, in order to assess changes in macroscale functional connectivity.

| Materials
All test chemicals are detailed in full, along with their Chemical Abstracts Services (CAS) registry numbers, in Figure 2 and its accompnaying legend. All test chemcials and reagents were obtained from Sigma-Aldrich (Gillingham, UK) or Tocris (Bristol, UK).

| The resting state network involves activation of hindbrain autonomic, sensory and motor centres
Analysis across all vehicle (water) control larvae was undertaken to describe the resting state network of animals analysed using our experimental setup. Analysis of the peak parameter data revealed low-magnitude, high-frequency neural activity, most frequently in hindbrain regions functionally associated with homeostatic control, sensory processing, arousal and motor activity (Figure 3a). Analysis of functional connectivity (Figure 3b) further supported this with the strongest functional connections and most connected regions also predominantly associated with these physiological processes.

| The 4dpf larval zebrafish is responsive to a range of molecular mechanisms known to induce neural network hyperexcitation in mammals
Peak parameter analysis revealed that many compound-treated larvae exhibited increases (net regional activation) or in some cases decreases (net regional inhibition) in GCaMP6s fluorescence in specific anatomical regions, when compared with the corresponding vehicle (water) controls (Figure 4). In the majority of cases, concentration-response relationships were also evident supporting compound-induced modification of neural activity ( Figure S3).
Exposure to GABA A receptor, glycine receptor, AChE and NA/5-HT/ DA modulating compounds generally resulted in a higher number of changes in activity within specific brain regions when compared with other pharmacological classes. Also of note, neural activity in hindbrain regions and in particular the cerebellum and its sub regions, the corpus cerebelli and valvula cerebelli, were more frequently altered after exposure to compounds associated with seizure induction ( Figure S4). In contrast, a number of potentially seizurogenic compounds showed no alteration of fluorescence intensity versus the water controls in any brain region. Of these DMCM, domoic acid, Cis-ACPD, TEA, bethanechol, theobromine and paraxanthine showed no observable effects in the separate maximum tolerated concentration or locomotor assessment of 7dpf animals and no differences in GCaMP6s fluorescence peak height or AUC. There was also no evidence of uptake for these compounds (<2% of the external concentration or not measurable) in 7dpf animals suggesting a lack of internal exposure might explain the inactivity observed in 4dpf larvae here (compound name highlighted in light grey in Figure S1). Low internal exposure, however, was less likely to explain an absence of altered peak height or AUC after exposure to 4-AP, PF-06767832, SNC80, meperidine, physostigmine, maprotiline, cocaine and theophylline. Exposure to clonidine and quinine, two compounds not normally associated with CNS excitation in mammals, conversely resulted in significantly elevated GCaMP-fluorescence in certain brain regions, especially in the case of clonidine When the number of brain regions exhibiting an increase or decrease in peak height or AUC were assessed (AUC data in Figure S5), strongly seizurogenic compounds generally showed more regions with elevated activity and fewer with reduced activity compared with the nonseizurogenic compounds. Notable exceptions were clonidine and quinine, as mentioned, and the 'positive' compounds fentanyl, carbachol and amitriptyline, which showed reductions in activity across multiple brain regions. Weaker seizurogenics and those categorised as showing seizure liability only under specific circumstances (category D in Figure 2) tended to locate away from the two extremes showing fewer significant changes in brain activity per se.

| Similarities in functional connectivity phenotypes of larvae exposed to monoaminergics and cholinergics
The functional connectomes resulting from exposure to each compound are shown in Figure S6 with an overview of published data on target distribution summarised in the Supporting Information. In addition, functional connections showing a difference compared with the resting state network are shown in Figure S7.
Cluster analysis revealed a number of groupings in which similar pharmacologies co-located ( Figure 5), most notably containing close groupings of monoaminergics (groups 3 and 8 in Figure 5) and cholinergics (both mAChR agonists and AChE inhibitors, groups 4 and 5). The majority of drugs that inhibit dopamine (DAT) norepinephrine (NET) and 5-HT (serotonin; SERT) transporters located in group F I G U R E 3 Summary of peak fluorescence intensity and functional connectivity data from vehicle (water) control larvae tested across 65 experimental datasets. Panel (a) shows the peak fluorescence intensity data per brain region of interest (ROI) as the average baseline-corrected fluorescence intensity (Δf/f = (f 1 − f 0 )/f 0 * 100, where f 1 = peak fluorescence intensity and f 0 = baseline fluorescence intensity within each ROI) of larvae within each control group (per column), presented as the % of the maximum peak intensity shown within that experimental dataset (n = 7-9 larvae per group). The most active brain region in each dataset is represented by 100%, with decreasing red/orange colour intensity representing ROIs with lower levels of comparative activity. Note the high consistency across experiments in brain regions exhibiting elevated basal activity, especially those that showed a mean peak intensity of >33% higher basal activity than the overall mean across all regions: intermediate hypothalamus, area postrema, inferior olive, locus coeruleus, Mauthner cells, medial vestibular nucleus, inferior raphe, spinal cord and spinal cord neuropil. Panel (b) shows the resting state functional connectome of control larvae presented as the average across 65 experimental datasets along with a table summarising the 10 most connected brain regions and the 10 strongest functional connections (n = 516 larvae). Abbreviated brain regions are as follows: area postrema (AP); cerebellum (C); corpus cerebelli (CC); dorsal thalamus (DT); eyes (E); eminentia granularis (EG); eminentia thalami (ET); habenulae (H); intermediate hypothalamus (I); inferior olive (IO); interpeduncular nucleus (IPN); lobus caudalis cerebelli (LC); locus coeruleus (Lco); lateral reticular nucleus (LRN); Mauthner (Mt); medial vestibular nucleus (MVN); noradrenergic neurons of the interfascicular and vagal areas (NIV); spinal cord neuropil (Np); olfactory bulb (OB); olfactory epithelium (OE); pineal (P); preoptic area (PA); pallium (pa); pretectum (PT); raphe inferior (RI); rostral hypothalamus (Ro); raphe superior (RS); spinal cord (SC); subpallium (Sp); posterior tuberculum (TB); torus longitudinaris (TL); tegmentum (tm); tectum neurophil (TN); torus semicircularis (TSC); tectum stratum periventriculare (TSP); tangential vestibular nucleus (TVN); Valvula cerebelli (V); vagal ganglia (VG); ventral thalamus (VT). Note: Positions of regions approximate to allow spacing of nodes and visualisation of connections 3, along with phencyclidine and two μ receptor agonists (fentanyl and meperidine). In contrast with all other antidepressants, amitriptyline exhibited a reduced activity across almost all brain regions versus the control state (Figure 4), despite widespread modulation of regional functional interconnectivity ( Figures S6 and S7, panel   33). The grouping of several muscarinic agonists and AChE inhibitors together also indicated a common phenotype amongst some cholinergics, although other class members located elsewhere (e.g. PF-06767832 and Donepezil, Figures S6 and S7, panels 16 and 23, respectively).
Notably, the GABA A receptor modulators ( Figure S6  3.4 | Alterations in brain region functional connectivity after exposure to cholinergic and monoaminergic compounds Figure 5b shows the number of functional connections for each brain region that differed after compound treatment, compared with the same connection in the resting state network. Furthermore, the blue and red colours, respectively, indicate whether the majority of those connections exhibited an increase or decrease in the r value compared with the resting state network and thus whether these connections were strengthened or weakened by compound treatment. The cholinergics generally showed a reduction in the strength of connectivity F I G U R E 4 Summary of peak GCaMP6s fluorescence intensity data measured in brain regions of larvae exposed to a range of neuroactive test compounds. Column 1 (left) shows the primary pharmacological mechanism of action; column 2, the abbreviated compound name; column 3, the strength of association with seizure induction in mammals according to the categories in Figure 2. Data (column 4 onwards) are expressed as % increase over the same regions of interest (ROI) from the corresponding water control group (ΔfT/fC = (f 1 − f 0 )/f 0 * 100, where f 1 = treated group and f 0 = control group, mean fluorescence peak height) (n = 6-10 larvae per group). Colour coding signifies the magnitude of effect with increasing red colour intensity signifying an increase, white no change and blue a decrease compared to the control. Within each box, statistically significant differences between compound-treated and control groups are shown (p < .017): + signifying a significant increase at one concentration, ++ at two and +++ at three. -, --and ---signifies a significant decrease for the same categories. O signifies an overall difference only, following Kruskal-Wallis analysis. Abbreviations are as follows: GABA = γ-aminobutyric acid; GLU = glutamate; ACh = acetylcholine; P1/PDE = adenosine/ phosphodiesterase; NA/DA/5-HT = noradrenaline/dopamine/5hydroxytryptamine (serotonin); Kv = potassium channels; OR = opioid receptor; αR = α-adrenoreceptor; H = histamine; CYP = cytochrome P450) across most brain regions, with the most commonly connected regions being the tectum stratum periventriculare, the locus coeruleus, medial vestibular nucleus and the inferior raphé. Donepezil, in stark contrast showed a widespread and large increase in functional connectivity across virtually all brain regions, whilst muscarine and PF-06767832 showed no effect. The main grouping of F I G U R E 5 Results of the hierarchical cluster and post-hoc analysis of the functional connectivity data. Panel (a) shows the dendrogram resulting from the cluster analysis on the Euclidean distances between mean r values for each of the compounds assessed. Identified groupings showing a mean Euclidean distance of no more than 3 (represented by the red dotted line in the dendrogram) are shown in the right hand table denoted by borders, along with the corresponding primary pharmacological mechanism of action and seizure liability category (as defined in Figure 2). Grey boxes represent the closest sub groupings. Numbers in the right hand column are references for each defined grouping that occurs in the corresponding dendrogram. Note in particular, the grouping of specific compounds that share a primary mechanism of action (e.g., especially groups 3, 4 and 5). Panel (b) shows the results of the post hoc assessment of changes in functional connectivity of each brain region in compound-treated animals, versus the same region in the resting state network. Abbreviated compound names, along with the primary pharmacological mechanism of action, seizure liability category and dendrogram group number are shown along the top of the table, with each registered brain region shown down the left hand side of the table. Each cell shows the number of connections to that region for which the mean r value had significantly changed versus the equivalent connection in the resting state network. Those boxes coloured red showed a majority of increased r values, blue a majority of decreases and orange an equal number of increases and decreases (see Figures S6 and S7 for all compound connectomes). These data reveal the brain regions in which functional connectivity was most significantly changed after compound treatment. Comparisons were undertaken using Mann-Whitney U tests with a Bonferroni correction applied (significance at p < .00006). A blank cell signifies there was no change in functional connectivity with any other region when compared to the resting state network. n = 12-19 larvae per drug exposed group and n = 516 for the resting state network monoaminergic drugs (group 3) generally exhibited increased r values compared with the resting state network suggesting an increase in the strength of functional connections after compound treatment.
Amongst these, the most consistently connected regions were in the hindbrain, except for the tegmentum and the intermediate hypothalamus which were also prominently connected within the close pairing of amoxipine and nomifensine (group 8). Amitriptyline was unusual amongst the monoaminergic reuptake inhibitors in that there was widespread reduction in functional connections particularly from posterior hind brain regions, although many of the most interconnected regions were in common with other mood enhancing drugs.
3.5 | Differences in patterns of brain functional connectivity between most nonseizurogenic and proseizurogenic compounds Many of the more active seizurogenic compounds exhibited an increase in the average number of connections and/or in the average r value per connection compared with the resting state network ( Figure 6), suggesting a tendency towards excitatory synchronisation between brain regions. Moreover, the weaker acting pro-seizurogenic and 7/8 nonseizurogenic compounds showed no altered network activity with respect to either metric used. Assessment of macroscale functional connectivity ( Figure S8) also suggested that seizurogenic compounds were more likely to exhibit changes in mesencephalic and, to a lesser degree, diencephalic, functional connectivity compared with nonseizure-associated compounds. Interestingly, however, the well-known seizurogenic compounds pilocarpine and bicuculline showed reductions in the average number and strength of functional connections.

| DISCUSSION
Analysis of the resting state network suggested stimulation from the imaging environment, with elevated activity apparent in regions involved in motor coordination and control (Heap et al., 2013), sensory processing (Yokogawa et al., 2012) and physiological homeostasis. Collectively, this mirrors other studies reporting activation of sensory-motor networks in conscious, immobilised larval zebrafish following external stimulation and subsequent functional imaging (Ahrens et al., 2012;Migault et al., 2018).
Analysis of GCaMP-indicated brain activity suggested that 4dpf zebrafish larvae are responsive to a wide range of pharmacologies associated with the induction of neural network hyperexcitation in mammals. The mechanistic classes that showed the most consistent induction of neural activity were the NMDA receptors and GABA A receptor agonists, glycine receptor agonist, AChE inhibitors and NA/5-HT/DA reuptake inhibitors. Although not a definitive measure of brain exposure in 4dpf larvae, the data for whole body compound uptake in 7dpf larvae provided an indication of where low compound uptake may be the reason for a relative lack of biological responsiveness (e.g. in the case of DMCM, domoic acid, Cis-ACPD, TEA, bethanechol, theobromine and paraxanthine). This highlights the importance of having some information on relative compound penetration for the interpretation of zebrafish chemical screens and for identifying false negatives that may occur as the consequence of low compound uptake (Berghmans et al., 2008;Cassar et al., 2020).
Another important observation is the frequency with which the cerebellum and related regions exhibited elevated neural activity after treatment with seizurogenic compounds. These changes in cerebellar subregions are likely reflective of activity in the region as a whole as full differentiation of cerebellar cellular layers is not thought to be complete until around 6dpf (Volkmann et al., 2008). Nevertheless, the frequent activation of these regions after exposure to seizurogenic compounds supports a role for cerebellar circuitry in neural network hyperexcitation. Previous work in mammals has suggested the cerebellum may play an important role in seizures and epilepsy (Rijkers et al., 2015), and this structure is a recognised target for deep brain stimulation to treat refractory epilepsy (Fisher, 2013). Furthermore, activation of the cerebellum has been reported in previous studies of PTZ-induced seizures in zebrafish (Baraban et al., 2005)  The absence of any changes in brain activity after exposure to several other seizurogenic compounds that showed evidence of uptake suggests that pharmacological insensitivity could be a factor.
Within-class individual compound-insensitivity could be due to subtle interspecies differences in receptor subtype homology, differences in the level of local tissue distribution, or temporal methodological insensitivity. For example, the larval zebrafish δ receptor has previously been shown to have a low affinity for SNC80 (Sanchez-Simon et al., 2012), PF6767832 is selective for the M 1 receptor alone and larval zebrafish have previously been shown to be insensitive to cocaine perhaps due to differential tissue partitioning resulting in reduced central stimulatory activity (Kirla et al., 2016). 4-AP, a wellknown seizurogenic compound, has previously been shown to induce generalised CNS hyperexcitation, but without consistent induction of seizures (Liu & Baraban, 2019). Moreover, Liu and Baraban (2019) reported that fast spiking was detected using electrophysiological, but not Ca 2+ imaging approaches and that 4-AP exposure results in dispersed small activated neuronal clusters. This perhaps suggests that the mesoscale segmentation and temporal resolution of our system is not sufficient for the consistent capture of fast ictal events induced by 4-AP and certain other compounds; therefore, faster, higher resolution systems may be needed for full neuropharmacological functional imaging coverage. In contrast, exposure to the centrally acting α 2 -adrenoceptor agonist clonidine resulted in a surprisingly high degree of CNS activation. Analysis of 7dpf larval locomotor activity, however, revealed a biphasic response comparable to that following Ca 2+ imaging and such a dose-response relationship has been reported previously in relation to anxiolysis versus anxiogenesis in rats (Söderpalm & Engel, 1988). Furthermore, although clonidine is a putative anticonvulsant, there are very limited data suggesting seizure induction under specific circumstances (Feron et al., 2008). Quinine F I G U R E 6 Assessment of average functional connectivity across the whole brain. Panel (a) shows a comparison of the functional connectivity across the whole brain of larvae exposed to each compound compared with the resting state network animals. Data are shown as the % difference in the number of functional connections between the two groups (where r > .7) along with the corresponding p value generated after Mann-Whitney U testing, along with the same parameters for the average r value of the two groups. Compounds are listed in order of the % size difference from the highest to lowest with values highlighted in red signifying increases and those in a blue decreases compared with the control value. p values highlighted in pale yellow were statistically significant (p < .05). Note that most compounds with no known association to seizure showed no change in functional connectivity versus the resting state network, whereas the stronger seizurogenic compounds (especially the precipitants) generally showed an increase in functional connectivity. Note also that some known seizurogenic compounds also exhibited a decrease in functional connectivity compared with the resting state network. Panel (b) shows representative connectomes from compounds across the full range of values from those showing a large increase in functional connectivity (top), to those showing a large decrease in functional connectivity (bottom). Data shown are connections that exhibited a different mean r value compared with the equivalent connection as a mean of all the vehicle (water) control animals tested. Nodes connected with a red line represent significant increases and blue significant decreases in the r value compared with the control. Nodes labelled in green showed 1 or more significantly different connections, in blue, 5 or more and in red 10 or more. Comparisons here were undertaken using Mann-Whitney U tests with a Bonferroni correction applied (significance at p < .00006). n = 12-19 larvae per drug exposed group and n = 516 for the resting state network. Abbreviated brain regions as detailed in Figure 3. Note: Positions of regions approximate to allow spacing of nodes and optimal visualisation of connectivity exposure also resulted in elevated CNS activity in a small number of brain regions perhaps due to its known blockade of several K + channels (Zou et al., 2018). The most strongly activated region was the medial vestibular nucleus, which is of particular interest given the known potent inhibitory effect of quinine on mechanosensory hair cells of the lateral line (Thomas et al., 2013) and may point to compensatory activation linked to vestibular inhibition in these animals.
Most of the D-classified compounds showed relatively low levels of CNS activation or suppression with the exception of nomifensine, galantamine and caffeine. Nomifensine is a NA/DA reuptake inhibitor which has been linked with seizures in the clinic (Edwards & Glen-Bott, 1987) and is known to increase locomotion in mice (Kitanaka et al., 2012) supporting the potential for neural-excitation.
Similarly, caffeine predominantly exerts its CNS stimulatory effects through adenosine receptor antagonism and has a well-established association with seizures but only at high dose levels (van Koert et al., 2018).
Amitriptyline and carbachol rather unexpectedly exhibited decreases in activity across multiple brain regions. Of these, amitriptyline also showed an inhibitory effect on larval locomotion and very high compound uptake. A sedative-like effect has previously been seen in adult zebrafish exposed to amitriptyline (Demin et al., 2017) and its known sedative effects in mammals are likely driven by activity at the histamine H 1 receptor (Lawson, 2017). Similarly, carbachol is also known to be a sedative in rats, an effect likely driven by action on the M 2 receptors of the brainstem (Ma et al., 2001).
To date, functional connectivity analysis of zebrafish neural activity has been limited to few studies. Rosch et al. (2018) for example used dynamic causal modelling to reveal the optic tectum may be central in driving network wide synchronisation and Diaz Verdugo et al. (2019) used pairwise correlation analysis (as used here) to show that glia may play an important role in the transition to a seizure. Our assessment of functional connectivity served two purposes:-first, we hypothesised that functional connectivity phenotypes would reflect mechanism-specific circuit activation or inhibition in 4dpf larvae and thus reveal differences between mechanisms of action, and second, that changes in overall functional connectivity would support a trend for the mesoscale and macroscale hypersynchronisation characteristic of seizure induction after activation via a range of pharmacological mechanisms of action. Our functional connectivity data revealed that larvae exposed to certain cholinergics and monoaminergic reuptake inhibitors showed common features within their phenotypes. Notably, both muscarinic agonists and AChE inhibitors grouped together suggesting a common cholinergic phenotype. The most prominent regions showing altered connectivity after exposure to cholinergics were the tectum stratum periventriculare, locus coeruleus, medial vestibular nucleus and the inferior raphe. The development of the cholinergic system in the zebrafish CNS was described by Arenzana et al. (2005), who reported that spatiotemporal development from 2dpf to 5dpf followed a caudorostral gradient. Our data also support preferentially altered functional connectivity in hind brain regions after treatment with cholinergic agents but of course reflect wider synaptic connectivity, as exemplified by the widespread activity of donepezil, rather than discrete distributions of cell bodies or subcellular receptor populations. Also of note was the close grouping of the Kv blocker XE991 and pilocarpine. Importantly, XE991 is Kv 7selective and inhibition of Kv 7 is known to be excitatory, functionally mimicking the direct agonism of muscarinic receptors, which in turn initiates channel closure and aspects of cholinergic stimulation (Brown, 2018).
The other class of compounds showing notable similarities amongst their functional connectivity phenotypes were the monoaminergic reuptake inhibitors. There was also subdivision of these compounds that may be reflective of the subtle balance in efficacy against each (zebrafish) transporter (Huot et al., 2015) and further work with more selective modulators would help to clarify this. Analysis of the brain regions showing the greatest changes in functional connectivity supported effects on noradrenaline and 5-HT circuitry. Aminergic neurons are also amongst the first to develop in zebrafish. Noradrenergic pathways mainly originate from the locus coeruleus; dopaminergic projections from the posterior tuberculum (Tay et al., 2011) and 5-HT neurons from the raphé nuclei (Maximino & Herculano, 2010).
Linking these regions is the important relay station, the interpeduncular nucleus (Bianco & Wilson, 2009), which was also prominently connected here. Altered functional connectivity with the tegmentum is also notably associated with compounds showing relative dopaminergic activity (e.g. cocaine, amphetamine and nomifensine). Given that the tegmentum is known to be relatively rich in zebrafish D 2A receptors (Maximino & Herculano, 2010), this may suggest that dopaminergic activity is influencing these phenotypes in contrast with other groupings where the more subtle balance between NET and SERT activity is the predominant driving force. The co-location of fentanyl, meperidine and phencyclidine with several monoaminergic reuptake inhibitors could be explained by interactions between these neurotransmitter systems, for example, the reported inhibition of the human SERT and NET by meperidine at low μM concentrations (Rickli et al., 2018) and the inhibition of NA and DA reuptake at low μM concentrations in rat hippocampus (Snell et al., 1988).
Using our analysis method, some pharmacological classes showed few similarities in their functional imaging phenotypes despite being extremely neuroactive in 4dpf larvae. Amongst the GABA A modulators tested, for example, this dissimilarity may be driven by differences in receptor subtype distribution and relative ligand efficacy (Bandara et al., 2020;Johnston, 2013). Alternatively, this may be reflective of the broad distribution of GABA A receptors across the CNS, meaning that targets with a more punctate spatial distribution, such as with the monoaminergics, are more easily differentiated using functional connectivity analysis. Similarly, despite being neuroactive neither NMDA nor (RS)-(tetrazol-5-yl)glycine exhibited any particularly strong functional connections. This could be the result of agonism-induced global depolarization resulting in blockade of synaptic transmission in which the Ca 2+ rise (and associated GCamp fluorescence) would remain due to Ca 2+ entry via the NMDAR itself and depolarization-mediated gating of voltage-gated Ca 2+ channels.
Extending previous findings using PTZ (Diaz Verdugo et al., 2019; Liu & Baraban, 2019), here we demonstrate altered functional connectivity resulting from compound treatment across a wide range of pharmacological mechanisms of action implicated in seizure generation, particularly amongst compounds used as precipitants in animal models of seizure and epilepsy. The presence of abnormal synchronisation is a key feature of seizure networks, as reported in experimental and clinical assessments of epileptic seizures (Jiruska et al., 2013).
Our data also revealed an increased prevalence of altered functional connectivity involving nodes within the mesencephalon, in particular, as well as in the diencephalon after exposure to compounds associated with seizures. This echoes the results of previous authors using broader scale segmentation (e.g. Diaz Verdugo et al., 2019;Rosch et al., 2018), who have suggested large midbrain structures, such as the optic tectum, play an important role as hubs in the transition of preictal to ictal network dynamics after exposure to the GABA A antagonist PTZ. Whole brain connectivity analysis also revealed a reduced overall functional connectivity after exposure to several compounds, most notably pilocarpine and bicuculline. Cymerblit-Sabba and Schiller (2012) found that overall network synchronicity decreased in the initial stages of a seizure after systemically applied pilocarpine or an intrahippocampal applied mixture of pilocarpine, kainate and picrotoxin in rats and that hypersynchronisation developed only as seizures progressed and intensified. It may, therefore, be the case here that longer term exposure of some compounds is required in order to reach the advanced stages of a seizure and allow a hypersynchronous phenotype to emerge.
Concluding, we present data demonstrating the power of functional brain imaging in larval zebrafish for assessing the action of neuroactive drugs in a highly relevant vertebrate model. These data help us to understand the relevance of the 4dpf larval zebrafish for neuropharmacological studies and reveal that even at this early development stage, these animals are highly responsive to a wide range of neuroactive compounds across multiple primary mechanisms of action. The 4dpf larva appears a particularly relevant model for the study of drugs affecting the monoaminergic and cholinergic systems and has potential utility for the assessment of new drugs targeting these circuits. Furthermore, the assessment of neural activity levels and overall network functional connectivity could help to identify unwanted seizure liability in drugs as part of a new drug CNS safety assessment cascade, as well as aiding in the identification and characterisation of new anti-epileptic drugs and therapeutic targets (Ghannad-Rezaie et al., 2019). This approach also has high potential as an in vivo model for the study of seizurogenic and epileptogenic mechanisms more widely, especially with advances in spatiotemporal resolution driving a move towards real time-imaging of large-scale neuronal assemblies at cellular, or even subcellular, resolutions.
the Replacement, Refinement and Reduction of Animals in Research

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. Some data may not be made available because of privacy or ethical restrictions.