Reconstruction scheme for excitatory and inhibitory dynamics with quenched disorder: application to zebrafish imaging

An inverse procedure is developed and tested to recover functional and structural information from global signals of brains activity. The method assumes a leaky-integrate and fire model with excitatory and inhibitory neurons, coupled via a directed network. Neurons are endowed with a heterogenous current value, which sets their associated dynamical regime. By making use of a heterogenous mean-field approximation, the method seeks to reconstructing from global activity patterns the distribution of in-coming degrees, for both excitatory and inhibitory neurons, as well as the distribution of the assigned currents. The proposed inverse scheme is first validated against synthetic data. Then, time-lapse acquisitions of a zebrafish larva recorded with a two-photon light sheet microscope are used as an input to the reconstruction algorithm. A power law distribution of the in-coming connectivity of the excitatory neurons is found. Local degree distributions are also computed by segmenting the whole brain in sub-regions traced from annotated atlas.


I. INTRODUCTION
The unmatched ability of the brain to cope with an extraordinarily large plethora of complex tasks, carried out in parallel, ultimately resides in the intricate web of interlinked connections which define the architecture of the embedding neurons network. Structural and functional information, as inferred from direct measurements of neuronal activity, under different experimental conditions, are fundamental pieces of a jigsaw puzzle of how the brains works, from simple organisms to more complicated creatures, across phylogenetic scales. Suitable methods have been developed which build on statistical mechanics tools, e.g. maximum entropy principles [1,2], to resolve the functional map that orchestrates the coordinated firing of neurons dislocated in different portions of the brain. However, other sources of heterogeneity in the brain should be accounted for as well. Neuronal excitability, namely the ability of neurons to respond to external inputs, is finely controlled through inhibitory/excitatory balance [3,4]. Furthermore, individual neurons can display a variable degree of inherent excitability, a source of spatial quenched disorder which reflects back in the ensuing activation patterns.
Motivated by this, in [5] we proposed and tested against both synthetic and real data, an inverse scheme to quantify the statistics of neurons' excitability, while inferring, from global activity measurements, the, a priori unknown, distribution of network connectivities. The * Corresponding author; chicchi.lorenzo@gmail.com method employs an extended model of Leaky-Integrate and Fire (LIF) neurons, with short-term plasticity. Only excitatory neurons are accounted for in [5]. These are assumed to be coupled via a directed network and display a degree of heterogeneity in the associated current, which sets the firing regime in which a neuron operates. The inverse scheme builds on the celebrated Heterogenous Mean-Field (HMF) approximation [6][7][8][9] and seeks to recover the distribution of the (in-degree) connectivity of (excitatory) neurons, concurrently with the distribution of the assigned currents, denoted by a. The HMF approximation was previously employed in [10][11][12][13] to reconstruct the topology of an underlying network from artificially generated data, meant to mimic neuronal signals. In [5] the approach was generalized so as to account for the dynamical heterogeneity, as stemming from the intrinsic degree of individual neurons excitability. As mentioned above, individual excitability acts as a key component of the dynamics and yields irregular patterns of activity like those displayed in real measurements. The reconstruction scheme was applied in [5] to longitudinal wide-field fluorescence microscopy data of cortical functionality in groups of awake mice and enabled us to identify altered distributions in neuron excitability immediately after the stroke, and in agreement with earlier observation [14][15][16]. Conversely, rehabilitation allowed to recover a distribution similar to pre-stroke conditions.
The goal of this work is to push forward the reconstruction algorithm by accounting for the simultaneous presence of both excitatory and inhibitory neurons, in a refined variant of the inversion scheme proposed in [5]. Notice that already in [12] intertangled families of exci-tatory and inhibitory neurons have been considered, but only in a simplified setting where currents were assumed to be homogeneous. Relaxing this ansatz proves however mandatory when aiming at bridging the gap between theory and experiments, a challenge that we shall hereafter tackle. In particular, we will recast the dynamics of the examined LIF model in a reduced setting by grouping in different classes (excitatory and inhibitory) neurons which bear distinct values of the current a and of the connectivity k.
Our extended inverse method aims at computing the distribution of the currents, as well as the distributions of the connectivities, for respectively excitatory and inhibitory neurons, via an iterative scheme which selfconsistently identifies the classes of neurons needed to interpolate the global activity field, supplied as an input to the algorithm. First, the performance of the method is evaluated in silico, against synthetically generated data. We then move forward to considering a direct application of the developed technique to custom-made two-photon Light Sheet (2PLS) microscope, optimized for high-speed (1 Hz) volumetric imaging of zebrafish larva (ZL, Danio rerio). Near infrared (NIR) light is used for excitation, covering a wavelength range that is not visible to the larva in order not to induce unwanted visual responses. Hence 2PLS microscopy allows to record whole-brain activity with high temporal and spatial resolution, by preventing undesired external bias [17][18][19]. The experimental input is processed with a properly devised methodology to return a spatially resolved raster plot for the spiking activity of neurons over time, which provides the ideal input for the reconstruction method to work. A power law distribution of the in-coming connectivity of excitatory neurons is found, which is robust over a significant range of the imposed fraction of inhibitory neurons. Local degree distributions are also recovered by partitioning the whole brain in bound sub-domains, traced from annotated atlas [20,21]. When manipulating experimental data one cannot distinguish among the contributions resulting from different neurons (excitatory vs. inhibitory). A procedure is however developed which allows for the degree distribution of the excitatory neurons to be determined, while accounting for the role exerted by the population of inhibitory ones.
The paper is organized as follows. Next section is devoted to introducing the reconstruction scheme and to challenge its performance against synthetically produced data. We will then turn to presenting the experimental platform and discuss the details that relate to data processing. Processed data are supplied as an input to the mathematical reconstruction scheme to yield the results which are presented in Section III B. Finally we will sum up and draw our conclusions.

