Emergence of time persistence in a data-driven neural network model

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

513 views

Version of Record published: May 10, 2023 (This version)
Accepted Manuscript published: March 14, 2023 (Go to version)
Accepted: March 13, 2023
Received: April 17, 2022
Preprint posted: February 4, 2022 (Go to version)

1. Of interest
Rodent ultrasonic vocal interaction resolved with millimeter precision using hybrid beamforming

Max L Sterling, Ruben Teunisse, Bernhard Englitz

Tools and Resources Sep 26, 2023
Further reading

Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Appendix 1
Appendix 2
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

Establishing accurate as well as interpretable models of network activity is an open challenge in systems neuroscience. Here, we infer an energy-based model of the anterior rhombencephalic turning region (ARTR), a circuit that controls zebrafish swimming statistics, using functional recordings of the spontaneous activity of hundreds of neurons. Although our model is trained to reproduce the low-order statistics of the network activity at short time scales, its simulated dynamics quantitatively captures the slowly alternating activity of the ARTR. It further reproduces the modulation of this persistent dynamics by the water temperature and visual stimulation. Mathematical analysis of the model unveils a low-dimensional landscape-based representation of the ARTR activity, where the slow network dynamics reflects Arrhenius-like barriers crossings between metastable states. Our work thus shows how data-driven models built from large neural populations recordings can be reduced to low-dimensional functional models in order to reveal the fundamental mechanisms controlling the collective neuronal dynamics.

Editor's evaluation

The authors show how high-dimensional neural signals can be reduced to low-dimensional models with variables that can be directly linked to behavior. The reduced model can account for long time scales of persistent activity that arise from transitions between metastable model states. The authors further show that the rate of these transitions is modulated by water temperature according to the classic Arrhenius law, although the results for different temperatures could not yet be unified into a single description based on real external temperature.

https://doi.org/10.7554/eLife.79541.sa0

Introduction

How computational capacities emerge from the collective neural dynamics within large circuits is a prominent question in neuroscience. Modeling efforts have long been based on top-down approaches, in which mathematical models are designed to replicate basic functions. Although they might be very fruitful from a conceptual viewpoint, these models are unable to accurately reproduce actual data and thus remain speculative. Recently, progress in large-scale recording and simulation techniques has led to the development of bottom-up approaches. Machine-learning models, trained on recorded activity, allow for the decoding or the prediction of neuronal activity and behavior (Glaser et al., 2020; Pandarinath et al., 2018). Unfortunately, the blackbox nature of these data-driven models often obscures their biological interpretation, for example, the identification of the relevant computational units (Butts, 2019). This calls for quantitative, yet interpretable approaches to illuminate the functions carried out by large neural populations and their neuronal substrate.

This work is an attempt to do so in the specific context of the anterior rhombencephalic turning region (ARTR), a circuit in the zebrafish larva that drives the saccadic dynamics and orchestrates the chaining of leftward/rightward swim bouts (Ahrens et al., 2013; Dunn et al., 2016; Wolf et al., 2017; Ramirez and Aksay, 2021; Leyden et al., 2021). The ARTR spontaneous activity exhibits temporal persistence, that is, the maintenance of activity patterns over long (∼ 10 s) time scales. This functional feature is ubiquitous in the vertebrate brain. It plays an essential role in motor control, as best exemplified by the velocity position neural integrator, a circuit that integrates neural inputs and allows for a maintenance of the eye position after an ocular saccade (Seung, 1996; Seung et al., 2000; Miri et al., 2011). Temporal persistence is also central to action selection (Wang, 2008) and short-term memory storage (Zaksas and Pasternak, 2006; Guo et al., 2017). As isolated neurons generally display short relaxation times, neural persistence is thought to be an emergent property of recurrent circuit architectures (Zylberberg and Strowbridge, 2017). Since the 1970s, numerous mechanistic network models have been proposed that display persistent activity. They are designed such as to possess attractor states, that is, stable activity patterns toward which the network spontaneously converges.

Although attractor models are conceptually appealing, assessing their relevance in biological circuits remains challenging. To this aim, recent advances in machine learning combined with large-scale methods of neural recordings may offer a promising avenue. We hereafter focus on energy-based network models, trained to replicate low-order data statistics, such as the mean activities and pairwise correlations, through the maximum entropy principle (Jaynes, 1957). In neuroscience, such models have been successfully used to explain correlation structures in many areas, including the retina (Schneidman et al., 2006; Cocco et al., 2009; Tkačik et al., 2015), the cortex (Tavoni et al., 2016; Tavoni et al., 2017; Nghiem et al., 2018), and the hippocampus (Meshulam et al., 2017; Posani et al., 2017) of vertebrates, and the nervous system of Caenorhabditis elegans (Chen et al., 2019). These models are generative, that is, they can be used to produce synthetic activity on short time scales, but whether they can reproduce long-time dynamical features of the biological networks remains an open question.

Here, we first report on spontaneous activity recordings of the ARTR network using light-sheet functional imaging at various yet ethologically relevant temperatures. These data demonstrate that the water temperature controls the persistence time scale of the ARTR network, and that this modulation is in quantitative agreement with the thermal dependence of the swimming statistics. We then infer energy-based models from the calcium activity recordings and show how these data-driven models not only capture the characteristics and probabilities of occurrence of activity patterns, but also reproduce the observed thermal dependence of the persistent time scale. We further derive a mathematically tractable version of our energy-based model, called mean-field approximation, whose resolution provides a physical interpretation of the energy landscape, of the dynamical paths there in, and of their changes with temperature. We finally extend the model to incorporate visual stimulation and correctly reproduce the previously reported visually driven ARTR dynamics (Wolf et al., 2017). This work establishes the capacity of data-driven network inference to numerically emulate persistent dynamics and to unveil fundamental network features controlling such dynamics.

Results

The water temperature controls behavioral and neuronal persistence time scales in zebrafish larvae

In this first section, we report on functional recordings of the ARTR dynamics performed at various temperatures (18–33°C). We show that the persistent time scale that characterizes the ARTR’s endogenous dynamics is thermally modulated. This dependence is reflected in the change in swimming statistics observed in freely swimming assays. We further characterize how the water temperature impacts the distribution of activity patterns.

ARTR endogeneous dynamics is thermally modulated

We used light-sheet functional imaging to record the ARTR activity in zebrafish larvae expressing a calcium reporter pan-neuronally (Tg(elavl3:GCaMP6)). The larvae, embedded in agarose, were placed in a water tank whose temperature was controlled in the range 18–33°C (see Appendix 2—figure 1A). ARTR neurons were identified using a combination of morphological and functional criteria, as detailed in Wolf et al., 2017. Their spatial organization is displayed in Figure 1A, for all recorded animals after morphological registration on a unique reference brain (145 ± 65 left neurons, 165 ± 69 right neurons, mean ± SD across 13 different fish, see Appendix 2—table 1). For each neuron, an approximate spike train $s (t)$ was inferred from the fluorescence signal using Bayesian deconvolution (Tubiana et al., 2020). A typical raster plot of the ARTR is shown in Appendix 2—figure 1B (recorded at 26°C), together with the mean signals of the left and right subcircuits, $m_{L, R} (t) = \frac{1}{N_{L, R}} \sum_{i \in L, R} s_{i} (t)$ .

Figure 1

Download asset Open asset

Temperature dependence of anterior rhombencephalic turning region (ARTR) dynamics and turn direction persistence.

(A) Morphological organization of the ARTR showing all identified neurons from 13 fish recorded with light-sheet calcium imaging. (B) Example of ARTR binarized signal sign ( $m_{L} - m_{R}$ ) (gray) along with the left ( $m_{L}$ , red) and right ( $m_{R}$ , blue) mean activities. (C) Averaged power spectra of the ARTR binarized signals for the five tested temperatures. The dotted vertical lines indicate the signal switching frequencies $ν$ as extracted from the Lorentzian fit (solid lines). (D) Temperature dependence of $ν$ . The lines join data points obtained with the same larva. (E) Swimming patterns in zebrafish larvae. Swim bouts are categorized into forward and turn bouts based on the amplitude of the heading reorientation. Example trajectory: each dot corresponds to a swim bout; the color encodes the reorientation angle. (F) The bouts are discretized as left/forward/right bouts. The continuous binary signal represents the putative orientational state governing the chaining of the turn bouts. (G) Power spectra of the discretized orientational signal averaged over all animals for each temperature (dots). Each spectrum is fitted by a Lorentzian function (solid lines) from which we extract the switching rate $k_{f l i p}$ . (H) Temperature dependence of $k_{f l i p}$ . Inset: relationship between $k_{f l i p}$ (behavioral) and $ν$ (neuronal) switching frequencies. Bar sizes represent SEM, and the dashed line is the linear fit.

To analyze the thermal dependence of the ARTR dynamics, we extracted from these recordings a binarized ARTR signal, $sign (m_{L} (t) - m_{R} (t))$ (see Figure 1B and Appendix 2—figure 1C for example signals from the same fish at different temperatures). The average power spectra of these signals for the five tested temperatures (average of 3–8 animals for each temperature, see Appendix 2—table 1) are shown in Figure 1C. We used a Lorentzian fit to further extract the alternation frequency $ν$ for each dataset (Figure 1C, solid lines). This frequency was found to increase with the temperature (Figure 1D). Although $ν$ could significantly vary across specimen at a given temperature, for a given animal, increasing the temperature induced an increase in the frequency in 87.5% of our recordings (28 out of 32 pairs of recordings).

In this analysis, we used the binarized ARTR activity to facilitate the comparison between behavioral and neural data, as described in the next section. However, the observed temperature dependence of the left-right alternation time scale was preserved when the spectra were computed from the ARTR activity, $m_{L} (t) - m_{R} (t)$ (see Appendix 2—figure 1D).

Impact of the water temperature on the turn direction persistence in freely swimming larvae

It has previously been shown that the ARTR governs the selection of swim bout orientations: turn bouts are preferentially executed in the direction of the most active (right or left) ARTR subcircuit (Dunn et al., 2016; Wolf et al., 2017), such that $sign (m_{L} (t) - m_{R} (t))$ constitutes a robust predictor of the turning direction of the animal; see Figure 5 - figure supplement 2E in Dunn et al., 2016. Therefore, the temporal persistence of the ARTR dynamics is reflected in a turn direction persistence in the animal’s swimming pattern, that is, the preferred chaining of similarly orientated turn bouts.

We thus sought to examine whether the thermal dependence of the ARTR endogenous dynamics could manifest itself in the animal navigational statistics. In order to do so, we used the results of a recent study (Le Goc et al., 2021), in which 5–7-day-old zebrafish larvae were video-monitored as they swam freely at constant and uniform temperature in the same thermal range (Figure 1E). We quantified the time scale of the turn direction persistence by assigning a discrete value to each turn bout: −1 for a right turn, +1 for a left turn (forward scouts were ignored). We then computed an orientational state signal continuously defined by the value of the last turn bout (Figure 1F). The power spectra of the resulting binary signals are shown in Figure 1G for various temperatures. We used a Lorentzian fit (‘Materials and methods,’ Equation 6) to extract, for each experiment, a frequency $k_{f l i p}$ . This rate, which defines the probability of switching orientation per unit of time, systematically increases with the temperature, from 0.1 to 0.6 s⁻¹ (Figure 1H). Increasing the temperature thus leads to a progressive reduction of the turn direction persistence time. The inset plot in Figure 1H establishes that the left/right alternation rates extracted from behavioral and neuronal recordings are consistent across the entire temperature range (slope = 0.81, $R = 0.99$ ).

ARTR activity maps are modulated by the temperature

We then investigated how the water temperature impacts the statistics of the ARTR activity defined by the mean activity of the left and right sub-populations, $m_{L}$ and $m_{R}$ . The probability maps in the $(m_{L}, m_{R})$ plane are shown in Figure 2A for two different temperatures, with the corresponding raster plots and time signals of the two subcircuits. At high temperature, the ARTR activity map is confined within an L-shaped region around $(m_{L} = 0, m_{R} = 0)$ and the circuit remains inactive for a large fraction of the time. Conversely, at lower temperature, the ARTR activity is characterized by long periods during which both circuits are active and shorter periods of inactivity. We quantified this thermal dependence of the activity distribution by computing the log-probability of the activity of either region of the ARTR at various temperatures (Appendix 2—figure 2A). The occupation rate of the inactive state ( $m_{L, R} \sim 0$ ) increases with temperature, with a corresponding steeper decay of the probability distribution of the activity. Consistently, we found that the mean activities $m_{L}$ and $m_{R}$ decreased with temperature (Appendix 2—figure 2B). Such a dependence might reflect varying levels of temporal coherence in the activity of the ARTR with the temperature. In order to test this, we computed the Pearson correlation at various temperature but we saw no clear dependency of the average correlation across ipsilateral or contralateral pairs of neurons (Appendix 2—figure 2C).

Figure 2

Download asset Open asset

Ising models reproduce characteristic features of the recorded activity.

(A) (Top) Probability densities $P (m_{L}, m_{R})$ , see Equation 2, of the activity state of the circuit (obtained from the spiking inference of the calcium data), in logarithmic scale, and for two different fish and water temperatures T = 20 and T = 30°C; Color encodes z-axis (same color bar for both). (Middle) 10-min-long raster plots of the activities of the left (red) and right (blue) subregions of the anterior rhombencephalic turning region (ARTR). (Bottom) Corresponding time traces of the mean activities $m_{L}$ and $m_{R}$ . (B) Processing pipeline for the inference of the Ising model. We first extract from the recorded fluorescence signals approximate spike trains using a Bayesian deconvolution algorithm (BSD). The activity of each neuron is then ‘0’ or ‘1.’ We then compute the mean activity and the pairwise covariance of the data, from which we infer the parameters $h_{i}$ and $J_{i j}$ of the Ising model. Finally, we can generate raster plot of activity using Monte Carlo sampling. (C) Same as (A) for the two corresponding inferred Ising models. The raster plots correspond to Monte Carlo-sampled activity, showing slow alternance between periods of high activity in the L/R regions. Here we show only two examples of a qualitative experimental vs. synthetic signals comparison. We provide in the supplementary materials the same comparison for every recording.

Our analysis thus indicates that the water temperature modulates both the endogenous dynamics and the activity distribution of the ARTR. For both aspects, we noticed a large variability between animals at a given temperature. This is not unexpected as it parallels the intra- and inter-individual variability in the fish exploratory kinematics reported in Le Goc et al., 2021. Nevertheless, we observed a strong positive correlation between the persistence time and the mean activity across animals and trials for a given temperature (Appendix 2—figure 2D and ‘Materials and methods’), indicating that both features of the ARTR may have a common drive.

A data-driven energy-based model reproduces the statistics of the ARTR dynamics

Our aim was to reproduce the ARTR spontaneous activity using an energy-based data-driven network model. The inference pipeline, going from raw fluorescence data to the model, is summarized in Figure 2B. We first reconstructed an estimated spike train for each ARTR neuron using a deconvolution algorithm (Tubiana et al., 2020). We divided the recording window ( $T_{r e c} \sim 1200 s$ for each session) in time bins whose width was set by the imaging frame rate ( $d t = 100 - 300 m s$ ). Each dataset thus consisted of a series of snapshots $s^{k} = (s_{1}^{k}, \dots, s_{N}^{k})$ of the ARTR activity at times $k$ , with $k = 1, \dots, T_{r e c} / d t$ ; here, $s_{i}^{k} = 1$ if cell $i$ is active or $s_{i}^{k} = 0$ if it is silent in time bin $k$ .

We then computed the mean activities, ${⟨ s_{i} ⟩}_{data}$ , and the pairwise correlations, ${⟨ s_{i} s_{j} ⟩}_{data}$ , as the averages of, respectively, $s_{i}^{k}$ and $s_{i}^{k} s_{j}^{k}$ over all time bins $k$ . We next inferred the least constrained model, according to the maximum entropy principle (Jaynes, 1957), that reproduced these quantities. This model, known as the Ising model in statistical mechanics (Ma, 1985) and probabilistic graphical model in statistical inference (Koller and Friedmann, 2009), describes the probability distribution over all $2^{N}$ possible activity configurations $s$ ,

P (s) = \frac{1}{Z} \exp (\sum_{i} h_{i} s_{i} + \sum_{i < j} J_{i j} s_{i} s_{j}),

where $Z$ is a normalization constant. The bias $h_{i}$ controls the intrinsic activity of neuron $i$ , while the coupling parameters $J_{i j}$ account for the effect of the other neurons $j$ activity on neuron $i$ (‘Materials and methods’). The set of parameters ${h_{i}, J_{i j}}$ were inferred using the Adaptative Cluster Expansion and the Boltzmann machine algorithms (Cocco and Monasson, 2011; Barton and Cocco, 2013; Barton et al., 2016). Notice that in Equation 1, the energy term in the parenthesis is not scaled by a thermal energy as in the Maxwell–Boltzmann statistics. We thus implicitly fix the model temperature to unity; of course, this model temperature has no relation with the water temperature $T$ . Although the model was trained to reproduce the mean activities and pairwise correlations (see Appendix 2—figure 3A–C and ‘Materials and methods’ for fourfold cross-validation), it further captured higher-order statistical properties of the activity such as the probability that $K$ cells are active in a time bin (Appendix 2—figure 3D; Schneidman et al., 2006).