II. GENERAL MATHEMATICAL FRAMEWORK
A. The model We use a Leaky Integrate and Fire (LIF) model to mimic the dynamics of individual neurons. Consider a pool of N LIF neurons and denote by v i (t) the membrane potential of neuron i. Further, label with I syn i (t) the synaptic current due to the incoming connections with other neurons of the collection.
The dynamics of the membrane potential is hence ruled by the following equation where a i stands for the external input current of neuron i. This is a crucial quantity, as it sets the dynamics of the corresponding neuron, in the uncoupled regime. The critical value a i = 1 separates between quiescent and active (spiking) regimes. When v i reaches the threshold value v th , neuron i emits a spike and the membrane potential v i is reset to the base value v r . Following [12], the membrane potential is rescaled by a suitable amount to have the spike threshold set at v th = 1 and the rest potential at v r = 0. The Tsodyks, Uziel, and Markram model [22,23] is assumed to describe the interactions among neurons, i.e. their coupling dynamics. More specifically, the dynamics of a synapse is expressed in terms of fractions of three different neurotransmitter states: active (y ij (t)), available (x ij (t)) and inactive (z ij (t)), where i and j stand respectively for post-synaptic and pre-synaptic neurons. The obvious constraint x ij (t) + y ij (t) + z ij (t) = 1 applies, for any time t.
During a spike of the pre-synaptic neuron, a fraction u ij of neurotransmitters in the available state is activated. The time evolution of the three stateṡ takes as input the spike train S j (t) = m δ(t − t * j (m)) of pre-synaptic neuron j emitting its n-th spike at time t j (n) [24]. The time is rescaled to the membrane time constant τ m = 30ms. The time constants are set to τ in = 0.2, τ i r = 3.4 if the post-synaptic neuron i is inhibitory, or τ i r = 26.6 if it is excitatory [12]. If the post-synaptic neuron i is excitatory, the fraction u ij (t) of neurotransmitters is set to the constant value of U = 0.5. Otherwise, u ij (t) evolves in timė where τ f = 33.25 and U f = 0.08 [12,24].
Equations (2) are coupled to Eq. (1) via the synaptic current where g is the coupling parameter and A ij stands for the elements of the network adjacency matrix A. The matrix entry is A ij = 1, if a link exists which goes from j to i, and provided j is an excitatory neuron. Conversely, A ij = −1 if the starting node j identifies inhibitory neuron. On the other hand, if A ij = 0 nodes i and j are not directly connected. Equations (2) can be cast in a more compact form so as to favour insight into the inspected processes and reduce the associated computational costs. To achieve this, we first notice that Eqs. (2) only depend on the pre-synaptic neuron j and the characteristics of the postsynaptic neuron i. Equation (2a) can hence be split into two distinct equations: where the apexes E (Excitatory) and I (Inhibitory) reflect the specificity of the target neuron. Similar arguments apply to Eqs. (2b), (2c), and (3). For each type of synapse, the average field, i.e. the fraction of neurotransmitters in the active state, is calculated as where E and I are the ensemble of excitatory and inhibitory neurons, respectively; f I stands for the fraction of inhibitory neurons in the network. The global fields are defined as B. The heterogeneous mean field ansatz.
The model described in the previous section is reformulated here in terms of a Heterogeneous Mean Field (HMF) approximation. The original neurons are classified according to their characteristics. More specifically, neurons of the same type (excitatory or inhibitory) and with the same incoming connectivity k and external current a are considered identical. Therefore, L × M equivalence classes are defined, where L (M ) is the number of sub-intervals in which the value range of k (a) has been divided. Moreover, we assume that neurons in the class k are subjected to a synaptic current proportional to their in-degree, i.e., for the excitatory and inhibitory neurons, respectively. Following this assumption, the model can be rewritten asv where ( †, * ) identify all possible pairs of post-synaptic and pre-synaptic neurons. Denote by P E (k) and P I (k) the in-degree distributions for the excitatory and inhibitory neurons, respectively, and P (a) the external current distribution. Equations (6) can be closed by the consistency relations Taking in account the discretization of the defined classes of equivalence, Eq. (10) turns intõ where we have implicitly introduced the discrete counterpart of the continuous probability distributions P * (k) = (P * (k 1 ), P * (k 2 ), ..., P * (k L )) and P(a) = (P (a 1 ), P (a 2 ), .., P (a M )).

C. Reconstruction scheme
In this section, we set up a general reconstruction scheme for recovering the a priori unknown distributions for the in-degree P E (k) and P I (k), as well as the external current P (a). This is achieved by interpolating the available global fields Y E (t) and Y I (t), under the simplified HMF descriptive framework. The main steps of the reconstruction algorithm are schematically depicted in Fig.  1A.
As already mentioned, we assume the global fields Y E (t) and Y I (t) to be given. Then we integrate Eqs. (6), by using these global fields Y E (t) and Y I (t) as inputs to the model. The equations are initialized with vari- k,a (t 0 )) randomly drawn from a uniform distribution, and generated so as to respect the constraints v † k,a (t 0 ) < 1 and y k,a (t 0 ) < 1. Forced by the external fields Y E (t) and Y I (t), the governing equations are integrated forward and the variables y † * k,a (t) are stored for each class (k, a), type of synapse ( †, * ), and time t. This process is repeated for H independent realizations of the initial conditions. The average fraction of neurotransmitters in the active state, for each class (k, a) and synapse type, is computed at any time of Then, the approximated global fieldsỸ † * are calculated, via Eq. (11), for an initial guess of the distributions P E (k), P I (k), and P(a). These latter are then recursively modified so as to improve the correspondence between the approximated fields and their true homologues Y † * .
Formally, for each ( †, * ) ∈ {EE, EI, IE, II}, we aim at minimizing the function (12) Note that the arguments of the above function are the target probability distributions P * (k) and P(a) that one aims at inferring.
The iterative algorithm operates as follows. The distribution P(a) is initially frozen to a given profile and the . . .
where the first row reflects the normalization condition. The problem is hence reduced to a linear system that can be readily solved to obtain a first estimate of the distributions P E (k) and P I (k).
As second step, the in-degree distributions P E (k) and P I (k) are fixed to the solutions found at the previous iteration. Similar to step one, we evaluate the quantities y † * am (t) = L l=1 P * (k l ) y † * k l ,am (t) and formulate the linear problem The above problem can be solved to obtained an updated estimate for the P(a). The overall procedure, consisting of two nested steps, is iterated until a maximum number of allowed cycles is reached, or, alternatively, the stopping criterion is eventually met.
Before proceeding in the analysis, we introduce a slightly modified notation. The in-degree k is normalized to the size of the network N . In formulae we will set k = k/N , withk ∈ [0, 1]. From hereon the distributions that constitute the target of the reconstruction scheme will be hence expressed as function ofk, instead of k.