Once inferred, the Ising model can be used to generate synthetic activity configurations $s$ . Here, we used a Monte Carlo (MC) algorithm to sample the probability distribution $P (s)$ in Equation 1. The algorithm starts from a random configuration of activity, then picks up uniformly at random a neuron index, say, $i$ . The activity $s_{i}$ of neuron $i$ is then stochastically updated to 0 or to 1, with probabilities that depend on the current states $s_{j}$ of the other neurons (see Equation 8 in ‘Materials and methods’ and code provided). The sampling procedure is iterated, ensuring convergence toward the distribution $P$ in Equation 1. This in silico MC dynamics is not supposed to reproduce any realistic neural dynamics, except for the locality in the activity configuration $s$ space.

Figure 2C shows the synthetic activity maps and temporal traces of Ising models trained on the two same datasets as in Figure 2A. For these synthetic signals, we use MC rounds, that is, the number of MC steps divided by the total number of neurons (‘Materials and methods’), as a proxy for time. Remarkably, although the Ising model is trained to reproduce the low-order statistics of the neuronal activity within a time bin only, the generated signals exhibit the main characteristics of the ARTR dynamics, that is, a slow alternation between the left and right subpopulations associated with long persistence times; see raster plots in Figure 2C.

Comparison of experimental and synthetic ARTR dynamics across recordings

We repeated the inference procedure described above for all our 32 recordings (carried out with $n = 13$ fish and 5 different water temperatures, see Appendix 2—table 2) and obtained the same number of sets of biases and couplings. We first compared the distributions of the left-right mean activity $m_{L} = \frac{1}{N_{L}} \sum_{i \in L} s_{i}$ and $m_{R} = \frac{1}{N_{R}} \sum_{i \in R} s_{i}$ extracted from the data and from the Ising model. In order to do so, we used the Kullback–Leibler (KL) divergence, a classical metrics of the dissimilarity between two probability distributions. The distribution of the KL divergences between the experimental test datasets (see ‘Materials and methods’) and their associated Ising models is shown in green in Figure 3A. The KL values were found to be much smaller than those obtained between experimental test datasets and Ising models trained from different recordings (red distribution). This result establishes that the Ising model quantitatively reproduces the ARTR activity distribution associated to each specimen and temperature.

Figure 3

Download asset Open asset

Comparison of model distributions and persistence times across fish and water temperatures.

(A) Distribution of the Kullback–Leibler divergences between test datasets and their corresponding Ising models (green), between test datasets and Ising models trained on different datasets (red) and between test datasets and their corresponding independent models that assume no connections between neurons (dark blue). Note that each dataset is divided in a training set corresponding to 75% of the time bins chosen randomly and a test set comprising the remaining 25%. (B) Average persistence times in simulations vs. experiments. Each dot refers to one fish at one water temperature; colors encode temperature.

This agreement crucially relies on the presence of inter-neuronal couplings in order to reproduce the pairwise correlations in the activity: a model with no connection (i.e., the independent model, see ‘Materials and methods’) fitted to reproduce the neural firing rates offers a very poor description of the data (see Figure 3A [dark blue distribution] and Appendix 2—figure 3E–G).

Finally, we examined to what extent the synthetic data could capture the neural persistence characteristics of the ARTR. The persistence times extracted from the data and from the MC simulations of the inferred models were found to be strongly correlated (Figure 3B, $R = 0.84$ ). The MC dynamics thus captures the inter-individual variability and temperature dependence of the ARTR persistent dynamics.

Spatial organization and temperature dependence of the Ising inferred parameters

In all recordings, inferred ipsilateral couplings are found to be centered around a positive value (std = 0.12, mean = 0.062), while contralateral couplings are distributed around 0 (mean = –0.001, std = 0.10); see Appendix 2—figure 4A–C. Still, a significant fraction of these contralateral couplings are strongly negative. We illustrated this point by computing the fraction of neuronal pairs ( $i, j$ ) that are contralateral for each value of the coupling $J_{i j}$ or the Pearson correlation (Appendix 2—figure 4D and E). Large negative values of couplings or correlations systematically correspond to contralateral pairs of neurons, whereas large positive values correspond to ipsilateral pairs of neurons.

In addition, we found that the ipsilateral couplings $J_{i j}$ decay, on average, exponentially with the distance between neurons $i$ and $j$ (Appendix 2—figure 4F), in agreement with findings in other neural systems (Posani et al., 2018). Spatial structure is also present in contralateral couplings (Appendix 2—figure 4G). Biases display a wide distribution ranging from –8 to 0 (std = 1.1, mean = –4.1, Appendix 2—figure 5A–C), with no apparent spatial structure.

We next examined the dependency of the Ising model parameters on the water temperature. To do so, for each fish, we selected two different water temperatures, and the corresponding sets of inferred biases and couplings, ${h_{i}, J_{i j}}$ . We then computed the Pearson correlation coefficient $R^{2}$ of the biases and of the coupling matrices at these two temperatures (inset of Appendix 2—figure 6). We saw no clear correlation between the model parameters at different temperatures, as shown by the distribution of $R^{2}$ computed across fish and across every temperatures (Appendix 2—figure 6).

Mean-field study of the inferred model unveils the energy landscape underlying the ARTR dynamics

Mean-field approximation to the data-driven graphical model

While our data-driven Ising model reproduces the dependence of the persistence time scale and activity distribution on the water temperature, why it does so remains unclear. To understand what features of the coupling and local bias parameters govern these network functional properties, we turn to mean-field theory. This powerful and mathematically tractable approximation scheme is commonly used in statistical physics to study systems with many strongly interacting components (Ma, 1985). In the present case, it amounts to deriving self-consistent equations for the mean activities $m_{L}$ and $m_{R}$ of the left and right ARTR subpopulations (Figure 4A and Appendix 1).

Figure 4

Download asset Open asset

Mean-field approximation of the inferred Ising model.

(A) Schematic view of the mean-field Ising model. (B) Examples of simulated $m_{L}$ and $m_{R}$ signals of the mean-field dynamical equations for two sets of parameters that correspond to fish ID #5 at two water temperatures (22°C and 30°C), see Appendix 2—table 2. (C) Free-energy landscapes in the ( $m_{L}$ , $m_{R}$ ) plane computed with the mean-field model. These data correspond to the same sets of parameters as in panel (B). Colored circles denote metastable states, and the line of black arrows indicates the optimal path between $(m^{l o w}, m^{l o w})$ and $(m^{l o w}$ , $m^{h i g h})$ states. (D) Schematic view of the free energy along the $m_{R}$ axes. The arrows denote the energy barriers $Δ$ $ℱ$ associated with the various transitions. The dark green arrow denotes $Δ F ((m^{h i g h}, m^{l o w}) \to (m^{l o w}, m^{l o w}))$ ; the purple arrow denotes $Δ F ((m^{l o w}, m^{l o w}) \to (m^{h i g h}, m^{l o w}))$ . (E) Values of the free-energy barriers as a function of temperature. Error bars are standard error of the mean (32 recordings, n = 13 fish at 5 different water temperatures). (F) Persistence time of the mean-field anterior rhombencephalic turning region (ARTR) model for all fish and runs at different experimental temperatures. Each dot refers to one fish at one temperature; colors encode temperature.

Within mean-field theory, each neuron $i$ is subject to (i) a local bias $H$ , (ii) an excitatory coupling $J > 0$ from the neurons in the ipsilateral region, and (iii) a weak coupling $I$ from the neurons in the contralateral side. These three parameters were set as the mean values of, respectively, the inferred biases $h_{j}$ and the inferred ipsilateral and contralateral interactions $J_{i j}$ . In addition, we introduce an effective size $K$ of each region to take into account the fact that mean-field theory overestimates interactions by replacing them with their mean value. This effective number of neurons was chosen, in practice, to best match the results of the mean-field approach to the full Ising model predictions (see Appendix 1, Appendix 2—table 2 and Appendix 2—figure 7A–C). It was substantially smaller than the number $N$ of recorded neurons. The selection method used to delineate the ARTR populations may yield different number of neurons in the $L$ and $R$ regions (see Appendix 2—table 1). This asymmetry was accounted for by allowing the parameters $H$ , $J$ , and $K$ defined above to take different values for the left and right sides.

Mean-field theory thus allowed us to reduce the data-driven Ising model, whose definition requires $\frac{1}{2} (N_{L} + N_{R}) (N_{L} + N_{R} + 1)$ parameters ${h_{i}, J_{i j}}$ , to a model depending on seven parameters ( $H_{L}, H_{R}, J_{L}, J_{R}, K_{L}, K_{R}, I$ ) only (Figure 4A), whose values vary with the animal and the experimental conditions, for example, temperature (Appendix 2—table 2).

Free energy and Langevin dynamics

The main outcome of the analytical treatment of the model is the derivation of the so-called free energy $ℱ (m_{L}, m_{R})$ as a function of the average activities $m_{L}$ and $m_{R}$ ; see Appendix 1. The free energy is a fundamental quantity as it controls the density of probability to observe an activation pattern $(m_{L}, m_{R})$ through

P (m_{L}, m_{R}) \propto e^{- F (m_{L}, m_{R})}

Consequently, the lower the free energy $ℱ$ , the higher the probability of the corresponding state $(m_{L}, m_{R})$ . In particular, the minima of the free energy define persistent states of activity in which the network can be transiently trapped.

The free energy landscape can be used to simulate dynamical trajectories in the activity space $(m_{L}, m_{R})$ . To do so, we consider a Langevin dynamics in which the two activities $m_{L} (t), m_{R} (t)$ evolve in time according to the stochastic differential equations,

τ \frac{d m_{L}}{d t} (t) = - \frac{\partial F}{\partial m_{L}} (m_{L} (t), m_{R} (t)) + ϵ_{L} (t),

τ \frac{d m_{R}}{d t} (t) = - \frac{\partial F}{\partial m_{R}} (m_{L} (t), m_{R} (t)) + ϵ_{R} (t),

where $τ$ is a microscopic time scale, and $ϵ_{L} (t), ϵ_{R} (t)$ are white noise ‘forces’, $⟨ ϵ_{L} (t) ⟩ = ⟨ ϵ_{R} (t) ⟩ = 0$ , independent and delta-correlated in time: $⟨ ϵ_{L} (t) ϵ_{R} (t^{'}) ⟩ = 0$ , $⟨ ϵ_{L} (t) ϵ_{L} (t^{'}) ⟩ = ⟨ ϵ_{R} (t) ϵ_{R} (t^{'}) ⟩ = 2 δ (t - t^{'})$ . This Langevin dynamical process ensures that all activity configurations $(m_{L}, m_{R})$ will be sampled in the course of time, with the expected probability as given by Equation 2.

Figure 4B shows the mean-field simulated dynamics of the left and right activities, $m_{L}$ and $m_{R}$ , with the parameters corresponding to two Ising models at two different temperatures in Figure 2C. We observe, at low temperatures, transient periods of self-sustained activity (denoted by $m^{h i g h}$ ) of one subcircuit, while the other has low activity ( $m^{l o w}$ ) (see time trace 1 in Figure 4B). At high temperature, high activity in either (left or right) area can be reached only transiently (see trace 2 in Figure 4B). These time traces are qualitatively similar to the ones obtained with the full inferred Ising model and in the data (Figure 2A and C, bottom).

Barriers in the free-energy landscape and dynamical paths between states

We show in Figure 4C the free-energy landscape in the $(m_{L}, m_{R})$ plane for the same two conditions as in Figure 4B. The minimization conditions $\frac{\partial ℱ}{\partial m_{L}} = \frac{\partial ℱ}{\partial m_{R}} = 0$ provide two implicit equations over the activities $m_{L}^{*}, m_{R}^{*}$ corresponding to the preferred states. For most datasets, we found four local minima: the low-activity minimum $(m_{L}^{*}, m_{R}^{*}) = (m^{l o w}, m^{l o w})$ , two asymmetric minima, $(m^{h i g h}, m^{l o w})$ and $(m^{l o w}, m^{h i g h})$ , in which only one subregion is strongly active, and a state in which both regions are active, $(m^{h i g h}, m^{h i g h})$ . The low-activity minimum $(m^{l o w}, m^{l o w})$ is the state of lowest free energy, hence with largest probability, while the high-activity state $(m^{h i g h}, m^{h i g h})$ has a much higher free energy and much lower probability. The free energies of the asymmetric minima $(m^{h i g h}, m^{l o w})$ and $(m^{l o w}, m^{h i g h})$ lie in between, and their values strongly vary with the temperature.

The Langevin dynamics defines the most likely paths (see ‘Materials and methods’) in the activity plane joining one preferred state to another, for example, from $(m^{h i g h}, m^{l o w})$ to $(m^{l o w}, m^{h i g h})$ as shown in Figure 4C. Along these optimal paths, the free energy $ℱ$ reaches local maxima, defining barriers to be overcome in order for the network to dynamically switchover (purple and green arrows in Figure 4C). The theory of activated processes stipulates that the average time to cross a barrier depends exponentially on its height $Δ ℱ$ :

t (Δ ℱ) \sim τ \times e^{Δ ℱ},

up to proportionality factors of the order of unity (Langer, 1969). Thus, the barrier $Δ F ((m^{h i g h}, m^{l o w}) \to (m^{l o w}, m^{l o w}))$ shown in dark green in Figure 4D controls the time needed for the ARTR to escape the state in which the left region is active while the right region is mostly silent, and to reach the all-low state. The barrier $Δ ℱ ((m^{l o w}, m^{l o w}) \to (m^{h i g h}, m^{l o w}))$ shown in purple is related to the rising time from the low-low activity state to the state where the right region is active, and the left one is silent.

Within mean-field theory, we estimated the dependence in temperature of these barriers height (Figure 4E and Appendix 2—figure 7D) and of the associated persistence times (Figure 4F). While substantial variations from animal to animal were observed, we found that barriers for escaping the all-low state and switching to either $L, R$ region increase with the water temperature. As a consequence, at high temperature, only the low–low activity state is accessible in practice to the system, and the mean activity remains low (see Appendix 2—figure 2D), with fluctuations within the low–low state. Conversely, at low water temperatures, barriers separating the low–low and the active high–low or low–high states are weaker, so the latter become accessible. As a first consequence, the mean activity is higher at low temperature (Appendix 2—figure 2D). Furthermore, the system remains trapped for some time in such an active state before switching to the other side, for example, from high–low to low–high. This is the origin of the longer persistence time observed at low temperature.

Ising and mean-field models with modified biases capture the ARTR visually driven dynamics

While the analyses above focused on the spontaneous dynamics of the ARTR, our data-driven approach is also capable of explaining activity changes induced by external and time-varying inputs. In order to illustrate this capacity, we decided to reanalyze a series of experiments, reported in Wolf et al., 2017, in which we alternatively illuminated the left and right eyes of the larva, for periods of 15–30 s, while monitoring the activity of the ARTR (Figure 5A) with a two-photon light-sheet microscope. During and after each stimulation protocol, 855 s of spontaneous activity was recorded on $n = 6$ fish. We found that the ARTR activity could be driven by this alternating unilateral visual stimulation: the right side of the ARTR tended to activate when the right eye was stimulated and vice versa (Figure 5B).

Figure 5

Download asset Open asset

Modified Ising model captures the behavior of anterior rhombencephalic turning region (ARTR) under visual stimulation.

(A) Scheme of the stimulation protocol. The left and right eyes are stimulated alternatively for periods of 15–30 s, after which a period of spontaneous (no stimulus) activity is acquired. (B) Example of the ARTR activity signals under alternated left–right visual stimulation. The small arrows indicate the direction of the stimulus. (C) Sketch of the modified Ising model, with additional biases $δ h_{i}$ to account for the local visual inputs. (D) Values of the additional biases averaged over the ipsilateral and contralateral (with respect to the illuminated eye) neural populations. n = 6 fish. (E) Monte Carlo activity traces generated with the modified Ising model. (F) Free-energy landscapes computed with the mean-field theory during spontaneous (left panel) and stimulated (right panel) activity for an example fish. (G) Free-energy along the optimal path as a function of $m_{L} - m_{R}$ during spontaneous (plain line) and stimulated (dotted line) activity. The model is the same as in panel (F).

To analyze these datasets, we first followed the approach described in Figure 2B, and inferred, for each fish, the sets of biases $h_{i}$ and interactions $J_{i j}$ using the spontaneous activity recording only (Appendix 2—table 3). In a Appendix 2—table 2second step, we exploited recordings of the visually driven activity to infer additional biases $δ h_{i}$ to the neurons, while keeping the interactions $J_{i j}$ fixed (Figure 5C); in practice, we defined two sets of additional biases, $δ {\overset{\leftarrow}{h}}_{i}$ and $δ {\vec{h}}_{i}$ , corresponding, respectively, to leftward and rightward illuminations. The underlying intuition is that biases encode inputs due to the stimulation, while the interactions between neurons can be considered as fixed over the experimental time scale. This simplified model reproduces the low order statistics of the data under stimulation (Appendix 2—figure 8A and B).

The inferred values of the additional biases, averaged over the entire subpopulation (right or left), are shown in Figure 5D for both ipsiversive or contraversive stimulation. The results show that light stimulation produces a strong increase of excitability for the ipsilateral neurons and a smaller one for contralateral neurons.

We then simulated the visual stimulation protocol by sampling the Ising model while alternating the model parameters, from ${h_{i} + δ {\vec{h}}_{i}, J_{i j}}$ to ${h_{i} + δ {\overset{\leftarrow}{h}}_{i}}, J_{i j}}$ , and back. The simulated dynamics of the model (Figure 5E) qualitatively reproduces the experimental traces of the ARTR activity (Figure 5B). In particular, the model captures the stabilizing effect of unilateral visual stimuli, which results in a large activation of the ipsilateral population, which in turn silences the contralateral subcircuit due to the negative $I$ coupling between both. This yields the anticorrelation between the left and right sides clearly visible in both the experimental and simulated traces, and much stronger in the case of spontaneous activity (Appendix 2—figure 8C–F).