A. Synthetic data
In this section we test the proposed reconstruction protocol against synthetic data. The reconstruction scheme was successfully validated on synthetic data for the case of homogeneous external current in [5,12]. Here, we test the reconstruction method on synthetic networks of excitatory and inhibitory neurons, assuming a quenched distribution of heterogeneous external currents. To this end we generate a random graph with N nodes whose structural characteristic is contained in the signed N ×N adjacency matrix A, which specifies the existence of (directed) links among pairs of adjacent nodes. Following the convention introduced above, negative entries (A ij = −1) indicate that the starting node (j) is of the inhibitory type, whereas for positive elements (A ij = 1) j belongs to the family of excitatory neurons. The network generation procedure is conceived so as to return a bell shaped distribution for both P E (k) and P I (k) (see Fig. 1B). Quenched disorder in the input currents is introduced, the assigned currents being distributed according to a uni-modal profile P (a) (see Fig. 1B). These are the exact distributions that we eventually seek to recover via the aforementioned reconstruction algorithm. Note that the domain of definition of P (a) includes the bifurcation value a = 1.
With this setting, Eqs. (1) and (2) are integrated and the fields Y E and Y I are calculated by using Eqs. (6) and (7). This is possible because we have access to all information which pertain to the network architecture and to the heterogenous collection of randomly generated currents.
The recorded global fields Y E and Y I are used as inputs to the reconstruction scheme presented in the previous section. Figure 1B shows the comparison between true and estimated distributions, at the end of the reconstruction procedure and for one generation of the synthetic network. By inserting the estimated distributions into Eq. (11), we obtain the global fieldsỸ E andỸ I . The comparison between estimated (Ỹ E ,Ỹ I ) and true (Y E , Y I ) global fields, as obtained by working on the index space, is presented in Fig. 1C. The agreement is excellent for both the inhibitory and the excitatory components.
When working with experimental data, however, one cannot isolate the contributions stemming from different neurons, grouped according to their specific traits (excitatory vs. inhibitory). This implies that the sums in Eq. (6), i.e. the input to the envisaged reconstruction scheme, cannot be in general accessed, as it was instead the case when working in the framework of the synthetic model considered above. To overcome this intrinsic limitation, we propose and test an alternative route, which performs adequately well when challenged against synthetic data. The idea is to propose an approximated version of the input Eqs. (6). To this end, we extend the sums which run on the excitatory neurons to all neurons and compute, under this approximation, the fields Y EE and Y IE . To validate this hypothesis we operate with synthetic networks with unclassified neurons. It can be shown that the approximation for the fields Y EE , Y IE correctly describes Y EE , but not Y IE . This conclusion is supported by systematic numerical investigations (data not shown), that we carried out by varying the relative proportion of excitatory and inhibitory neurons. Building on the above, we therefore write: which enables us to compute the approximated field Y EI as Finally, we can estimate Y E as: Remark that the above expression for Y E is obtained without grouping the neurons in excitatory and inhibitory classes, but provided f I , the fraction of inhibitory neurons, is eventually known. As we lack information on the corresponding field Y I , we can run the reconstruction scheme in a simplified setting which is solely targeted to reconstructing the probability distributions P (a) and P E (k). In Fig. 1D, the results of the revisited inversion method are displayed, having set f I to the correct value, i.e. assuming the relative proportion of excitatory and inhibitory neurons which has been effectively employed in generating the synthetic dataset. The reconstruction algorithm is still capable of returning a faithful representation of both P (a) and P E (k).
Conversely, when f I is set to zero, the reconstruction scheme compensates for the missing inhibitory component by predicting a reduced average connectivity of the excitatory population, as compared to the correct value assumed in the data generation scheme, see Fig. 1E. In the following, we will apply the reconstruction scheme in this latter version to the analysis of 2PLS microscope images of living zebrafish larva.

B. Experimental data
In this section we apply our reconstruction framework to calcium fluorescence microscopy data of zebrafish larva brain. Indeed, each time a neuron fires an action potential, the strong depolarization occurring triggers a rapid and transient increase of intracellular calcium concentration [25]. Thus, following calcium-dependent fluorescence dynamics represents an indirect measurement of neuronal spiking activity [26]. A description of the experimental set up is provided in the Methods section.

Data processing
In order to apply the inverse scheme to real data, it is necessary to pre-process the wide field calcium images with the purpose of first identifying the location of the neurons. From individual traces of each spotted neuron, we will single out the spiking events, record the times of occurrence and build up the corresponding raster plot. This will serve as input to the reconstruction algorithm.
To this aim, data are first downsampled to 2 × 2 pixels so that the new pixel size is comparable to the neurons' nuclear size ( Fig. 2A). For every new pixel, the maximum value of the calcium fluorescence is calculated and only the pixels with maximum intensity projection (MIP) above a threshold are identified as neurons. The threshold is fixed to the average MIP value over all the image pixels. Furthermore, we operate a moving average to remove low frequency fluctuation (Fig. 2B) and select as presumed neurons the pixels that show large asymmetry in the recorded traces. To implement this step, we compute the skewness from individual time series of the calcium activity and classify as neurons the pixels that yield a sufficiently skewed signal 1 . Figure 2C outlines the different phases of the process for one of the layers of the collection. More specifically, in Fig. 2C.i the MIP of the down-sampled pixels are depicted, in grey scale. In Fig. 2C.ii a binary representation of the whole brain is displayed where only pixels with MIP above the fixed threshold are highlighted. Lastly, in Fig. 2C.iii only the pixels which exhibit a strong asymmetric signal, i.e. the neurons with skewness above the imposed threshold, are shown. At the end of the selection process the number of identified neurons is around (1 ÷ 3) × 10 3 for each layer, which correspond to a total of 49 × 10 3 .
Once the neurons are identified, we proceed to construct the raster plot. To this aim, for each selected neuron we analyze the time series of the calcium fluorescence to remove the background noise and detect events, which we call spikes. More specifically, a spike is defined as the time of threshold crossing. The thresholds are set to the mean value of the recorded time series plus two times the associated standard deviation. In addition, in order to avoid double detections due to noise, we discarded all events that succeeded the previous event by less than a minimum inter-event interval of 5 data points (5 seconds) 2 [28]. The general overview of the spike trains emitted by neurons in a sample layer results in a raster plot (Fig. 2D). Time is on the horizontal axis, whereas the vertical axis displays the neuron indices. Each spike of neuron i is associated with a red dot in the row i, at the corresponding time of spiking.