To better understand the Ising dynamics under visual stimulation, we resort, as previously, to mean-field theory. For asymmetric stimulation, our mean-field model includes, during the periods of stimulation, extra biases $Δ H_{L}$ and $Δ H_{R}$ over neurons in, respectively, the left and right areas (Figure 5C), while the couplings $J$ and $I$ remain unchanged. We show in Figure 5F the free-energy $ℱ$ as a function of $m_{L}, m_{R}$ for an example fish. Due to the presence of the extra bias, the landscape is tilted with respect to its no-stimulation counterpart (Figure 5G), entailing that the left- or right-active states are much more likely, and the barrier separating them from the low–low state is much lower. As a consequence, the time necessary for reaching the high-activity state is considerably reduced with respect to the no-stimulation case (see Equation 5). These results agree with the large probability of the high-activity states and the fast rise to reach these states in the Ising traces in Figure 5E compared with Figure 2C.

Discussion

Modeling high-dimensional data, such as extensive neural recordings, imposes a trade-off between accuracy and interpretability. Although highly sophisticated machine-learning methods may offer quantitative and detailed predictions, they might in turn prove inadequate to elucidate fundamental neurobiological mechanisms. Here, we introduced a data-driven network model, whose biologically grounded architecture and relative simplicity make it both quantitatively accurate and amenable to detailed mathematical analysis. We implemented this approach on functional recordings performed at various temperature of a key population of neurons in the zebrafish larvae brain, called ARTR, that drives the orientation of tail bouts and gaze (Dunn et al., 2016; Wolf et al., 2017; Ramirez and Aksay, 2021; Leyden et al., 2021).

First, we demonstrate that the persistent time scale of the ARTR endogenous dynamics decreases with the temperature, mirroring the thermal modulation of turn direction persistence in freely swimming behavioral assays. We then demonstrate that our energy-based model not only captures the statistics of the different activity patterns, but also numerically reproduces the endogenous pseudo-oscillatory network dynamics, and their thermal dependence. The inferred Ising model is then analyzed within the so-called mean-field formulation, in which the coupling and bias parameters are replaced by their values averaged over the left and right subpopulations. It yields a two-dimensional representation of the network energy landscape where the preferred states and associated activation barriers can be easily evaluated. We show how this combined data-driven and theoretical approach can be applied to analyze the ARTR response to transient visual stimulation. The latter tilts the energy landscape, strongly favoring some states over others.

Origin and functional significance of the temperature dependence of the ARTR dynamics

The brains of cold-blooded animals need to operate within the range of temperature that they experience in their natural habitat, for example, 18–33°C for zebrafish (Gau et al., 2013). This is a peculiarly stringent requirement since most biophysical processes are dependent on the temperature. In some rare instances, regulation mechanisms might stabilize the circuit dynamics in order to preserve its function, as best exemplified by the pyloric rhythm of the crab whose characteristic phase relationship is maintained over an extended temperature range (Tang et al., 2010). Yet in general, an increase in temperature tends to increase the frequency of oscillatory processes (Robertson and Money, 2012). The observed acceleration of the ARTR left/right alternation with increasing temperature could thus directly result from temperature-dependent cellular mechanisms. Furthermore, one cannot rule out the possibility that the ARTR dynamics could also be indirectly modulated by temperature via thermal-dependent descending neuromodulatory inputs.

As a result of this thermal modulation of the neuronal dynamics, many cold-blooded animals also exhibit temperature dependence of their behavior (Long and Fee, 2008; Neumeister et al., 2000; Stevenson and Josephson, 1990). Here, we were able to quantitatively relate the two processes (neuronal and motor) by demonstrating that an increase in temperature consistently alters the pattern of spontaneous navigation by increasing the left/right alternation frequency. Interpreting the functional relevance of this modification of the swimming pattern is tricky since many other features of the animal’s navigation are concurrently impacted by a change in temperature, such as the bout frequency, turning rate, turn amplitude, etc. Nevertheless, we were able to show in a recent study that this thermal dependence of the swimming kinematic endows the larva with basic thermophobic capacity, thus efficiently protecting them from exposure to the hottest regions of their environment (Le Goc et al., 2021).

Ising model is not trained to reproduce short-term temporal correlations, but is able to predict long-term dynamics

The graphical model we introduced in this work was trained to capture the low-order statistics of snapshots of activity. Because graphical models are blind to the dynamical nature of the population activity, it is generally believed that they cannot reproduce any dynamical feature. Nevertheless, here we demonstrate that our model can quantitatively replicate aspects of the network long-term dynamics such as the slow alternation between the two preferred states. To better understand this apparent paradox, it is necessary to distinguish short and long time scales. At short time scale, defined here as the duration of a time bin (of the order of a few 100 ms), the model cannot capture any meaningful dynamics. The MC algorithm we used to generate activity is an abstract and arbitrary process, and the correlations it produces between successive time bins cannot reproduce the ones in the recording data. Capturing the short-term dynamics would require a biologically grounded model of the cell–cell interactions, or, at the very least, to introduce parameters capturing the experimental temporal correlations over this short time scale (Marre et al., 2009; Mézard and Sakellariou, 2011).

Yet, the inability of the Ising model to reproduce short time dynamical correlations does not hinder its capacity to predict long-time behavior. The separation between individual neuronal processes (taking place over time scales smaller than 100 ms) and network-scale activity modulation, which happens on time scales ranging from 1 to 20 s, is here essential. The weak dependence of macroscopic processes on microscopic details is in fact well known in many fields outside neuroscience. A classic example is provided by chemical reactions, whose kinetics are often controlled by a slow step due to the formation of the activated complex and to the crossing of the associated energy barrier $Δ E$ , requiring a time proportional to $e^{Δ E / (k T)}$ . All fast processes, whose modeling can be very complex, contribute an effective microscopic time scale $τ$ in Arrhenius’ expression for the reaction time (see Equation 5). In this respect, what really matters to predict long time dynamical properties is a good estimate of $Δ E$ or, equivalently, of the effective energy landscape felt by the system. This is precisely what the Ising model is capable of doing. This explains why, even if temporal information are not explicitly included in the training process, our model may still be endowed with a predictive power over the long-term network dynamics.

Energy-landscape-based mechanism for persistence

In a preceding article (Wolf et al., 2017), we developed a mathematical model of the ARTR in which the left and right ARTR population were represented by a single unit. To account for the ARTR persistent dynamics, an intrinsic adaptation time scale had to be introduced in an ad hoc fashion. While the mean-field version of the inferred Ising model shows some formal mathematical similarity with this two-unit model, it differs in a fundamental aspect. Here, the slow dynamics reflects the itinerant exploration of a two-dimensional energy landscape (Figure 4C), for which the barriers separating metastable states scale linearly with the system size. The time to cross these barriers in turn grows exponentially with the system size, as prescribed by Arrhenius law, and can be orders of magnitude larger than any single-neuron relaxation time. Persistence is therefore an emerging property of the neural network.

Mean-field approximation and beyond

The mean-field approach, through a drastic simplification of the Ising model, allows us to unveil the fundamental network features controlling its coarse-grained dynamics. Within this approximation, the distributions of couplings and of biases are replaced by their average values. The heterogeneities characterizing the Ising model parameters (Appendix 2—figure 4 and Appendix 2—figure 5), and ignored in the mean-field approach, may, however, play an important role in the network dynamics.

In the Ising model, the ipsilateral couplings are found to be broadly distributed such as to possess both negative and positive values. This leads to the presence of so-called frustrated loops, that is, chains of neurons along which the product of the pairwise couplings is negative. The states of activities of the neurons along such loops cannot be set in a way that satisfies all the excitatory and inhibitory connections, hence giving rise to dynamical instabilities in the states of the neurons. The absence of frustrated loops in the network (Figure 4A) stabilizes and boosts the activity, an artifact we had to correct for in our analytical treatment by introducing an effective number of neurons $K$ , much smaller than the total numbers of neurons $N$ s. Neglecting the variability of the contralateral couplings also constitutes a drastic approximation of the mean-field approach. This is all the more true that the average contralateral coupling $I$ happens to be small compared to its standard deviation.

Couplings are not only broadly distributed but also spatially organized. Ipsilateral couplings $J_{i j}$ decay with the distance between neurons $i$ and $j$ (Appendix 2—figure 4F). Similarly, contralateral couplings show strong correlations for short distances between the contralateral neurons (Appendix 2—figure 4G). The existence of a local spatial organization in the couplings is not unheard of in computational neuroscience and can have important functional consequences. It is, for instance, at the basis of ring-like attractor models and their extensions to two or three dimensions (Tsodyks and Sejnowski, 1995). Combined with the presence of variable biases $h_{i}$ , short-range interactions can lead to complex propagation phenomena, intensively studied in statistical physics in the context of the Random Field Ising Model. (Schneider and Pytte, 1977; Kaufman et al., 1986). As the most excitable neurons (with the largest biases) fire, they excite their neighbors, who in turn become active, triggering the activation of other neurons in their neighborhood. Such an avalanche mechanism could explain the fast rise of activity in the left or right region, from low- to high-activity state.

Interpretation of the functional connectivity

The inferred functional couplings $J_{i j}$ ’s are not expected to directly reflect the corresponding structural (synaptic) connectivity. However, their spatial distribution appears to be in line with the known ARTR organization (Dunn et al., 2016; Kinkhabwala et al., 2011) characterized by large positive (excitatory) interactions within the left and right population, and by the presence of negative (inhibitory) contralateral interactions. Although the contralateral couplings are found to be, on average, almost null, compared to the ipsilateral excitatory counterparts, they drive a subtle interplay between the left and right regions of the ARTR.

Our neural recordings demonstrate a systematic modulation of the ARTR dynamics with the water temperature, in quantitative agreement with the thermal dependence of the exploratory behavior in freely swimming assays. The model correctly captures this thermal modulation of the ARTR activity, and in particular the decay of the persistence time with the temperature. This owes to a progressive change in the values of both the couplings and the biases, which together deform the energy landscape and modulate the energy barriers between metastable states. The fact that the inferred functional connectivity between neurons does not display simple temperature dependence is not unexpected as different membrane currents can have different temperature dependence (Partridge and Connor, 1978).

In addition, as shown in Appendix 2—table 2, the inferred parameters largely vary across datasets. This variability is partially due to the difficulty to separately infer the interactions $J_{i j}$ and the biases $h_{i}$ , a phenomenon not specific to graphical model but also found with other neural, for example, Integrate-and-Fire network models (Monasson and Cocco, 2011). This issue can be easily understood within mean-field theory. For simplicity, let us neglect the weak contralateral coupling $I$ . The mean activity $m$ of a neuron then depends on the total ‘input’ $J m + H$ it receives, which is the sum of the bias $H$ and of the mean ipsilateral activity $m$ , weighted by the recurrent coupling $J$ . Hence, the combination $J m + H$ is more robustly inferred than $H$ and $J$ taken separately (Appendix 2—figure 7E).

The capacity to quantitatively capture subtle differences in the spontaneous activity induced by external cues is an important asset of our model. Recent studies have shown that spontaneous behavior in zebrafish larvae is not time-invariant but exhibits transitions between different regimes, lasting over minutes and associated with specific brain states. These transitions can have no apparent cause (Le Goc et al., 2021) or be induced by external (e.g., stimuli; Andalman et al., 2019) or internal cues (e.g., hunger states; Marques et al., 2019). Although they engage brain-wide changes in the pattern of spontaneous neural dynamics, they are often triggered by the activation of neuromodulatory centers such as the habenula-dorsal raphe nucleus circuit (Corradi and Filosa, 2021). Training Ising models in various conditions may help decipher how such neuromodulation impacts the network functional couplings leading to distinct dynamical regimes of spontaneous activity.

Data-driven modeling and metastability

With its slow alternating activity and relatively simple architecture, the ARTR offers an ideally suited circuit to test the capacity of Ising models to capture network-driven dynamics. The possibility to experimentally modulate the ARTR persistence time scale further enabled us to evaluate the model ability to quantitatively represent this slow process. The ARTR is part of a widely distributed hindbrain network that controls the eye horizontal saccadic movements, and which includes several other neuronal populations whose activity is tuned to the eye velocity or position (Joshua and Lisberger, 2015; Wolf et al., 2017). A possible extension of the model would consist in incorporating these nuclei in order to obtain a more complete representation of the oculomotor circuit. Beyond this particular functional network, a similar data-driven approach could be implemented to capture the slow concerted dynamics that characterize numerous neural assemblies in the zebrafish brain (van der Plas et al., 2021).

The importance of metastable states in cortical activity in mammals has been emphasized in previous studies as a possible basis for sequence-based computation (Harvey et al., 2012; Brinkman et al., 2022). Our model suggests that these metastable states are shaped by the connectivity of the network and are naturally explored during ongoing spontaneous activity. In this respect, the modification of the landscape resulting from visual stimulation, leading to a sharp decrease in the barrier separating the states, is reminiscent of the acceleration of sensory coding reported in Mazzucato et al., 2019. Our principled data-driven modeling could be useful to assess the generality of such metastable-state-based computations and of their modulation by sensory inputs in other organisms.

Materials and methods

All data and new codes necessary to reproduce the results reported in this work can be accessed at https://gin.g-node.org/Debregeas/ZF_ARTR_thermo and https://github.com/SebastWolf/ZF_ARTR_thermo.

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers
Strain, strain background (Danio rerio)	Tg(elavl3:H2B-GCaMP6s)	Vladimirov et al., 2014
Strain, strain background (D. rerio)	Tg(elavl3:H2B-GCaMP6f)	Quirin et al., 2016
Software, algorithm	Blind Sparse Deconvolution	Tubiana et al., 2020	BSD
Software, algorithm	Computational Morphometry Toolkit	https://www.nitrc.org/projects/cmtk/	CMTK
Software, algorithm	Adaptive Cluster Expansion	Barton and Cocco, 2013	ACE

Zebrafish lines and maintenance

Request a detailed protocol

All animals subjects were zebrafish (Danio rerio), aged 5–7 days post-fertilization (dpf). Larvae were reared in Petri dishes in embryo medium (E3) on a 14/10 hr light/dark cycle at 28°C, and were fed powdered nursery food (GM75) every day from 6 dpf.

Calcium imaging experiments were conducted on nacre mutants that were expressing either the calcium indicator GCaMP6f (12 fish) or GCaMP6s (1 fish) in the nucleus under the control of the nearly pan-neuronal promoter Tg(elavl3:H2B-GCaMP6). Both lines were provided by Misha Ahrens and published in Vladimirov et al., 2014 (H2B-GCaMP6s) and Quirin et al., 2016 (H2B-GCaMP6f).

All experiments were approved by the Le Comité d’Éthique pour l’Expérimentation Animale Charles Darwin (02601.01).

Behavioral assays

Request a detailed protocol

The behavioral experiments and preprocessing have been described in detail elsewhere (Le Goc et al., 2021). Shortly, it consists of a metallic pool regulated in temperature with two Peltier elements, recorded in uniform white light from above at 25 Hz. A batch of 10 animals experienced 30 min in water at either 18, 22, 26, 30, or 33°C (10 batches of 10 fish, involving 170 different individuals, were used). Movies were tracked with FastTrack (Gallois and Candelier, 2021), and MATLAB (The MathWorks) was used to detect discrete swim bouts from which the differences of orientation between two consecutive events are computed, referred to as turn or reorientation angles $δ θ$ .

Turn angles distributions could be fitted as the sum of two distributions (Gaussian and Gamma), whose intersection was used to define an angular threshold to categorize events into forward (F), left turn (L), or right turn (R, Figure 1E). This threshold was found to be close to 10° for all tested temperatures.

Then we ternarized $δ θ$ values, based on F, L, or R classification (Figure 1F), and computed the power spectrum of the binary signals defined from symbols L and R only, with the periodogram MATLAB function and averaged by temperature (Figure 1G). The outcome was fitted to the Lorentzian expression corresponding to a memory-less equiprobable two-state process (Odde and Buettner, 1998):

S (f) \propto \frac{2 k_{f l i p}}{4 k_{f l i p}^{2} + {(2 π f)}^{2}},

where $k_{f l i p}$ is the rate of transition from one state to another. The inverse of the fitted flipping rate $k_{f l i p}$ represents the typical time spent in the same orientational state, that is, the typical time taken to switch turning direction.

Light-sheet functional imaging of spontaneous activity

Request a detailed protocol

Volumetric functional recordings were carried out using custom-made one-photon light-sheet microscopes whose optical characteristics have been detailed elsewhere (Panier et al., 2013). Larvae were mounted in a 1 mm diameter cylinder of low melting point agarose at 2% concentration.

Imaged volume corresponded to 122 ± 46 μm in thickness, split into 16 ± 4 slices (mean ± SD). Recordings were of length 1392 ± 256 s with a brain volume imaging frequency of 6 ± 2 Hz (mean ± SD).

Image preprocessing, neurons segmentation, and calcium transient ( $Δ F / F$ ) extraction were performed offline using MATLAB, according to the workflow previously reported (Panier et al., 2013; Wolf et al., 2017; Migault et al., 2018).