Results
As described in the previous section, we process 3D calcium fluorescence data so as to identify pixels containing neurons. Figure 3 shows the results of this identification for eight different layers of the zebrafish brain. Colors reflect the average cross-correlation 3 at lag zero of 1 The skewness threshold is here put to 0.4. In doing so we select a number of putative neurons which is comparable to the known size of a zebrafish larva brain [27] 2 An additional analysis has been carried out using a minimum inter-event interval of 2 data points. From the corresponding raster plot, the global fields have been calculated and they appear indistinguishable from the ones obtained using a minimum interevent interval of 5 data points. 3 Cross-correlation measures the similarity between two series at different time shifts, or lags. In formulae, the cross-correlation between two vectors xt and yt at lag τ is defined as the expected each neuron with all other selected neurons of the brain. The higher the correlation, the more reddish the color displayed. The patterns of correlations are rather symmetric, an observation which can be interpreted as an a posteriori validation of the implemented procedure for automatic neurons selection. A movie which allows to navigate across successive layers of the whole 3D stack can be found in the SM. The processing of data explained above allows us to obtain an experimental raster plot describing the events, or spikes, associated to each neuron. Indeed, the raster plot contains information about the spike train function S i (t), for all neurons i, and can be readily employed to recover the global field Y E to be supplied as an input to the reconstruction scheme. More explicitly, the experimentally determined S i (t) is used to integrate Eqs. (2) and (3), by breaking the coupling with Eq. (1) which sets the evolution of the membrane potential 4 . This is of great advantage since Eq. (1) contains the specific information about the network connections, i.e., the adjacency matrix elements, which are a priori unknown. In other words, the raster plot provides us a way to compute the global field (input of the reconstruction process described in section II C) without knowing the underlying structure of the network.
Since the true fraction f I of inhibitory neurons is unknown, the global field Y E in the approximated form [Eq. (17)] is computed for different values of f I . In particular, we first reconstruct the in-degree distribution P E (k) for the excitatory neurons and the external current distribution P (a) for different values of f I and we store the results 5 . Secondly, we computed the reconstruction error M SE = 1/T , where the asterisk denotes the complex conjugation. 4 The reactions parameters are set to the nominal values as declared above [12]. The same parameters are assumed in the reconstruction, and, in this respect, the model equations acts as a filter to transform the supplied fluorescence data in the ideal input for the inverse procedure to run. 5 The coupling strength g is set to the (experimentally justified) value of 30 [12,24,29] adopted in the forward simulations of the model.
interval, the reconstructed distributions show common features. In particular, the external current distribution P (a) is peaked in the vicinity of the critical value a = 1 (Fig. 4C). The neurons associated to a > 1 get self-excited and promote the activation of other neurons which would be instead quiescent in the uncoupled limit. The small bumps that are found for relatively large values of the intensities a can be traced back to the high frequency component of the signal to be interpolated, and, as such, bear limited fundamental interests. The largescale dynamics of the recorded time-series, including the heterogeneous modulation of the macroscopic field oscillations, is instead encoded in the distribution of intensities that define the bulk of P (a), i.e. the limited excerpt of curve which is found in correspondence of the leftmost portion of the support in a.
The reconstructed in-degree distributions P E (k) for the excitatory neurons, at different values of f I , are depicted in Fig. 4E, in log-log scale. Although over a limited support ink, the obtained distributions seem to display a power law decay, P E (k) ∝k −α , the characteristic exponent α being only modestly influenced by the chosen value of f I . Our findings suggest that excitatory neurons are organized in a network with few hubs and many more peripheral nodes.
As an additional test, we partition the full set of identified nodes in 10 different populations, reflecting distinct anatomical regions, as follows available atlases [21,27]. The reconstruction algorithm is then applied to each of the selected region, treated as independent from the surrounding context, so as to access the local degree distribution. Results, displayed in Fig. 5, are in line with those reported in Fig. 4. Moreover, neurons characterized by a significant connectivity, the above referenced hubs, seem unevenly distributed across different anatomical regions.

IV. CONCLUSIONS
Reconstructing structural and functional information from brain activity represents a topic of outstanding importance, which can in principle trigger applied and fundamental fallout. In [5] we presented, and successfully tested, an inverse scheme which aimed at inferring the distributions of both firing rates and networks connectivity, from global activity fields. The method builds on the Leaky-Integrate and Fire (LIF) model which we modified by the inclusion of quenched disorder, in the assigned individual currents. The imposed degree of heterogeneity in the currents yields non trivial a-periodic patterns, which resemble those recorded in vivo. The dynamical model considered in [5] solely included the population of excitatory neurons. Starting from these background, we here have generalized the reconstruction procedure of [5] so as to account for the simultaneous presence of both excitatory and inhibitory neurons, while still dealing with the effect of the current heterogeneity. The dynamics of the examined multi-species LIF model is recast in a simplified framework, by grouping together neurons that belong to the same class (inhibitory vs. excitatory), while sharing the similar currents and in-degree. The output of the reduced model, driven by the excitatory and inhibitory global fields, is self-consistently used to seed an iterative scheme which seeks at fitting the supplied fields, via suitably adjusting the unknown distributions. These latter are the distributions of the incoming degrees for, respectively, excitatory and inhibitory neurons, as well as the distribution of the imposed currents. The method is tested on synthetic data and yields satisfying performances. Having in mind applications to real data, we also dealt with a setting where it is not possible to separate the contribution that pertains to the excitatory component from that stemming from the inhibitory counterpart. In this case, we propose and test a procedure which enables to recover the distribution of incoming degrees for the network of excitatory neurons (assumed predominant), while gauging the role exerted by inhibitors.
The devised protocol is then applied to whole-brain functional data resulting from light-sheet calcium imaging of a zebrafish larva. The experimental input is processed with an automatic procedure which allows us to identify putative neurons , and to extract their fluorescence signal. Remarkably, the cross-correlation maps produced show a high grade of clusterization, which faithfully matches the anatomical boundaries of multiple brain regions identified using zebrafish brain atlases [20,21]. From the calcium signal displayed by each selected neuron we build up an experimental representation of the raster plot of the spiking activity of the zebrafish brain, which forms the input to the reconstruction scheme. A power law distribution of the in-coming connectivity of excitatory neurons is found, which is only modestly affected by the imposed fraction of inhibitors. Local degree distributions are also reconstructed by analysing the signal from specific regions, which correspond to distinct anatomical areas. Interestingly, the anatomical districts considered in the analysis can be divided into two different groups, according to the reconstructed probability distributions of both their excitatory incoming connections P E (k) and excitability P (a). The first group, including dorsal thalamus, medial tegmentum, superior raphe, hindbrain and spinal cord, is characterized by higher connections and lower excitability. Conversely, the second group, comprising telencephalon, habenulae, optic tectum and cerebellum, is described by lower connections and higher excitability. The reconstruction scheme reflects the specific functional connectivity of the larval brain during spontaneous activity precisely under these experimental conditions. Indeed, during measurements the zebrafish larva is embedded in agarose, and thus exposed to a diffused tactile stimulation, which could explain the higher incoming connections calculated for the dorsal thalamus, a sensory relay station [30,31]. Furthermore, despite mechanically and pharmacologically immobilized, larval at-tempts to escape the restrained condition could account for the higher incoming excitatory connections calculated for dorsal raphe (whose activity has been correlated with arousal state, vigilance and responsiveness [32]) and for the most caudal regions, namely hindbrain and spinal cord, responsible for the initiation of motor behaviours [33,34]. Moreover, in this scenario we observe a lower probability distribution of incoming excitatory connections for cerebellum. This result may be associated to the function of motor coordination and refinement [35,36] of this region, typically related to the actual execution of a movement. Finally, since the measurement is performed in complete darkness, the larva is not exposed to any visual cue (the IR laser used for two-photon excitation is not perceived by the larval visual system) and this could explain the lower incoming connections calculated for the optic tectum, the main retinorecipient structure in the zebrafish brain [37,38]. The large-scale oscillations in the recorded time series reflect back in the recovered distribution of currents: a significant fraction of neurons appear to operate in the quiescent state, while a minority self-excite to orchestrate the dynamics of the ensemble. The existence of possible correlations between individual connectivities and associated neurons currents cannot be resolved within the proposed approach and defines an interesting target for future investigations.

A. Validity of the HMF approximation
We here challenge the predictive ability of the HMF approximation. To this end, we first calculate the average inter-spike interval (ISI) -the average distance in time between successive spikes -for the system in its original formulation, i.e. in the space of the nodes. The computed ISI is confronted to the homologous quantity obtained under the HMF scenario. The comparison is drawn in Fig. 6, for both excitatory and inhibitory neurons, and confirm the accuracy of the reduced HMF scheme.