A Peltier module is attached to the lower part of the pool (made of tin) with thermal tape (3M). A type T thermocouple (Omega) is placed near the fish head (<5 mm) to record the fish surrounding temperature. The signal from a thermocouple amplifier (Adafruit) is used in a PID loop implemented on an Arduino board, which mitigate the Peltier power to achieve the predefined temperature target, stable at ±0.5°C. The temperature regulation software and electronics design are available on Gitlab under a GNU GPLv3 license (https://gitlab.com/GuillaumeLeGoc/arduino-temperature-control copy archived at Le Goc, 2022).

The ARTR neurons were selected using a method described elsewhere (Wolf et al., 2017). First, a group of neurons was manually selected on a given slice based on a morphological criterion such that the ARTR structure (ipsilateral correlations and contralateral anticorrelation) is revealed. Then, neurons showing Pearson’s correlation (anti-correlation) higher than 0.2 (less than –0.15, respectively) are selected, manually filtering them on a morphological criterion. Those neurons are then added to the previous ones, whose signals are used to find neurons from the next slice and so on until all slices are treated.

For fish that were recorded at different temperatures, to ensure that the same neurons are selected, we used the Computational Morphometry Toolkit (CMTK, https://www.nitrc.org/projects/cmtk/) to align following recordings onto the first one corresponding to the same individual. Resulting transformations are then applied to convert neurons coordinates in a consistent manner through all recordings involving the same fish.

Visually driven recordings

Request a detailed protocol

Volumetric functional recordings under visual stimulation were carried using our two-photon light-sheet microscope described in Wolf et al., 2015. The stimulation protocol was previously explained in Wolf et al., 2017: two LEDs were positioned symmetrically outside of the chamber at 45° and 4.5 cm from the fish eyes, delivering a visual intensity of 20 μW/cm². We alternately illuminated 17 times each eye for 10 s, 15 s, 20 s, 25 s, and 30 s while performing two-photon light-sheet brain-wide functional imaging. Synchronization between the microscope and the stimulation set-up was done using a D/A card (NI USB-6259 NCS, National Instruments) and a LabVIEW program. Brain volume image frequency was of 1 Hz on the six recorded fish. Recordings last for 4500 s, 856 s of which is spontaneous activity. We extracted the ARTR neurons following the same procedure described above, yielding 89 ± 54 neurons (mean ± SD).

Time constants definitions

Request a detailed protocol

For the flipping rates (Figure 1D), we defined the time-dependent signed activity of the ARTR (Figure 1B) through

σ (t) = sign (m_{L} (t) - m_{R} (t)),

where $m_{L, R} (t) = \frac{1}{N_{L, R}} \sum_{i \in L, R} s_{i} (t)$ are the average activities in the L, R regions. A power spectrum density is estimated for each signal with the Thomson’s multitaper method through the pmtm MATLAB function (time-halfbandwidth product set to 4). The power spectrum densities were then fitted with a Lorentzian spectrum see Equation 6 and Figure 1G.

ARTR left and right persistence times (Figure 3B) are defined as the time $m_{L}$ and $m_{R}$ signals spend consecutively above an arbitrary threshold set at 0.1. Left and right signals are treated altogether. Changing the threshold does induce a global offset but does not change the observed effect of temperature, the relation with $m_{L}$ and $m_{R}$ mean signals, nor the relation with the persistence times of the synthetic signals. The persistence times of the synthetic signals, generated with the Ising models, are computed using the same procedure: we compute the time $m_{L}$ and $m_{R}$ synthetic signals spend consecutively above an arbitrary threshold set at 0.1, we then normalize these durations by the corresponding experimental frame rate in order to compare the different recordings (Figure 3B). For the mean-field simulated dynamics of the left and right activities, we also follow the same strategy in order to compute the persistence times displayed in Figure 4F.

Inference of Ising model from neural activity

From spontaneous activity to spiking data, to biases and connectivity

Request a detailed protocol

For each recording (animal and/or temperature), approximate spike trains were inferred from the fluorescence activity signal using the Blind Sparse Deconvolution algorithm (Tubiana et al., 2020). This algorithm features automatic (fully unsupervised) estimation of the hyperparameters, such as spike amplitude, noise level, and rise and decay time constants, but also an automatic thresholding for binarizing spikes such as to maximize the precision-recall performance. The binarized activity of the $N$ recorded neurons was then described for each time bin $t$ , into a $N$ -bit binary configuration $s_{t}$ , with, $s_{i} (t) = 1$ if neuron $i$ is active in bin $t$ , 0 otherwise.

The functional connectivity matrix $J_{ij}$ and the biases $h_{i}$ defining the Ising probability distribution over neural configurations (see Equation 1) were determined such that the pairwise correlations and average activities computed from the model match their experimental counterparts. In practice, we approximately solved this hard inverse problem using the Adaptative Cluster Expansion and the MC learning algorithms described in Cocco and Monasson, 2011 and in Barton and Cocco, 2013. The full code of the algorithms can be downloaded from the GitHub repository: https://github.com/johnbarton/ACE/ ( Barton, 2019).

Monte Carlo sampling

Request a detailed protocol

In order to generate synthetic activity, we resorted to Gibbs sampling, a class of Monte Carlo Markov Chain method, also known as Glauber dynamics. At each time step $k$ , a neuron, say, $i$ , is picked up uniformly at random, and the value of its activity is updated from $s_{i}^{k}$ to $s_{i}^{k + 1} = 0, 1$ according to the probability

P (s_{i}^{k + 1} ∣ s_{j \neq i}^{k}) = \frac{\exp (s_{i}^{k + 1} (h_{i} + \sum_{j} J_{i j} s_{j}^{k}))}{1 + \exp (h_{i} + \sum_{j} J_{i j} s_{j}^{k})}

which depends on the current activities of the other neurons. As this updating fulfills detailed balance, the probability distribution of $s^{k}$ eventually converges to $P$ in Equation 1. A Monte Carlo round is defined as the number of Monte Carlo steps divided by the total number of neurons, $N$ . The code used can be accessed from the link provided at the beginning of the ‘Materials and methods’ section.

Cross-validation and independent model

Request a detailed protocol

We cross-validated the Ising models (see Appendix 2—figure 3) dividing the datasets in two parts: for each experiment, 75% of each dataset is used as a training set and the remaining 25% is used as a test set. Each training set is used to infer an Ising model. We then compare the mean activity and covariance of the test set with the one computed from the simulated data generated by the models (Appendix 2—figure 3A and B). We also show the relative variation of the models’ log likelihood computed on the training data and the test data (Appendix 2—figure 3C). In addition, as a null hypothesis, we decided to compare the Ising models fitted on the data with the independent model. The independent model depends on the mean activities ${⟨ s_{i} ⟩}_{data}$ only and reads

P (s) = \frac{1}{Z} \exp (\sum_{i} h_{i} s_{i}),

We demonstrate in Appendix 2—figure 3E–F the inefficiency of the independent models, comparing the mean activity and covariance of the test set with the one computed from the simulated data generated by the independent models. We also show the relative variation, between the Ising and the independent models, of the log likelihood computed on the training data and the test data (Appendix 2—figure 3G).

Real data and models comparison

Request a detailed protocol

To quantify the quality of the log-probability landscapes reproduction by the Ising models (Figure 3A), we used the Kullback–Leibler divergence between (1) a dataset $i$ and the synthetic signals generated with the model trained on that dataset $i$ (green) and (2) the dataset $i$ with synthetic signals generated with every other models (red). With c_i the count in the two-dimensional bin $i$ (10 × 10 bins used) and $α$ a pseudocount (set to 1), the probability in bin $i$ is defined as $P_{i} = \frac{c_{i} + α}{\sum_{j} (c_{j} + α)}$ . The Kullback–Leibler divergence between a data/model pair is then defined as

D_{K L} = \sum_{i} P_{d a t a, i} \log_{10} (\frac{P_{d a t a, i}}{P_{m o d e l, i}}) .

We follow the exact same procedure in order to compare the independent model and their corresponding datasets (Figure 3A in blue). In this case, we use synthetic signals generated with the independent model to define $P_{m o d e l, i}$ .

Inference of additional biases from visually driven activity recordings

Request a detailed protocol

For the visually driven activity recordings, we infer the additional biases $δ {\overset{\leftarrow}{h}}_{i}$ from the recordings of the ARTR activity (Figure 5D) during, for example, the leftward light stimulations as follows. Let $\overset{\leftarrow}{B}$ the number of time bins $t = 1, 2, \dots, \overset{\leftarrow}{B}$ in the recording, and $s_{t}$ the corresponding binarized activity configurations. We define, for each neuron $i$ ,

ρ_{i} (δ h) = \sum_{t = 1}^{\overset{\leftarrow}{B}} \frac{\exp (h_{i} + \sum_{j} J_{i j} s_{j} (t) + δ h)}{1 + \exp (h_{i} + \sum_{j} J_{i j} s_{j} (t) + δ h)} .

$ρ_{i} (δ h)$ represents the mean activity of neuron $i$ , when subject to a global bias summing $h_{i}$ , the other neurons activities $s_{j} (t)$ weighted by the couplings $J_{i j}$ , and an additional bias $δ h$ , averaged over all the frames $t$ corresponding to left-sided light stimulation. It is a monotonously increasing function of $δ h$ , which matches the experimental average activity $\frac{1}{\overset{\leftarrow}{B}} \sum_{t = 1}^{\overset{\leftarrow}{B}} s_{i} (t)$ for a unique value of its argument. This value defines $δ {\overset{\leftarrow}{h}}_{i}$ . The same procedure was followed to infer the additional biases $δ {\vec{h}}_{i}$ associated to rightward visual stimulations.

Appendix 1

Mean-field theory for the ARTR activity

Derivation of the free energy

We consider an Ising model with N_L and N_R neurons in, respectively, the left and right regions. Each neuron activity variable can take two values, i = 0, 1, corresponding to silent and active states (within a time window). The ‘energy’ of the system reads

E (s_{1}, \dots, s_{N_{L}}, s_{N_{L} + 1}, \dots, s_{N_{L} + N_{R}}) = - {\tilde{H}}_{L} \sum_{i = 1}^{N_{L}} s_{i} - {\tilde{H}}_{R} \sum_{i = N_{L} + 1}^{N_{L} + N_{R}} s_{i} - \frac{1}{2} \sum_{i \neq j} {\tilde{J}}_{i j} s_{i} s_{j},

where ${\tilde{H}}_{L}, {\tilde{H}}_{R}$ are biases acting on the neurons, and the coupling matrix is defined through

\begin{array}{ll} {\tilde{J}}_{i j} = {\begin{cases} {\tilde{J}}_{L} & i f & 1 \leq i, j \leq N_{L}, \\ {\tilde{J}}_{R} & i f & N_{L} + 1 \leq i, j \leq N_{L} + N_{R}, \\ \tilde{I} & o t h e r w i s e \end{cases} \end{array}

We now introduce the left and right average activities:

m_{L} = \frac{1}{N_{L}} \sum_{i = 1}^{N_{L}} s_{i}, m_{R} = \frac{1}{N_{R}} \sum_{i = N_{L} + 1}^{N_{L} + N_{R}} s_{i} .

The energy $E$ of a neural activity configuration in Equation 12 can be expressed in terms of these average activities:

\begin{array}{ll} E (m_{L}, m_{R}) & = - N_{L} ({\tilde{H}}_{L} - \frac{{\tilde{J}}_{L}}{2}) m_{L} - N_{R} ({\tilde{H}}_{R} - \frac{{\tilde{J}}_{R}}{2}) m_{R} \\ - \frac{(N_{L})^{2}}{2} {\tilde{J}}_{L} m_{L}^{2} - \frac{(N_{R})^{2}}{2} {\tilde{J}}_{R} m_{R}^{2} - \tilde{I} N_{L} N_{R} m_{L} m_{R} . \end{array}

We may now compute the partition function normalizing the probability of configurations,

Z = \sum_{{s_{i} = 0, 1}} e^{- E (s_{1}, \dots, s_{N_{L} + N_{R}})} = \sum_{m_{L}, m_{R}} M_{L} (m_{L}) M_{R} (m_{R}) e^{- E (m_{L}, m_{R})},

where the sums run over fractional values of the average left and right activities, from 0 to 1 with steps equal to, respectively, $2 / N_{L}$ and $2 / N_{R}$ , and the multiplicities $ℳ_{L}$ and $ℳ_{R}$ measure the numbers of neural configurations with prescribed average activities. We approximate these multiplicities with the standard entropy-based expressions, which are exact in the limit of large sizes $K_{L}$ , $K_{R}$ :

M_{L} (m_{L}) ≃ e^{N_{L} S (m_{L})}, M_{R} (m_{R}) ≃ e^{N_{R} S (m_{R})},

where

S (m) = - m \ln m - (1 - m) \ln (1 - m)

is the entropy of a 0 − 1 variable with mean $m$ . As a consequence, the activity-dependent free energy is given by

\begin{array}{ll} F (m_{L}, m_{R}) & = E (m_{L}, m_{R}) - N_{L} S (m_{L}) - N_{R} S (m_{R}) \\ = - \frac{N_{L} J_{L}}{2} {m_{L}}^{2} - \frac{N_{R} J_{R}}{2} {m_{R}}^{2} - I \sqrt{N_{L} N_{R}} m_{L} m_{R} - N_{L} H_{L} m_{L} - N_{R} H_{R} m_{R} \\ + N_{L} (m_{L} \ln m_{L} + (1 - m_{L}) \ln (1 - m_{L})) + N_{R} (m_{R} \ln m_{R} + (1 - m_{R}) \ln (1 - m_{R})) \end{array}

where the bias and coupling parameters are, respectively, $H_{L} = {\tilde{H}}_{L} - \frac{{\tilde{J}}_{L}}{2}$ , $H_{R} = {\tilde{H}}_{R} - \frac{{\tilde{J}}_{R}}{2}$ , $J_{L} = N_{L} {\tilde{J}}_{L}$ , $J_{R} = N_{R} {\tilde{J}}_{R}$ , $I = \sqrt{N_{L} N_{R}} \tilde{I}$ .

The sizes $N_{L}, N_{R}$ enter formula (19) for the free energy in two ways:

Implicitly, through the biases $H_{L}, H_{R}$ and the couplings $J_{L}, J_{R}, I$ . These parameters are equal to, respectively, the average bias and the total ipsilateral and contralateral couplings acting on each neuron in the $L$ and $R$ regions. They are effective parameters defining the mean-field theory.
Explicitly, as multiplicative factors to the free energy contributions coming from the left and right regions. The sizes then merely act as effective inverse ‘temperatures,’ in the Boltzmann factor $e^{- F (m_{L}, m_{R})}$ associated to the probability of the $L, R$ activities.

Mean-field theory generally overestimates the collective effects of interactions; a well-known illustration of this artifact is the prediction of the existence of a phase transition in the unidimensional ferromagnetic Ising model with short-range interactions, while such a transition is rigorously known not to take place (Ma, 1985). We expect these effects to be strong here due to the wide distribution of inferred Ising couplings (Appendix 2—figure 4A). Many pairs of neurons carry close to zero couplings, and the interaction neighborhood of a neuron is effectively much smaller than $N_{L}$ and $N_{R}$ . To compensate for the overestimation of interaction effects, we thus propose to keep Equation 19 for the free energy, but with effective sizes $K_{L}, K_{R}$ replacing the numbers $N_{L}, N_{R}$ of recorded neurons (see Equation 2), leading to the expression of the free energy:

\begin{array}{ll} F (m_{L}, m_{R}) & = - \frac{K_{L} J_{L}}{2} m_{R}^{2} - \frac{K_{R} J_{R}}{2} m_{R}^{2} - I \sqrt{K_{L} K_{R}} m_{L} m_{R} - K_{L} H_{L} m_{L} - K_{R} H_{R} m_{R} \\ + K_{L} (m_{L} 1 n m_{L} + (1 - m_{L}) 1 n (1 - m_{L})) + K_{R} (m_{R} 1 n m_{R} + (1 - m_{R}) 1 n (1 - m_{R})) \end{array}

These effective sizes $K_{L}, K_{R}$ are expected to be smaller than $N_{L}, N_{R}$ . Their values are fixed through the comparison of the Langevin dynamical traces with the traces coming from the data; see below.

Langevin dynamical equations

The dynamical Langevin equations read

τ \frac{d m_{L}}{d t} = K_{L} (J_{L} m_{L} + H_{L}) + I \sqrt{K_{L} K_{R}} m_{R} - K_{L} \log (\frac{m_{L}}{1 - m_{L}}) + ϵ_{L} (t),

τ \frac{d m_{R}}{d t} = K_{R} (J_{R} m_{R} + H_{R}) + I \sqrt{K_{L} K_{R}} m_{L} - K_{R} \log (\frac{m_{R}}{1 - m_{R}}) + ϵ_{R} (t),

where $ϵ_{L}, ϵ_{R}$ denote white-noise processes; see main text.

Fit of the effective sizes $K_{L}$ and $K_{R}$

The effective sizes $K_{L} = N_{L} / A$ and $K_{R} = N_{R} / A$ were fitted generating Langevin trajectories of the activities ( $m_{L}$ , $m_{R}$ ) for a large set of values of $A$ (i.e., $K_{L}$ and $K_{R}$ ), and with fixed parameters ( $H_{L}$ , $H_{R}$ , $J_{L}$ , $J_{R}$ , $τ$ ). For each value of $K_{L}$ and $K_{R}$ , we computed the KL divergence between the experimental and the Langevin distributions of ( $m_{L}$ , $m_{R}$ ) (see Appendix 2—figure 7A–C). The effective sizes $K_{L}$ and $K_{R}$ are the ones that minimize the value of the KL divergence. For low values of $A$ , the KL divergence can be noisy and creates artifacts. To avoid these artifacts, we assume that $A > 2$ .

Appendix 2

Appendix 2—figure 1

Download asset Open asset

Temperature dependence of the anterior rhombencephalic turning region (ARTR) activity.

(A) Schematic of the experimental setup used to perform brain-wide calcium imaging of a zebrafish larva at controlled water temperature. (B) Raster plot of the ARTR spontaneous dynamics showing alternating right/left activation. The top and bottom traces are the ARTR average signal of the left and right subcircuits. (C) Example ARTR sign( $m_{L} - m_{R}$ ) binarized signals measured at three different temperatures (same larva). (D) Averaged power spectrum of the ARTR signals $m_{R} - m_{L}$ for the five tested temperatures. Lorentzian fits are shown as solid lines.

Appendix 2—figure 2

Download asset Open asset

Effect of temperature on the anterior rhombencephalic turning region (ARTR) time persistence and activity.

(A) Pdf of activities of both sides of the ARTR. Color encodes temperature. (B) Temperature-averaged mean activity of ARTR left and right neuronal subpopulations. Error bars are standard error of the mean. (C) Temperature-averaged Pearson correlation for left/right ispilateral pairs (yellow line) or for contralateral pairs of neurons (purple line). Error bars are standard deviations (32 recordings), n = 13 fish at 5 different water temperatures. (D) ARTR persistence time vs. mean activity; note the quasi-linear dependence of these quantities ( $R = 0.91$ ). Each dot is the mean persistence time computed for one fish at one temperature; colors encode temperature.

Appendix 2—figure 3

Download asset Open asset

Inference of the anterior rhombencephalic turning region (ARTR) Ising model.

(**A, B**) Comparison between the mean activities (A) and pairwise correlations (B) computed from experimental test data and from synthetic (Ising model-generated) data (32 recordings, n = 13 fish). Ising models were trained on a distinct subset of the experimental data. (C) Relative variation of the log–likelihoods of the Ising models between training and test data, showing the absence of overfitting. (D) Probability that $K$ of the $N$ neurons in the ARTR are simultaneously active in the data (black dots) and in the model (yellow line) configurations. (**E, F**) In order to demonstrate the need for effective connections in our model, we generated synthetic data with independent models of the training dataset. Here, we compare the mean activity (E) and the pairwise covariance (F) computed on the experimental test dataset and using independent models. (G) Excess log likelihood of the Ising models compared to the independent model for training and test data set (see ‘Materials and methods’).

Appendix 2—figure 4

Download asset Open asset

Correlation structure within the anterior rhombencephalic turning region (ARTR) and properties of the inferred couplings.

(A) Probability density function of the functional connectivity for the ipsilateral (gold line) and the contralateral (purple line) couplings. These pdf were obtained by averaging across all animals. (B) Probability density function of the functional Pearson correlation for the ipsilateral (gold line) and the contralateral (purple line) couplings. (C) Box plot across experiments of the average value of the ipsilateral and contralateral couplings. (D) Probability to have an ipsilateral (gold line) or a contralateral (purple line) pair of neuron given its effective connectivity. For a given range of the effective connectivity, we compute the number of ipsilateral and contralateral pairs of neurons. (E) Probability to have an ipsilateral (gold line) or a contralateral (purple line) pair of neuron given its Pearson correlation. (F) Functional connectivity $J_{i j}$ as a function of the distance between neurons $i, j$ . (G) Correlation between the couplings $J_{i j}$ and $J_{k p}$ , between one neuron $i$ and one neuron $k$ as a function of their distance $d_{i k}$ for every possible pair $(i, k)$ .

Appendix 2—figure 5

Download asset Open asset

Distribution of biases in the inferred anterior rhombencephalic turning region (ARTR) Ising model.

(A) Bias parameter distribution for an example fish. (B) Box plot across experiments of the average value of the biases for the left and right subpopulations of the ARTR. (C) Box plot across animals of the standard deviation of the biases for the left and right subpopulations of the ARTR.

Appendix 2—figure 6

Download asset Open asset

Correlation of Ising parameters at different temperatures.

For each fish (n = 13), we extract from the scatter plots of the coupling $J_{i j}$ and bias h_i inferred from activity recordings at two different temperatures, the Pearson correlation coefficients $R_{p e a r s o n}$ . The distribution of $R_{p e a r s o n}^{2}$ values are shown for all fish and pairs of temperature. Inset: Example scatter plots of the inferred biases h_i (left) and effective couplings $J_{i j}$ (right) for the same fish at two different temperature T = 22 and T = 30°C.

Appendix 2—figure 7

Download asset Open asset

Mean-field model of the anterior rhombencephalic turning region (ARTR).

(**A, B**) Kullback–Leibler divergence between the experimental and the Langevin distributions as a function of $N / K$ , where $N$ is the total number of neurons of the left or right subpopulation, and $K$ is the effective extent of neuronal interaction (see ‘Materials and methods’) for two datasets. (C) Probability density function of $K_{R}$ (blue line) and $K_{L}$ (red line) across all recordings. (D) Free-energy difference between stationary sates of the landscape as a function of the temperature (32 recordings, n=13 fish). (E) Average values (for all experiments and regions) of $K (H + J M)$ as a function of the temperature of the water. Error bars are standard error of the mean.

Appendix 2—figure 8

Download asset Open asset

A modified Ising model explains visually driven properties of the anterior rhombencephalic turning region (ARTR).

(**A, B**) To assess the performance of the model for visually driven experiments, we compare the mean activity (A) and the pairwise covariance (B) computed on the spontaneous part of the recordings to synthetic data. (C) Scatter plot of the correlation between contralateral pairs of neurons under visual stimulation vs. spontaneous activity on n = 6 fish. (D) Scatter plot of the correlation between ipsilateral pairs of neurons under visual stimulation vs. spontaneous activity. (E) Average Pearson correlation in the experimental recordings between contralateral (the p-value of a paired sampled t-test is provided) and ipsilateral pairs of cells during stimulated and spontaneous activity (n = 6 fish). (F) Average Pearson correlation in the simulated activity of the ARTR between contralateral and ipsilateral pairs of cells during stimulated and spontaneous activity (n = 6 fish).

Appendix 2—table 1

Datasets properties.

Temperature (°C)	ID	Line	Age (dpf)	$N_{L}$	$N_{R}$	Acquisition rate (Hz)	Duration (s)
18	12	NucFast	6	146	180	5	1200
18	13	NucFast	7	37	96	8	1200
18	14	NucFast	6	179	174	8	1200
22	2	Nuc slow	7	177	212	3	1106
22	3	NucFast	5	152	85	3	1812
22	5	NucFast	5	158	123	5	1500
22	6	NucFast	5	98	134	5	1500
22	7	NucFast	6	122	221	5	1500
22	11	NucFast	6	295	320	5	1200
22	13	NucFast	7	37	96	8	1200
22	14	NucFast	6	179	174	8	1200
26	2	Nuc slow	7	177	212	3	1812
26	3	NucFast	5	152	85	3	1812
26	4	NucFast	5	110	76	3	1812
26	5	NucFast	5	158	123	5	1500
26	6	NucFast	5	98	134	5	1500
26	7	NucFast	6	122	221	5	1500
26	11	NucFast	6	295	320	5	1200
26	13	NucFast	7	37	96	8	1200
26	14	NucFast	6	179	174	8	1200
30	2	Nuc slow	7	177	212	3	1812
30	4	NucFast	5	110	76	3	1812
30	5	NucFast	5	158	123	5	1500
30	6	NucFast	5	98	134	5	1500
30	7	NucFast	6	122	221	5	1500
30	13	NucFast	7	37	96	8	1200
30	14	NucFast	6	179	174	8	1200
30	15	NucFast	7	202	252	8	1200
33	14	NucFast	6	179	174	8	1200
33	15	NucFast	7	202	252	8	1200
33	16	NucFast	6	127	123	7	1200
33	17	NucFast	5	62	170	10	1200

Appendix 2—table 2

Parameters of mean-field models.

Temperature (°C)	ID	$J_{L}$	$J_{R}$	$I$	$H_{L}$	$H_{R}$	$K_{L}$	$K_{R}$
18	12	7.06	7.23	–0.6	–3.66	–3.63	6.51	8.03
18	13	6.2	7.84	0.6	–3.53	–4.34	3.18	8.27
18	14	7.27	7.24	0.31	–3.88	–3.99	11.04	10.74
22	2	8.2	8.28	0.12	–4.24	–4.23	6.65	7.96
22	3	8.18	7.14	0.55	–4.26	–4.13	9.38	5.24
22	5	7.59	7.01	0.4	–4.03	–3.8	5.56	4.33
22	6	7.13	8.69	1.1	–4.49	–4.64	5.21	7.12
22	7	7.09	7.46	0.43	–3.73	–3.95	6.28	11.39
22	11	7.82	7.59	–0.1	–4.07	–3.91	8.28	8.98
22	13	6.54	7.82	1.45	–4.29	–4.5	7.11	18.46
22	14	7.41	8.03	0.47	–4.28	–4.43	10.91	10.6
26	2	8.37	8.22	–0.49	–4.47	–4.31	9.72	11.64
26	3	8.42	7.49	0.53	–4.56	–4.62	8.26	4.61
26	4	8.63	6.44	0.85	–4.83	–4.79	10.37	7.16
26	5	7.29	7.59	0.48	–3.92	–4.14	9.08	7.06
26	6	7.43	7.86	0.41	–3.99	–4.1	8.59	11.75
26	7	7.55	7.96	0.32	–4.08	–4.22	4.45	8.06
26	11	7.27	7.45	0.37	–3.89	–3.92	10.31	11.18
26	13	6.99	7.3	0.6	–3.99	–3.94	6.37	16.55
26	14	7.91	7.35	0.5	–4.34	–4.16	11.32	11.01
30	2	7.54	7.96	–0.12	–4.54	–4.56	7.02	8.41
30	4	8.36	7.73	0.11	–4.52	–4.18	9.64	6.66
30	5	6.77	6.42	0.66	–3.8	–3.87	9.18	7.15
30	6	7.35	7.38	0.45	–3.91	–3.97	7.53	10.3
30	7	7.43	8.07	0.42	–3.93	–4.38	7.09	12.84
30	13	6.91	7.41	0.73	–4.13	–4.03	5.78	15
30	14	7.51	7.45	0.11	–3.87	–3.89	9.42	9.15
30	15	8.01	8.33	0.58	–4.45	–4.46	13.83	17.26
33	14	6.74	7.02	0.76	–3.8	–3.97	9.32	9.06
33	15	6.99	7.47	–0.02	–3.68	–3.91	14.85	18.52
33	16	7.53	8.25	–0.11	–4.16	–4.43	14.43	13.97
33	17	6.66	7.36	0.45	–3.69	–3.89	11.92	32.69

Appendix 2—table 3

Parameters of mean-field models.

ID	$J_{L}$	$J_{R}$	$I$	$H_{L}$	$H_{R}$	$K_{L}$	$K_{R}$
1	7.54	7.35	–0.67	–3.75	–3.44	5.60	3.43
2	7.10	7.42	0.64	–3.69	–4.02	7.91	12.82
3	7.51	7.92	–0.28	–3.96	–4.08	4.98	3.90
4	8.38	6.25	–0.04	–3.68	–3.18	13.33	4.44
5	8.73	8.24	0.01	–4.38	–4.13	6.11	6.89
6	7.87	7.71	0.51	–4.17	–4.09	16.19	15.52

Data availability

All data necessary to reproduce the results reported in this work can be accessed from https://gin.g-node.org/Debregeas/ZF_ARTR_thermo. Codes can be accessed at https://github.com/SebastWolf/ZF_ARTR_thermo copy archived at Wolf, 2023.

The following data sets were generated

1. Wolf S
2. Le Goc G
3. Debregeas G
4. Cocco S
5. Monasson R
(2023) G-Node GIN
ID Debregeas/ZF_ARTR_thermo. ZF_ARTR_thermo.

https://gin.g-node.org/Debregeas/ZF_ARTR_thermo

References

1. Ahrens MB
2. Orger MB
3. Robson DN
4. Li JM
5. Keller PJ
(2013) Whole-brain functional imaging at cellular resolution using light-sheet microscopy
Nature Methods 10:413–420.
https://doi.org/10.1038/nmeth.2434
- PubMed
- Google Scholar
1. Andalman AS
2. Burns VM
3. Lovett-Barron M
4. Broxton M
5. Poole B
6. Yang SJ
7. Grosenick L
8. Lerner TN
9. Chen R
10. Benster T
11. Mourrain P
12. Levoy M
13. Rajan K
14. Deisseroth K
(2019) Neuronal dynamics regulating brain and behavioral state transitions
Cell 177:970–985.
https://doi.org/10.1016/j.cell.2019.02.037
- PubMed
- Google Scholar
1. Barton J
2. Cocco S
(2013) Ising models for neural activity inferred via selective cluster expansion: structural and coding properties
Journal of Statistical Mechanics 2013:03002.
https://doi.org/10.1088/1742-5468/2013/03/P03002
- Google Scholar
(2016) ACE: adaptive cluster expansion for maximum entropy graphical model inference
Bioinformatics 32:3089–3097.
https://doi.org/10.1093/bioinformatics/btw328
- PubMed
- Google Scholar
Software
1. Barton JP
(2019) Ace, version 8032305
GitHub.

https://github.com/johnbarton/ACE/
1. Brinkman BAW
2. Yan H
3. Maffei A
4. Park IM
5. Fontanini A
6. Wang J
7. La Camera G
(2022) Metastable dynamics of neural circuits and networks
Applied Physics Reviews 9:011313.
https://doi.org/10.1063/5.0062603
- PubMed
- Google Scholar
1. Butts DA
(2019) Data-Driven approaches to understanding visual neuron activity
Annual Review of Vision Science 5:451–477.
https://doi.org/10.1146/annurev-vision-091718-014731
- PubMed
- Google Scholar
1. Chen X
2. Randi F
3. Leifer AM
4. Bialek W
(2019) Searching for collective behavior in a small brain
Physical Review. E 99:052418.
https://doi.org/10.1103/PhysRevE.99.052418
- PubMed
- Google Scholar
(2009) Neuronal couplings between retinal ganglion cells inferred by efficient inverse statistical physics methods
PNAS 106:14058–14062.
https://doi.org/10.1073/pnas.0906705106
- PubMed
- Google Scholar
1. Cocco S
2. Monasson R
(2011) Adaptive cluster expansion for inferring Boltzmann machines with noisy data
Physical Review Letters 106:090601.
https://doi.org/10.1103/PhysRevLett.106.090601
- PubMed
- Google Scholar
1. Corradi L
2. Filosa A
(2021) Neuromodulation and behavioral flexibility in larval zebrafish: from neurotransmitters to circuits
Frontiers in Molecular Neuroscience 14:718951.
https://doi.org/10.3389/fnmol.2021.718951
- PubMed
- Google Scholar
1. Dunn TW
2. Mu Y
3. Narayan S
4. Randlett O
5. Naumann EA
6. Yang CT
7. Schier AF
8. Freeman J
9. Engert F
10. Ahrens MB
(2016) Brain-Wide mapping of neural activity controlling zebrafish exploratory locomotion
eLife 5:e12741.
https://doi.org/10.7554/eLife.12741
- PubMed
- Google Scholar
1. Gallois B
2. Candelier R
(2021) FastTrack: an open-source software for tracking varying numbers of deformable objects
PLOS Computational Biology 17:e1008697.
https://doi.org/10.1371/journal.pcbi.1008697
- PubMed
- Google Scholar
1. Gau P
2. Poon J
3. Ufret-Vincenty C
4. Snelson CD
5. Gordon SE
6. Raible DW
7. Dhaka A
(2013) The zebrafish ortholog of TRPV1 is required for heat-induced locomotion
The Journal of Neuroscience 33:5249–5260.
https://doi.org/10.1523/JNEUROSCI.5403-12.2013
- PubMed
- Google Scholar
(2020) Machine learning for neural decoding
ENeuro 7:ENEURO.0506-19.2020.
https://doi.org/10.1523/ENEURO.0506-19.2020
- PubMed
- Google Scholar
1. Guo ZV
2. Inagaki HK
3. Daie K
4. Druckmann S
5. Gerfen CR
6. Svoboda K
(2017) Maintenance of persistent activity in a frontal thalamocortical loop
Nature 545:181–186.
https://doi.org/10.1038/nature22324
- PubMed
- Google Scholar
1. Harvey CD
2. Coen P
3. Tank DW
(2012) Choice-specific sequences in parietal cortex during a virtual-navigation decision task
Nature 484:62–68.
https://doi.org/10.1038/nature10918
- PubMed
- Google Scholar
1. Jaynes ET
(1957) Information theory and statistical mechanics
Physical Review 106:620–630.
https://doi.org/10.1103/PhysRev.106.620
- Google Scholar
1. Joshua M
2. Lisberger SG
(2015) A tale of two species: neural integration in zebrafish and monkeys
Neuroscience 296:80–91.
https://doi.org/10.1016/j.neuroscience.2014.04.048
- PubMed
- Google Scholar
(1986) Multicritical points in an Ising random-field model
Physical Review. B, Condensed Matter 34:4766–4770.
https://doi.org/10.1103/physrevb.34.4766
- PubMed
- Google Scholar
1. Kinkhabwala A
2. Riley M
3. Koyama M
4. Monen J
5. Satou C
6. Kimura Y
7. Higashijima S-I
8. Fetcho J
(2011) A structural and functional ground plan for neurons in the hindbrain of zebrafish
PNAS 108:1164–1169.
https://doi.org/10.1073/pnas.1012185108
- PubMed
- Google Scholar
Book
1. Koller D
2. Friedmann N
(2009)
Probabilistic Graphical Models: Principles and Techniques

MIT Press.
- Google Scholar
1. Langer JS
(1969) Statistical theory of the decay of metastable states
Annals of Physics 54:258–275.
https://doi.org/10.1016/0003-4916(69)90153-5
- Google Scholar
(2021) Thermal modulation of zebrafish exploratory statistics reveals constraints on individual behavioral variability
BMC Biology 19:208.
https://doi.org/10.1186/s12915-021-01126-w
- PubMed
- Google Scholar
Software
1. Le Goc G
(2022) Arduino temperature control, version swh:1:rev:9327f2fb3dd1f1cd844c35712299a8929a084b6e
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:da1f130cf0dee06a154abbe18703802994b56570;origin=https://gitlab.com/GuillaumeLeGoc/arduino-temperature-control.git;visit=swh:1:snp:a1187e6e48df7ddb93de7817883c4c1a545a62c4;anchor=swh:1:rev:9327f2fb3dd1f1cd844c35712299a8929a084b6e
(2021) A distributed saccade-associated network encodes high velocity conjugate and monocular eye movements in the zebrafish hindbrain
Scientific Reports 11:12644.
https://doi.org/10.1038/s41598-021-90315-2
- PubMed
- Google Scholar
1. Long MA
2. Fee MS
(2008) Using temperature to analyse temporal dynamics in the songbird motor pathway
Nature 456:189–194.
https://doi.org/10.1038/nature07448
- PubMed
- Google Scholar
Book
1. Ma SK
(1985) Statistical Mechanics
World Scientific.
https://doi.org/10.1142/0073
- Google Scholar
1. Marques JC
2. Li M
3. Schaak D
4. Robson DN
5. Li JM
(2019) Internal state dynamics shape brainwide activity and foraging behaviour
Nature 577:239–243.
https://doi.org/10.1038/s41586-019-1858-z
- PubMed
- Google Scholar
(2009) Prediction of spatiotemporal patterns of neural activity from pairwise correlations
Physical Review Letters 102:138101.
https://doi.org/10.1103/PhysRevLett.102.138101
- PubMed
- Google Scholar
(2019) Expectation-induced modulation of metastable activity underlies faster coding of sensory stimuli
Nature Neuroscience 22:787–796.
https://doi.org/10.1038/s41593-019-0364-9
- PubMed
- Google Scholar
1. Meshulam L
2. Gauthier JL
3. Brody CD
4. Tank DW
5. Bialek W
(2017) Collective behavior of place and non-place neurons in the hippocampal network
Neuron 96:1178–1191.
https://doi.org/10.1016/j.neuron.2017.10.027
- PubMed
- Google Scholar
1. Mézard M
2. Sakellariou J
(2011) Exact mean-field inference in asymmetric kinetic Ising systems
Journal of Statistical Mechanics 2011:L07001.
https://doi.org/10.1088/1742-5468/2011/07/L07001
- Google Scholar
(2018) Whole-Brain calcium imaging during physiological vestibular stimulation in larval zebrafish
Current Biology 28:3723–3735.
https://doi.org/10.1016/j.cub.2018.10.017
- PubMed
- Google Scholar
1. Miri A
2. Daie K
3. Arrenberg AB
4. Baier H
5. Aksay E
6. Tank DW
(2011) Spatial gradients and multidimensional dynamics in a neural integrator circuit
Nature Neuroscience 14:1150–1159.
https://doi.org/10.1038/nn.2888
- PubMed
- Google Scholar
1. Monasson R
2. Cocco S
(2011) Fast inference of interactions in assemblies of stochastic integrate-and-fire neurons from spike recordings
Journal of Computational Neuroscience 31:199–227.
https://doi.org/10.1007/s10827-010-0306-8
- PubMed
- Google Scholar
(2000) Effects of temperature on escape jetting in the squid Loligo opalescens
The Journal of Experimental Biology 203:547–557.
https://doi.org/10.1242/jeb.203.3.547
- PubMed
- Google Scholar
(2018) Maximum-entropy models reveal the excitatory and inhibitory correlation structures in cortical neuronal activity
Physical Review. E 98:012402.
https://doi.org/10.1103/PhysRevE.98.012402
- PubMed
- Google Scholar
1. Odde DJ
2. Buettner HM
(1998) Autocorrelation function and power spectrum of two-state random processes used in neurite guidance
Biophysical Journal 75:1189–1196.
https://doi.org/10.1016/S0006-3495(98)74038-X
- PubMed
- Google Scholar
1. Pandarinath C
2. O’Shea DJ
3. Collins J
4. Jozefowicz R
5. Stavisky SD
6. Kao JC
7. Trautmann EM
8. Kaufman MT
9. Ryu SI
10. Hochberg LR
11. Henderson JM
12. Shenoy KV
13. Abbott LF
14. Sussillo D
(2018) Inferring single-trial neural population dynamics using sequential auto-encoders
Nature Methods 15:805–815.
https://doi.org/10.1038/s41592-018-0109-9
- PubMed
- Google Scholar
1. Panier T
2. Romano SA
3. Olive R
4. Pietri T
5. Sumbre G
6. Candelier R
7. Debrégeas G
(2013) Fast functional imaging of multiple brain regions in intact zebrafish larvae using selective plane illumination microscopy
Frontiers in Neural Circuits 7:65.
https://doi.org/10.3389/fncir.2013.00065
- PubMed
- Google Scholar
1. Partridge LD
2. Connor JA
(1978) A mechanism for minimizing temperature effects on repetitive firing frequency
The American Journal of Physiology 234:C155–C161.
https://doi.org/10.1152/ajpcell.1978.234.5.C155
- PubMed
- Google Scholar
1. Posani L
2. Cocco S
3. Ježek K
4. Monasson R
(2017) Functional connectivity models for decoding of spatial representations from hippocampal CA1 recordings
Journal of Computational Neuroscience 43:17–33.
https://doi.org/10.1007/s10827-017-0645-9
- PubMed
- Google Scholar
(2018) Integration and multiplexing of positional and contextual information by the hippocampal network
PLOS Computational Biology 14:e1006320.
https://doi.org/10.1371/journal.pcbi.1006320
- PubMed
- Google Scholar
1. Quirin S
2. Vladimirov N
3. Yang C-T
4. Peterka DS
5. Yuste R
6. Ahrens MB
(2016) Calcium imaging of neural circuits with extended depth-of-field light-sheet microscopy
Optics Letters 41:855–858.
https://doi.org/10.1364/OL.41.000855
- PubMed
- Google Scholar
1. Ramirez AD
2. Aksay ERF
(2021) Ramp-to-threshold dynamics in a hindbrain population controls the timing of spontaneous saccades
Nature Communications 12:4145.
https://doi.org/10.1038/s41467-021-24336-w
- PubMed
- Google Scholar
1. Robertson RM
2. Money TGA
(2012) Temperature and neuronal circuit function: compensation, tuning and tolerance
Current Opinion in Neurobiology 22:724–734.
https://doi.org/10.1016/j.conb.2012.01.008
- PubMed
- Google Scholar
1. Schneider T
2. Pytte E
(1977) Random-field instability of the ferromagnetic state
Physical Review B 15:1519–1522.
https://doi.org/10.1103/PhysRevB.15.1519
- Google Scholar
(2006) Weak pairwise correlations imply strongly correlated network states in a neural population
Nature 440:1007–1012.
https://doi.org/10.1038/nature04701
- PubMed
- Google Scholar
1. Seung HS
(1996) How the brain keeps the eyes still
PNAS 93:13339–13344.
https://doi.org/10.1073/pnas.93.23.13339
- PubMed
- Google Scholar
1. Seung HS
2. Lee DD
3. Reis BY
4. Tank DW
(2000) Stability of the memory of eye position in a recurrent network of conductance-based model neurons
Neuron 26:259–271.
https://doi.org/10.1016/s0896-6273(00)81155-1
- PubMed
- Google Scholar
1. Stevenson RD
2. Josephson RK
(1990) Effects of operating frequency and temperature on mechanical power output from moth flight muscle
Journal of Experimental Biology 149:61–78.
https://doi.org/10.1242/jeb.149.1.61
- Google Scholar
1. Tang LS
2. Goeritz ML
3. Caplan JS
4. Taylor AL
5. Fisek M
6. Marder E
(2010) Precise temperature compensation of phase in a rhythmic motor pattern
PLOS Biology 8:e1000469.
https://doi.org/10.1371/journal.pbio.1000469
- PubMed
- Google Scholar
(2016) Neural assemblies revealed by inferred connectivity-based models of prefrontal cortex recordings
Journal of Computational Neuroscience 41:269–293.
https://doi.org/10.1007/s10827-016-0617-5
- PubMed
- Google Scholar
(2017) Functional coupling networks inferred from prefrontal cortex activity show experience-related effective plasticity
Network Neuroscience 1:275–301.
https://doi.org/10.1162/NETN_a_00014
- PubMed
- Google Scholar
1. Tkačik G
2. Mora T
3. Marre O
4. Amodei D
5. Palmer SE
6. Berry MJ
7. Bialek W
(2015) Thermodynamics and signatures of criticality in a network of neurons
PNAS 112:11508–11513.
https://doi.org/10.1073/pnas.1514188112
- PubMed
- Google Scholar
1. Tsodyks M
2. Sejnowski T
(1995)
Associative memory and hippocampal place cells

International Journal of Neural Systems 6:S81–S86.
- Google Scholar
1. Tubiana J
2. Wolf S
3. Panier T
4. Debregeas G
(2020) Blind deconvolution for spike inference from fluorescence recordings
Journal of Neuroscience Methods 342:108763.
https://doi.org/10.1016/j.jneumeth.2020.108763
- PubMed
- Google Scholar
Preprint
1. van der Plas TL
2. Tubiana J
3. Le Goc G
4. Migault G
5. Kunst M
6. Baier H
7. Bormuth V
8. Englitz B
9. Debrégeas G
(2021) Neural Assemblies Uncovered by Generative Modeling Explain Whole-Brain Activity Statistics and Reflect Structural Connectivity
bioRxiv.
https://doi.org/10.1101/2021.11.09.467900
- Google Scholar
1. Vladimirov N
2. Mu Y
3. Kawashima T
4. Bennett DV
5. Yang CT
6. Looger LL
7. Keller PJ
8. Freeman J
9. Ahrens MB
(2014) Light-Sheet functional imaging in fictively behaving zebrafish
Nature Methods 11:883–884.
https://doi.org/10.1038/nmeth.3040
- PubMed
- Google Scholar
1. Wang XJ
(2008) Decision making in recurrent neuronal circuits
Neuron 60:215–234.
https://doi.org/10.1016/j.neuron.2008.09.034
- PubMed
- Google Scholar
1. Wolf S
2. Supatto W
3. Debrégeas G
4. Mahou P
5. Kruglik SG
6. Sintes JM
7. Beaurepaire E
8. Candelier R
(2015) Whole-Brain functional imaging with two-photon light-sheet microscopy
Nature Methods 12:379–380.
https://doi.org/10.1038/nmeth.3371
- PubMed
- Google Scholar
1. Wolf S
2. Dubreuil AM
3. Bertoni T
4. Böhm UL
5. Bormuth V
6. Candelier R
7. Karpenko S
8. Hildebrand DGC
9. Bianco IH
10. Monasson R
11. Debrégeas G
(2017) Sensorimotor computation underlying phototaxis in zebrafish
Nature Communications 8:651.
https://doi.org/10.1038/s41467-017-00310-3
- PubMed
- Google Scholar
Software
1. Wolf S
(2023) ZF_ARTR_thermo, version swh:1:rev:fe1d2670dc3eebb1c0d810ee7add067595229d82
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:a476955c71411d6271791238fa85f14cee582def;origin=https://github.com/SebastWolf/ZF_ARTR_thermo;visit=swh:1:snp:869b0237ab4804caf34f7883d0114c339c52b732;anchor=swh:1:rev:fe1d2670dc3eebb1c0d810ee7add067595229d82
1. Zaksas D
2. Pasternak T
(2006) Directional signals in the prefrontal cortex and in area MT during a working memory for visual motion task
The Journal of Neuroscience 26:11726–11742.
https://doi.org/10.1523/JNEUROSCI.3420-06.2006
- PubMed
- Google Scholar
1. Zylberberg J
2. Strowbridge BW
(2017) Mechanisms of persistent activity in cortical circuits: possible neural substrates for working memory
Annual Review of Neuroscience 40:603–627.
https://doi.org/10.1146/annurev-neuro-070815-014006
- PubMed
- Google Scholar

Decision letter

Tatyana O Sharpee

Reviewing Editor; Salk Institute for Biological Studies, United States
Joshua I Gold

Senior Editor; University of Pennsylvania, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Emergence of time persistence in a data-driven neural network model" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Joshua Gold as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

(1) Please follow suggestions from reviewing on cross-validating model fitting

(2) Fitting the model across temperature values and discussing the relationship between the physical and model temperatures

(3) Analysis of temporal alignments between changes in the swim direction and the onset of sign change in the difference between the mean population activities of left and right hemispheres

Reviewer #1 (Recommendations for the authors):

1. Temporal alignment between neural activity and fish behavior.

Line 101: "Impact of the bath temperature on the orientational persistence in freely swimming larvae".

The authors compare the statistics of their sign(m_L-m_R) to the statistics of swim bouts in a previous paper and show that the two variables are highly correlated. However, this is just an indirect observation since they didn't establish whether changes flips in sign(m_L-m_R) are temporally aligned to changes in the fish swimming direction when analyzed simultaneously. I believe this is an important missing link in the manuscript and I suggest that the authors show that this relationship holds at least in one example – although it'd be ideal to show this at different temperature, it is sufficient to show it for one particular temperature. The dataset in Figure 4 is appropriate for this extra analysis, which should be performed both during spontaneous and stimulated epochs.

Related to this issue: what is the probability of having a change in swim direction, without a flip of sign(m_L-m_R)? And vice versa, what is the probability of having a flip of sign(m_L-m_R) without a change in swim direction?

Also related, what is the distribution of intervals between changes in the swim bout directions? If these behavioral observable is faster than the temporal resolution of the calcium indicator dynamics, this might introduce some statistical biases in the analysis.

2. Cross-validation

Line 143: "Inference of an Ising model from functional recordings of the ARTR." The authors fit an Ising model to the data showing that it captures various observables in the data. However, the model is not cross-validated. The authors should use standard cross-validation procedures to establish the goodness of fit of the various observables on held-out data. In all plots in Figure 2C-F and Appendix 2 Figure 2B-E results should be replaced with a comparison with held-out data not used to fit the model.

Moreover, what is the null hypothesis here? The authors should compare the held-out performance of their model to shuffled models trained on surrogate datasets obtained eg by destroying pairwise correlations.

3. Definition of temperature in the model

It is not explained at all how the water temperature appears in the model. Is it a missing 1/T factor in Equation (1)? This is a central issue the authors should be very clear about. What's a bit confusing is that the authors define "inverse temperature" as the effective size K, but there's no mention of the actual water temperature at all.

Line 234-236: In the Langevin formalism, the diffusion constant D usually appears in the noise autocorrelation function = 2*D*δ(t-t'). This is usually taken to be equal to the temperature (inverse β) in the partition function (which is missing from its definition!). How are the statistical mechanical temperature and the actual water temperature entering here? This is very confusing, I strongly urge the authors to clarify this point.

Reviewer #2 (Recommendations for the authors):

This is a beautiful study that combines physics with biology, treating a biological network as an inorganic system, by minimizing its free energy and maximizing entropy. It answers a lot of questions and hints at a valuable wider applicability of this approach. At the same time I am somewhat confused about the approach the authors have taken to fitting the models, explained below, and it would be good to have a conversation about that and possibly adjust some details of the approach, or to clarify. With appropriate revisions or clarifications this work would be a strong candidate for the journal.

Comments, questions, suggestions, and concerns:

1. I am puzzled why the authors fit the Ising model (and the mean field model) separately for every temperature. There are ways to include temperature directly into the model – for example, Ising models have an intrinsic dependence on (an abstracted) temperature, that is not included in equation 1 (line 157), where in typical formulations P(σ) = exp(-β*H(σ)) / Z, where β = 1/(k_b * T), so that the couplings J and biases h get naturally scaled by the temperature T, and the distribution over configurations σ becomes flatter with higher T. Of course the correspondence between the real ARTR and the real temperature, and the Ising model and the model temperature, should not be taken too literally. But qualitatively there may be an analogous effect, namely a reduction of effective couplings and a reduction in biases with increasing T, and a resulting flattening of the distribution over configurations (σ), producing a more "disordered" network with less coherence and shorter time constants. An alternative would be to treat temperature, or a function of temperature, as an input to the model, for example through the modified bias term h in the energy hamiltonian. Either of these seems more systematic than re-estimating the couplings and biases for every temperature, can the authors explain why they did it that way?

(I could not spot a systematic dependence on T of the parameters in Appendix 2 table 2 but as the authors explain, fitting the couplings and the biases can be hard to disentangle, however their comment "This owes to a progressive change in the values of both the couplings and the biases" suggests a systematic change, which I could not spot. However it could still work to incorporate a 1/T dependence explicitly.)

2. Then the traces in Figure 2A-D seem consistent with a loss of correlation amongst individual neurons in ARTR clusters, because the average becomes smaller (suggesting less coherence), but it is unclear because the Y units are "a.u.". We would need to know that the Y axis units in A is the same as that in B, and that the Y axis in C is the same as in D. Since the Y axes are labeled "a.u." so I am not sure whether this is the case. There is no reason not to provide true units (∆F/F for A,B, and average value for C,D) so that this comparison is possible.

3. Assuming my understanding of the above is correct, then, if the temperature dependence is incorporated in the model through for example the 1/T dependence or through the biases, it should be possible to fit the model at one temperature and then predict model behavior at all other temperatures. With this interpretation, the inclusion of multiple temperatures serves more as a "cross validation" of the model. I would find this a more convincing demonstration of the utility of the Ising model for network dynamics than what is currently demonstrated (i.e. more convincing than training on all temperatures and testing on all temperatures). This is analogous to the other strength of the model touted by the authors, namely that training on low-order statistics reproduces the higher-order statistics.

I realize that the temperature parameter T may not only reflect bath temperature, but also other factors like neuromodulatory inputs from the rest of the brain that depend on temperature, as the authors describe in the Discussion. Nevertheless, I think this should be attempted, even if the model T ends up scaling for example supralinearly with bath temperature, for example, T = bath-temperature + other-influences-that-are-a-function-of-T.

The same comments generalize to the mean field model. That is, I think it should be attempted to fit both the Ising and the mean field models to all temperatures and have the couplings and biases scale with 1/T (or 1/f(T) with f describing a combined dependence on bath temperature and other influences) or incorporate T into the biases as an external drive, instead of fitting the parameters separately for every temperature (of course individual fish should be fit separately). Then the held-out temperatures can serve for cross validation.

4. The way the Ising model is turned into something from which temporal dynamics can be read off is through the Monte Carlo sampler. This seems reasonable and elegant given what I understand about MC samplers and I trust it makes sense, but the paper explains the workings of the MC sampler poorly. There is a reference to some code online, but text providing the reader with a good intuition, or pseudocode, is missing. Given the importance of the MC sampler, the explanation of it, and why it can be interpreted as time evolution on long timescales, should be much more central.

5. Thus the interpretation of the Ising-ARTR correspondence seems to be first, that in both cases, less "coherence" between the units leads to longer time constants (under MC sampling in the model), and second, that increasing temperatures lead to less "coherence", thus linking temperature to time constants. If this is right, it would be useful to spell this out more clearly.

6. The Ising model is proposed to be a good model of ARTR network dynamics, and the model neurons (let's call these particles, for clarity, and call ARTR neurons, neurons) are fit to individual ARTR neurons. For readers, it would be very useful to see a graphical comparison of the dynamics of the particles, at the population level, and the dynamics of individual neurons, also at the population level. The Methods state that the mean and variance of each is fit but readers need a more personalized intuition for what the model looks like in comparison to the biological network. Although Appendix 2 Figure 1 and Appendix 2 Figures 2A(5) goes some way to depict the Ising particles and the real neurons, I think that having an elaborated population representation comparison figure in the main text would be useful, just for readers to visualize.

7. Line 165 "contralateral couplings vanish on average", also mentioned in the Discussion, seems strange because I would expect that contralateral couplings are on average negative. How does this fit in with contralateral inhibition? (As an aside, for biologists, the word vanish might be unclear and a different way of describing it might be better.) The discussion also says that the couplings are almost null but drive a subtle interplay. It seems to me that readers need to know if they are small but overall negative, or whether they can be small and overall positive; the latter would be hard to interpret but the former would be more straightforward to interpret. Why is the significance not quantified statistically?

8. The ARTR consists of excitatory and inhibitory clusters, as the authors describe, and refer to the Ising couplings, but I did not see a clear depiction or quantification of whether these couplings segregate, only a passing mention in the Discussion (lines 409-415). Can the authors depict the network connectivity accordingly, as particles with signed weights between them, and compared to their location in the brain?

9. In some ways, it is surprising that the behavior changes so much with temperature, because the field has learned from Eve Marder's group that network dynamics tend to be surprisingly robust to "crude" perturbations of the entire network. It is possible that it was difficult for evolution to compensate for the temperature dependence of ARTR, but it is also worth thinking about whether there might be an ethological purpose for the temperature dependence. Could it be useful for the animal to alternate turn direction more often when it's warm?

10. The phrasing "orientational persistence" came across as a bit of a misnomer to me, though I might be wrong – orientation typically means the angle of the animal in its environment, whereas this is about the change in angle, so orientation persistence might be interpreted by a casual reader as the fish swimming in a straight line at a fixed angle. I am not sure what else to call it though – "turn direction persistence" is probably more accurate but also clumsier. Is there a better term? If you decide to keep using "orientational persistence" it would be useful to very carefully and explicitly explain what you mean by it.

11. Cold-blooded animals have a body temperature that quickly converges to that of their surroundings (especially tiny fish). Their neural populations will likely experience a wider range of temperature fluctuation and therefore follow certain thermal dynamics laws unlike any warm-blooded animals (such that the flipping rate increases with temperature following Arrhenius law). It is important to discuss if the time persistence that emerges from the network model is a special case for cold-blooded animals experiencing different temperatures (through the 1/(kB*T) scaling of the couplings and biases), a direct thermal effect, rather than "purposeful" neural network computation (through external inputs that could be modeled by making the biases depend on T).

12. It is therefore good to discuss how applicable this free energy landscape description of a neural network is to other neural networks with persistent activities in other animal models, with different sensory inputs, in order to reach a larger audience and have wider potential applications. What types of time persistence can this interpretation address and what are its limitations?

13. Line 9: "networks" should probably be "network".

14. Figure 1 C and G: labeling the colors by their temperature in the figure would be helpful.

15. Line 246, it is Figure 3B showing the energy landscape, not 3C.

https://doi.org/10.7554/eLife.79541.sa1

Author response

Essential revisions:

(1) Please follow suggestions from reviewing on cross-validating model fitting

(2) Fitting the model across temperature values and discussing the relationship between the physical and model temperatures

(3) Analysis of temporal alignments between changes in the swim direction and the onset of sign change in the difference between the mean population activities of left and right hemispheres

We have followed these recommendations to write our amended manuscript. In particular, we now:

(1) include a systematic cross-validation of the inferred models following the comments by Reviewer 1. We show the predictions of the model (for observables such as the mean activity and the pairwise correlations, and for log-likelihoods) in the figures on the main text, e.g. figure 3A, as well as in the appendix figure 3.

(2) explain why we have fitted a model for each water temperature and each fish, and carry out a systematic comparison of the models inferred at different water temperatures for the same fish. We have added a new figure 3B to show the absence of correlation between models corresponding to different water temperatures. We have also included a new section in the Discussion, entitled “Origin and functional significance of the temperature dependence of the ARTR dynamics” to better discuss the role of temperature on the activity. We have also slightly rewritten the introduction of the model in the Results section to better distinguish the notions of water temperature and model temperature (implicitly set by the amplitude of the J,h parameters), as our first manuscript could be confusing on this point.

(3) point out the exact figure in Dunn et al., eLife 2016 in which the authors demonstrate that the sign of the difference in activity of the right and left ARTR populations dictates the swim direction. Replicating these observations would thus be redundant but would also require a different and more complex experimental setup. Indeed, they were obtained by recording fictive turns using electrical recordings from peripheral motor nerves in paralysed larvae, as actual tail bouts tend to occur at very low frequency in tethered configurations.

Reviewer #1 (Recommendations for the authors):

1. Temporal alignment between neural activity and fish behavior.

Line 101: "Impact of the bath temperature on the orientational persistence in freely swimming larvae".

The authors compare the statistics of their sign(m_L-m_R) to the statistics of swim bouts in a previous paper and show that the two variables are highly correlated. However, this is just an indirect observation since they didn't establish whether changes flips in sign(m_L-m_R) are temporally aligned to changes in the fish swimming direction when analyzed simultaneously. I believe this is an important missing link in the manuscript and I suggest that the authors show that this relationship holds at least in one example – although it'd be ideal to show this at different temperature, it is sufficient to show it for one particular temperature. The dataset in Figure 4 is appropriate for this extra analysis, which should be performed both during spontaneous and stimulated epochs.

Related to this issue: what is the probability of having a change in swim direction, without a flip of sign(m_L-m_R)? And vice versa, what is the probability of having a flip of sign(m_L-m_R) without a change in swim direction?

We agree with the reviewer that the relation between the change in swim turn direction and the neuronal switch between the left and right ARTR subpopulations, is central to the proposed scenario. This relationship has been previously established by Dunn et al. in their eLife article (Dunn 2016). Below we reproduce the panel E of figure 5 —figure supplement 2 from this article, which precisely demonstrates the robust correlation that exists between m_L-m_R and the turning direction. Notice that to obtain these results, the authors recorded fictive turns using electrical recordings from peripheral motor nerves on both sides of the tail of paralysed larvae, as actual tail bouts tend to occur at very low frequency in tethered configurations. In 2017, we also indirectly confirmed this relationship by establishing that m_L-m_R was a strong predictor of the gaze direction, and that the latter was strongly correlated with the tail bout orientations (Wolf et al. 2017).

We modified the manuscript to clarify this point by explicitly referring to this figure:

line 111-114

“It has previously been shown that the ARTR governs the selection of swim bout orientations: turn bouts are preferentially executed in the direction of the most active (right or left) ARTR subcircuit (Dunn et al., 2016, Wolf et al., 2017), such that sign $(m_{L} (t) - m_{R} (t))$ constitutes a robust predictor of the turning direction of the animal, see figure 5 —figure supplement 2E in (Dunn et al., 2016).”

Also related, what is the distribution of intervals between changes in the swim bout directions? If these behavioral observable is faster than the temporal resolution of the calcium indicator dynamics, this might introduce some statistical biases in the analysis.

The distributions of orientational persistence times, as measured experimentally, are shown in Author response image 1. Since two separate turning swim bouts are needed to detect a change in turn direction, this distribution shows a rapid decay for T<1s, which is the mean interbout interval. The left-right alternation process thus occurs on time scales that are particularly well adapted to calcium imaging, whose temporal resolution after spike inference is of the order of a few hundreds of ms.

Author response image 1

Download asset Open asset

2. Cross-validation

Line 143: "Inference of an Ising model from functional recordings of the ARTR." The authors fit an Ising model to the data showing that it captures various observables in the data. However, the model is not cross-validated. The authors should use standard cross-validation procedures to establish the goodness of fit of the various observables on held-out data. In all plots in Figure 2C-F and Appendix 2 Figure 2B-E results should be replaced with a comparison with held-out data not used to fit the model.

Moreover, what is the null hypothesis here? The authors should compare the held-out performance of their model to shuffled models trained on surrogate datasets obtained eg by destroying pairwise correlations.

In order to answer this important point, we include in the revised version of the paper a systematic cross-validation of the inferred models summarized in Appendix 2 Figure 3.

Indeed for each experiment, we divided the data sets in two parts: 75% of each data set is used as a training set and the remaining 25% is used as a test set. Each training set is used to infer an Ising model. We then compare the mean activity and covariance of the test set with the one computed from the simulated data generated by the models (see Appendix 2 Figure 3A-B). We also computed the relative variation of the models' log likelihood computed on the training data and the test data (see Appendix 2 Figure 3C). The very weak value of this relative variation confirms the quality of the fits.

Moreover, as suggested by the reviewer, we decided to use the independent model as a null hypothesis. We thus compared the Ising models fitted on the data with the independent model. The independent model only depends on the firing rates of the neurons. We then compare the mean activity and covariance of the test set with the one computed from the simulated data generated by the independent models (see Appendix 2 Figure 3E-G). As expected, these independent models with no connections poorly reproduce the data. The excess log likelihood confirms this result (Appendix 2 Figure 3G). Indeed, the relative variation, between the Ising and the independent models, of the log likelihood computed on the training data and the test data is around 50%.

In addition, we also followed the reviewer’s suggestion and computed the KL divergence in Figure 3 using the test sets. Indeed, the distribution of the KL divergence between the experimental test datasets and their associated Ising models is much smaller than those obtained between experimental test datasets and Ising models trained on different recordings (red distribution). We also added, for comparison, the distribution of KL divergence obtained using independent models, which demonstrate the importance of interneuronal connections to reproduce the data.

line 206–209

“This agreement crucially relies on the presence of inter-neuronal couplings in order to reproduce the pairwise correlations in the activity: a model with no connection (i.e. the independent model, see Methods) fitted to reproduce the neural firing rates, offers a very poor description of the data, see Figure 3A (dark blue distribution) and Appendix 2 Figure 3E-G.”

3. Definition of temperature in the model

It is not explained at all how the water temperature appears in the model. Is it a missing 1/T factor in Equation (1)? This is a central issue the authors should be very clear about. What's a bit confusing is that the authors define "inverse temperature" as the effective size K, but there's no mention of the actual water temperature at all.

Line 234-236: In the Langevin formalism, the diffusion constant D usually appears in the noise autocorrelation function = 2*D*δ(t-t'). This is usually taken to be equal to the temperature (inverse β) in the partition function (which is missing from its definition!). How are the statistical mechanical temperature and the actual water temperature entering here? This is very confusing, I strongly urge the authors to clarify this point.

We agree that the previous version of the manuscript was not clear enough in the definition and introduction of temperature(s). In the new version, we make clear that temperature refers to the temperature of water, and denote this quantity by T.

Regarding the Ising distribution defined in Equation (1), we implicitly set the model temperature T_model to 1, which amounts to omitting it in Equation (1). Considering a different value for T_model would simply amount to rescaling the biases and couplings parameters by the same quantity, and would not affect the model distribution.

The fact that the model temperature T_model is equal to unity is consistent with the distribution of left and right activities in Equation (3) and the Langevin dynamics in Equation (4), see covariance of the noise epsilon, which is 2* δ(t-t’). As stressed by the referee, the covariance should be equal to 2 * T_model * δ(t-t’), hence our expression corresponds to T_model=1.

We now clearly state after Equation 1 that the model temperature is set to one throughout the paper.

line 172–174

“Notice that in Equation 1, the energy term in the parenthesis is not scaled by a thermal energy as in the Maxwell–Boltzmann statistics. We thus implicitly fix the model temperature to unity; of course, this model temperature has no relation with the water temperature T.”

The use of the vocable inverse temperature (which was done with quotation marks) in Appendix right after Equation (18) to design the sizes K_L, K_R was very confusing, and we apologize for that. We simply meant that, as the sizes K enter multiplicatively the free-energy F they could be formally interpreted as effective inverse temperatures in the Boltzmann factor e^-F. We have removed this confusing part of the sentence, and hope that the current formulation is more clear.

Reviewer #2 (Recommendations for the authors):

This is a beautiful study that combines physics with biology, treating a biological network as an inorganic system, by minimizing its free energy and maximizing entropy. It answers a lot of questions and hints at a valuable wider applicability of this approach. At the same time I am somewhat confused about the approach the authors have taken to fitting the models, explained below, and it would be good to have a conversation about that and possibly adjust some details of the approach, or to clarify. With appropriate revisions or clarifications this work would be a strong candidate for the journal.

Comments, questions, suggestions, and concerns:

1. I am puzzled why the authors fit the Ising model (and the mean field model) separately for every temperature. There are ways to include temperature directly into the model – for example, Ising models have an intrinsic dependence on (an abstracted) temperature, that is not included in equation 1 (line 157), where in typical formulations P(σ) = exp(-β*H(σ)) / Z, where β = 1/(k_b * T), so that the couplings J and biases h get naturally scaled by the temperature T, and the distribution over configurations σ becomes flatter with higher T. Of course the correspondence between the real ARTR and the real temperature, and the Ising model and the model temperature, should not be taken too literally. But qualitatively there may be an analogous effect, namely a reduction of effective couplings and a reduction in biases with increasing T, and a resulting flattening of the distribution over configurations (σ), producing a more "disordered" network with less coherence and shorter time constants. An alternative would be to treat temperature, or a function of temperature, as an input to the model, for example through the modified bias term h in the energy hamiltonian. Either of these seems more systematic than re-estimating the couplings and biases for every temperature, can the authors explain why they did it that way?

(I could not spot a systematic dependence on T of the parameters in Appendix 2 table 2 but as the authors explain, fitting the couplings and the biases can be hard to disentangle, however their comment "This owes to a progressive change in the values of both the couplings and the biases" suggests a systematic change, which I could not spot. However it could still work to incorporate a 1/T dependence explicitly.)

In the new manuscript we have made clearer that β is a rescaling parameter, which is integrated in the inferred value of couplings, see answer to Major point 3 of Reviewer 1.

line 172–174

In addition, we now show in Appendix 2 Figure 6 some scatter plots of the couplings and biases at different temperatures for one fish, as well as statistics for all fish. These results clearly indicate the absence of any strong correlation between the Ising parameters at different temperatures (for the same animal). As a consequence, there is no simple global rescaling between the Ising parameters attached to two water temperatures. This is why we have fitted a model for each data set.

line 228–233

“We next examined the dependency of the Ising model parameters on the water temperature. To do so, for each fish, we selected two different water temperatures, and the corresponding sets of inferred biases and couplings, {h_i, J_ij}. We then computed the Pearson correlation coefficient R² of the biases and of the coupling matrices at these two temperatures (inset of Appendix 2 Figure 6). We saw no clear correlation between the model parameters at different temperatures, as shown by the distribution of R² computed across fish and across every temperatures (Appendix 2 Figure 6).”

2. Then the traces in Figure 2A-D seem consistent with a loss of correlation amongst individual neurons in ARTR clusters, because the average becomes smaller (suggesting less coherence), but it is unclear because the Y units are "a.u.". We would need to know that the Y axis units in A is the same as that in B, and that the Y axis in C is the same as in D. Since the Y axes are labeled "a.u." so I am not sure whether this is the case. There is no reason not to provide true units (∆F/F for A,B, and average value for C,D) so that this comparison is possible.

The notation a.u. was misleading, as we are showing the time trace of the mean (over each region) of the binarized activity. The axis labels have been corrected in Figure 2 (beware that the numbers of the panels have changed with respect to the first version of the manuscript).

3. Assuming my understanding of the above is correct, then, if the temperature dependence is incorporated in the model through for example the 1/T dependence or through the biases, it should be possible to fit the model at one temperature and then predict model behavior at all other temperatures. With this interpretation, the inclusion of multiple temperatures serves more as a "cross validation" of the model. I would find this a more convincing demonstration of the utility of the Ising model for network dynamics than what is currently demonstrated (i.e. more convincing than training on all temperatures and testing on all temperatures). This is analogous to the other strength of the model touted by the authors, namely that training on low-order statistics reproduces the higher-order statistics.

In the new manuscript, we include in the revised version of the paper a systematic cross-validation of the inferred models summarized in Appendix 2 Figure 3, see Point 2 in the Answer to Reviewer 1. Note that this cross validation is done at each temperature separately. We discuss the putative relationship between models at different temperatures below.

line 644–658

“We cross-validated the Ising models (see Appendix Figure 3) dividing the data sets in two parts: for each experiment, 75% of each data set is used as a training set and the remaining 25% is used as a test set. Each training set is used to infer an Ising model. We then compare the mean activity and covariance of the test set with the one computed from the simulated data generated by the models…” (see manuscript)

I realize that the temperature parameter T may not only reflect bath temperature, but also other factors like neuromodulatory inputs from the rest of the brain that depend on temperature, as the authors describe in the Discussion. Nevertheless, I think this should be attempted, even if the model T ends up scaling for example supralinearly with bath temperature, for example, T = bath-temperature + other-influences-that-are-a-function-of-T.

The same comments generalize to the mean field model. That is, I think it should be attempted to fit both the Ising and the mean field models to all temperatures and have the couplings and biases scale with 1/T (or 1/f(T) with f describing a combined dependence on bath temperature and other influences) or incorporate T into the biases as an external drive, instead of fitting the parameters separately for every temperature (of course individual fish should be fit separately). Then the held-out temperatures can serve for cross validation.

Varying the water temperature was done here to investigate how the ARTR circuit persistent dynamics could be modeled across different behavioral conditions. The use of Ising is obviously not limited to this particular experimental protocol, which is of course adequate for cold-blooded animals only.

In general, we do not expect that the Ising models corresponding to different conditions can be related to each other in a simple way, e.g. through a global rescaling of the biases and/or of the couplings. However, we agree with Reviewer 2 that it is legitimate to attempt at relating the models inferred from the recordings at different water temperatures. We have done so in two ways:

– We have directly compared the couplings and biased inferred from the recordings of activity at different water temperatures for the same fish. Results are shown in Appendix 2 Figure 6, and show no clear interrelation between the parameters corresponding to different temperatures, as would be the case if parameters would transform according to a simple global rescaling.

– To confirm this finding, following the recommendations of Reviewer 2, we have carried out the following computation. We first infer the parameters h_i , J_ij of an Ising model from the recording of a fish at water temperature T. Then we rescale these parameters as a x h_i , b x J_ij and we generate synthetic activity with these rescaled ising models. We then compute the KL divergence between the distribution of activity defined by this rescaled Ising model and the recorded data at another temperature T’ (for the same fish). We then minimize this KL divergence over the global rescaling factors a and b in the [0.5;1.5] ranges, and compare this lowest divergence to the one between the same recorded data (at temperature T’) and the Ising model inferred from those data (as shown in Figure 3A of the main text). Ratios of these KL divergences are reported in Author response image 2. These results show that the global rescalings of the parameters inferred from data at temperature T cannot reproduce correctly the data at another temperature T’.

Author response image 2

Download asset Open asset

4. The way the Ising model is turned into something from which temporal dynamics can be read off is through the Monte Carlo sampler. This seems reasonable and elegant given what I understand about MC samplers and I trust it makes sense, but the paper explains the workings of the MC sampler poorly. There is a reference to some code online, but text providing the reader with a good intuition, or pseudocode, is missing. Given the importance of the MC sampler, the explanation of it, and why it can be interpreted as time evolution on long timescales, should be much more central.

We give now the details of the MC sampler in the main text and in methods of the revised version.

line 179-186

“The algorithm starts from a random configuration of activity, then picks up uniformly at random a neuron index, say, i. The activity s_i of neuron i is then stochastically updated to 0 or to 1, with probabilities that depend on the current states s_i of the other neurons (see Equation 8} in Methods, and code provided). The sampling procedure is iterated, ensuring convergence towards the distribution P in Equation 1. This in silico MC dynamics is not supposed to reproduce any realistic neural dynamics, except for the locality in the activity configuration s space.”

line 634-643

“In order to generate synthetic activity, we resorted to Gibbs sampling, a class of Monte Carlo Markov Chain method, also known as Glauber dynamics…” (see manuscript)

5. Thus the interpretation of the Ising-ARTR correspondence seems to be first, that in both cases, less "coherence" between the units leads to longer time constants (under MC sampling in the model), and second, that increasing temperatures lead to less "coherence", thus linking temperature to time constants. If this is right, it would be useful to spell this out more clearly.

The origin of the dependence of the persistence time on water temperature is rather subtle. We think that this dependence is not directly related to a change of the correlation structure in the activity. To support this point we provide a new panel (Appendix 2 figure 2C) showing that no significant relation between the average correlations and the temperature exists.

However, as our mean-field analysis shows, barriers increase with temperature, see Figure 4E, a fact that may appear counterintuitive at first sight. As a consequence, at high temperature, only the low-low activity state is accessible in practice to the system, and the mean activity remains low, see Appendix 2 Figure 2B, with fluctuations within the low-low state. Conversely, at low water temperatures, barriers separating the low-low and the high-low (or low-high) states are weaker, so the active states are accessible. This has two consequences. First, the mean activity is higher at low temperature (Appendix 2 Figure 2B). Second, the system may remain trapped for some time in such an active state before switching to the other side, and reversing activity. This is the origin of the longer persistence at low temperature.

We have improved the manuscript to make these points clearer, see end of section “Barriers in the free-energy landscape and dynamical paths between states”.

line 306-313

“As a consequence, at high temperature, only the low-low activity state is accessible in practice to the system, and the mean activity remains low, see Appendix 2 Figure 2D, with fluctuations within the low-low state. Conversely, at low water temperatures, barriers separating the low-low and the active high-low or low-high states are weaker, so the latter become accessible. As a first consequence, the mean activity is higher at low temperature (Appendix 2 Figure 2D). Furthermore, the system remains trapped for some time in such an active state before switching to the other side, e.g. from high-low to low-high. This is the origin of the longer persistence time observed at low temperature.”

6. The Ising model is proposed to be a good model of ARTR network dynamics, and the model neurons (let's call these particles, for clarity, and call ARTR neurons, neurons) are fit to individual ARTR neurons. For readers, it would be very useful to see a graphical comparison of the dynamics of the particles, at the population level, and the dynamics of individual neurons, also at the population level. The Methods state that the mean and variance of each is fit but readers need a more personalized intuition for what the model looks like in comparison to the biological network. Although Appendix 2 Figure 1 and Appendix 2 Figures 2A(5) goes some way to depict the Ising particles and the real neurons, I think that having an elaborated population representation comparison figure in the main text would be useful, just for readers to visualize.

We agree with Reviewer 2 that such a representation is useful. We now show in Figure 2 of the main text raster plots of the activity in the recording (panel A) and in the simulated Ising dynamics (panel C). We also provide in the supplementary materials the same comparison for every recording.

7. Line 165 "contralateral couplings vanish on average", also mentioned in the Discussion, seems strange because I would expect that contralateral couplings are on average negative. How does this fit in with contralateral inhibition? (As an aside, for biologists, the word vanish might be unclear and a different way of describing it might be better.) The discussion also says that the couplings are almost null but drive a subtle interplay. It seems to me that readers need to know if they are small but overall negative, or whether they can be small and overall positive; the latter would be hard to interpret but the former would be more straightforward to interpret. Why is the significance not quantified statistically?

We now provide a detailed analysis on both contralateral correlation and couplings in the Appendix 2 Figure 4. Indeed we show that strong negative correlation and couplings systematically correspond to pairs of contralateral neurons. We provide in the main text a paragraph on this point.

line 218-222

“… computing the fraction of neuronal pairs (i,j) that are contralateral for each value of the coupling J_ij or the Pearson correlation (Appendix 2 Figure 4D-E). Large negative values of couplings or correlations systematically correspond to contralateral pairs of neurons, whereas large positive values correspond to ipsilateral pairs of neurons.”

8. The ARTR consists of excitatory and inhibitory clusters, as the authors describe, and refer to the Ising couplings, but I did not see a clear depiction or quantification of whether these couplings segregate, only a passing mention in the Discussion (lines 409-415). Can the authors depict the network connectivity accordingly, as particles with signed weights between them, and compared to their location in the brain?

We do not see a spatial organization in the inferred couplings that would reflect the neuronal type organization. We stress that the relationship between the inferred couplings J_ij and the true physiological interactions is subtle. For instance, in our model, couplings are symmetric, which makes the distinction between inhibitory and excitatory neurons uneasy. In addition, the sampling time bin is long compared to synaptic integrations times, which makes causal dependencies hard to infer even with models where connections are not a priori required to be symmetric, see for instance Figure 7b and related discussion in a previous publication by two of us (SC,RM): J Comput Neurosci DOI 10.1007/s10827-010-0306-8

9. In some ways, it is surprising that the behavior changes so much with temperature, because the field has learned from Eve Marder's group that network dynamics tend to be surprisingly robust to "crude" perturbations of the entire network. It is possible that it was difficult for evolution to compensate for the temperature dependence of ARTR, but it is also worth thinking about whether there might be an ethological purpose for the temperature dependence. Could it be useful for the animal to alternate turn direction more often when it's warm?

The reviewer is right that some networks have evolved to maintain robust dynamics across a large thermal range. However, this is rather the exception than the rule, and even the frequency of the pyloric rhythm of the crab, studied by Marder, increases with temperature (although the phase is maintained).

Regarding whether or not the thermal dependence in ARTR dynamics is favorable for the animal, it is hard to conclude as many other kinematic features are also modified, as we have quantified recently. These different points are now discussed in a new section of the discussion.

line 379-400

“The brains of cold-blooded animals need to operate within the range of temperature that they experience in their natural habitat, e.g. 15-33°C for zebrafish [Gau et al., 2013]. This is a peculiarly stringent requirement since most biophysical processes are dependent on the temperature. In some rare instances, regulation mechanisms might stabilize the circuit dynamics in order to preserve its function, as best exemplified by the pyloric rhythm of the crab whose characteristic phase relationship is maintained over an extended temperature range [Tang et al., 2010]. Yet in general, an increase in temperature tends to increase the frequency of oscillatory processes [Robertson et al., 2012]. The observed acceleration of the ARTR left/right alternation with increasing temperature, could thus directly result from temperature-dependent cellular mechanisms. Furthermore, one cannot rule out the possibility that the ARTR dynamics could also be indirectly modulated by temperature via thermal-dependent descending neuromodulatory inputs.

As a result of this thermal modulation of the neuronal dynamics, many cold-blooded animals also exhibit temperature-dependence of their behavior [Long et al., 2008, Neumeister et al., 2000, Stevenson et al., 1990]. Here we were able to quantitatively relate the two processes (neuronal and motor) by demonstrating that an increase in temperature consistently alters the pattern of spontaneous navigation by increasing the left/right alternation frequency. Interpreting the functional relevance of this modification of the swimming pattern is tricky, since many other features of the animal’s navigation are concurrently impacted by a change in temperature, such as the bout frequency, turning rate, turn amplitude, etc. Nevertheless, we were able to show in a recent study that this thermal dependence of the swimming kinematic endows the larva with basic thermophobic capacity, thus efficiently protecting them from exposure to the hottest regions of their environment [Le Goc et al., 2021].”

10. The phrasing "orientational persistence" came across as a bit of a misnomer to me, though I might be wrong – orientation typically means the angle of the animal in its environment, whereas this is about the change in angle, so orientation persistence might be interpreted by a casual reader as the fish swimming in a straight line at a fixed angle. I am not sure what else to call it though – "turn direction persistence" is probably more accurate but also clumsier. Is there a better term? If you decide to keep using "orientational persistence" it would be useful to very carefully and explicitly explain what you mean by it.

We thank the reviewer for this suggestion: turn direction persistence is indeed better. We modified the text accordingly.

11. Cold-blooded animals have a body temperature that quickly converges to that of their surroundings (especially tiny fish). Their neural populations will likely experience a wider range of temperature fluctuation and therefore follow certain thermal dynamics laws unlike any warm-blooded animals (such that the flipping rate increases with temperature following Arrhenius law). It is important to discuss if the time persistence that emerges from the network model is a special case for cold-blooded animals experiencing different temperatures (through the 1/(kB*T) scaling of the couplings and biases), a direct thermal effect, rather than "purposeful" neural network computation (through external inputs that could be modeled by making the biases depend on T).

The observed speeding up of the ARTR alternating dynamics with the water temperature T is consistent with a direct thermal effect. We agree that the biases could also exhibit some variation with T. This is not unexpected, as different membrane currents may have different temperature dependence, see our response to point 9. However we do not see a strong dependence of the biases with T, see Author response image 3 showing the mean value of h_i as a function of the water temperature.

Author response image 3

Download asset Open asset

12. It is therefore good to discuss how applicable this free energy landscape description of a neural network is to other neural networks with persistent activities in other animal models, with different sensory inputs, in order to reach a larger audience and have wider potential applications. What types of time persistence can this interpretation address and what are its limitations?

Time persistence emerges in our work through the sampling of (generally) two metastable states. However, our approach could be applied to more complex landscapes, with more metastable states. Examples of such situations relevant for computation can be found in several works, e.g. Harvey et al. (2012), Mazzucato et al. (2019), Brinkman et al. (2022). It would be interesting to test our approach in the case of continuous attractors, e.g. ring attractors that serve as a theoretical basis for head direction networks (and their higher dimensional extension for place/grid cells circuits). Metastability in this context would be different as a continuum of states are expected to coexist. However we expect motion along the attractor direction, in particular drift, to happen on longer time scales than the damping of fluctuations along transverse directions. Some of us have already inferred graphical models on multi-electrode recordings of hippocampal place cells, see for instance Posani et al. (2018) PLoS Comp Bio 14:e1006320. However the limited number of recorded cells (compared to the recordings reported in our present manuscript) did not allow us to obtain sufficiently accurate generative performance with MC dynamics.

13. Line 9: "networks" should probably be "network".

This has been corrected.

14. Figure 1 C and G: labeling the colors by their temperature in the figure would be helpful.

The figure has been modified.

15. Line246, it is Figure 3B showing the energy landscape, not 3C.

This has been corrected.

https://doi.org/10.7554/eLife.79541.sa2

Article and author information

Author details

Sebastien Wolf

Laboratory of Physics of the Ecole Normale Supérieure, CNRS UMR 8023 & PSL Research, Sorbonne Université, Université de Paris, Paris, France

Contribution
Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - review and editing

Contributed equally with
Guillaume Le Goc

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-9394-3291
Guillaume Le Goc

Institut de Biologie Paris-Seine (IBPS), Laboratoire Jean Perrin, Sorbonne Université, CNRS, Paris, France

Contribution
Data curation, Validation, Investigation, Visualization, Methodology, Writing - review and editing

Contributed equally with
Sebastien Wolf

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-6946-1142
Georges Debrégeas

Institut de Biologie Paris-Seine (IBPS), Laboratoire Jean Perrin, Sorbonne Université, CNRS, Paris, France

Contribution
Conceptualization, Funding acquisition, Investigation, Methodology, Writing – original draft, Project administration

Contributed equally with
Simona Cocco and Rémi Monasson

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-3698-4497
Simona Cocco

Laboratory of Physics of the Ecole Normale Supérieure, CNRS UMR 8023 & PSL Research, Sorbonne Université, Université de Paris, Paris, France

Contribution
Conceptualization, Supervision, Methodology, Writing – original draft

Contributed equally with
Georges Debrégeas and Rémi Monasson

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-1852-7789
Rémi Monasson

Laboratory of Physics of the Ecole Normale Supérieure, CNRS UMR 8023 & PSL Research, Sorbonne Université, Université de Paris, Paris, France

Contribution
Conceptualization, Formal analysis, Supervision, Investigation, Methodology, Writing – original draft, Project administration

Contributed equally with
Georges Debrégeas and Simona Cocco

For correspondence
remi.monasson@phys.ens.fr

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4459-0204

Funding

Agence Nationale de la Recherche (Locomat)

Rémi Monasson

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

SW acknowledges support by a fellowship from the Fondation pour la Recherche Médicale (SPF 201809007064), and GLG by the Systems Biology network of Sorbonne Université. We thank the IBPS fish facility staff for the fish maintenance. We are grateful to Carounagarane Dore for his contribution to the design of the experimental setups. We thank Misha Ahrens for providing the GCaMP line.

Ethics

All experiments were approved by Le Comité d'Éthique pour l'Expérimentation Animale Charles Darwin (02601.01).

Senior Editor

Joshua I Gold, University of Pennsylvania, United States

Reviewing Editor

Tatyana O Sharpee, Salk Institute for Biological Studies, United States

Version history

Preprint posted: February 4, 2022 (view preprint)
Received: April 17, 2022
Accepted: March 13, 2023
Accepted Manuscript published: March 14, 2023 (version 1)
Version of Record published: May 10, 2023 (version 2)

Copyright

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