B. Experimental setup
The experimental optical setup employed is a modified version of the setup described in detail in [39]. Briefly, 930 nm NIR light is generated by a pulsed titaniumsapphire oscillator (Chameleon Ultra II, Coherent) and conditioned by a pulse compressor (PreComp, Coherent). After being attenuated by a combination of a halfwave plate and a GlanThompson polarizer, the beam passes through an electro-optical modulator that periodically switches its polarization plane orientation be-tween two states: parallel or orthogonal to the optical table surface. Then a pair of retarders are used to pre-compensate for polarization distortions. After that, the beam is routed to the galvanometric mirror assembly and scanned by a resonant galvanometric mirror (CRS-8 kHz, Cambridge Technology) along larval rostro-caudal direction, to generate the digitally-scanned LS, and by a closed-loop galvanometric mirror (6215H, Cambridge Technology) along larval dorso-ventral direction. Finally, the beam is relayed to an objective (XLFLUOR4X/340/0,28,Olympus), placed at the lateral side of the larva, by a scan-lens (50 mm focal length), a tube-lens (75 mm focal length) and a pair of relay lenses (250 mm and 200 mm focal lengths).
Differently from the setup described in [39], a polarizing beam splitter (PBS) is present between the tube-lens and the first relay lens and the optical components downstream the PBS are duplicated on the opposite side of the larva. In this way the excitation light is steered on one optical arm or on the other depending on its instantaneous polarization orientation by the PBS.
The detection arm of the microscope is identical to what described in [39]. A water-immersion objective (XLUMPLFLN20XW, Olympus) placed on the dorsal side of the larva collects the emitted green fluorescent light while being scanned along the axial dimension by an objective scanner (PIFOC P-725.4CD,Physik Instrumente) synchronously with the closed-loop galvanometric mirror oscillations. The optical image is then spectrally filtered (FF01-510/84-25 nm BrightLine,Semrock) and projected on a camera (ORCA-Flash4.0 V3, Hamamatsu) sensor by a sequence of two tube lenses (300 mm and 200 mm focal lengths) and an objective (UP-LFLN10X2, Olympus).
The larvae were imaged with volumetric acquisitions (frequency: 1 Hz) composed by 31 planes spaced by 5 µm and with a pixel size of about 2 × 2 µm 2 .
We employed a transgenic strain of zebrafish larvae Tg(elavl3:H2B-GCaMP6s) [40,41] in homozygous albino background that express the fluorescent calcium indicator GCaMP6s under a pan-neuronal promoter and with a nuclear localization. Sample mounting was performed as described in [42]. Briefly, before the acquisition each larva was immersed in a solution of the paralyzing agent d-tubocurarine (2 mM; 93750, Sigma-Aldrich), included in 1.5% (w/v) low gelling temperature agarose (A9414, Sigma-Aldrich) in fish water (150 mg/L Instant Ocean, 6.9 mg/L NaH2PO4, 12.5 mg/L Na2HPO4, pH 7.2) and mounted on a custom-made glass support immersed in thermostated fish water. The animals were maintained according to standard procedures [43] and observed at 4 days post fertilization. Fish maintenance and handling were carried out in accordance with European and Italian law on animal experimentation (D.L. 4 March 2014, no. 26).  1. A) Schematic outline of the reconstruction procedure. The global fields Y E (t) and Y I (t) constitute the inputs of the model in the HMF approximation (i). Different choices for the probability distributions P E (k), P I (k), and P (a) are iteratively tested in order to find the best match between the input fields and the reconstructed fields, as obtained by using the equations displayed in the red box (ii). B) Outcome of the reconstruction procedure: the true probability distributions of a synthetic network are compared with those obtained with the proposed reconstruction method. A random network with N = 5000 is considered here. The fraction of inhibitory neurons is set to f I = 0.05. The number of classes defined in the HMF approximation for the in-degree and the external current is L = 250 and M = 250 respectively. C) Comparison between the true global fields and the ones obtained via the reconstructed distributions. The plot in the inset is a zoom in of a peak. D-E) Outputs of the reconstruction are compared with the true external current probability distribution P (a) and the true in-degree distribution P E (k) for the excitatory neurons of the same network; the network is made of N = 1000 neurons of which a fraction f I = 0.2 are inhibitors. In the HMF approximation one hundred classes have been defined for both the in-degree and external current, namely, L = 100 and M = 100. In D) the correct fraction of inhibitory neurons is taken into account, while in E) the inhibitory neuron effects are not considered.