 Research article
 Open Access
 Published:
Thermal modulation of Zebrafish exploratory statistics reveals constraints on individual behavioral variability
BMC Biology volume 19, Article number: 208 (2021)
Abstract
Background
Variability is a hallmark of animal behavior. It contributes to survival by endowing individuals and populations with the capacity to adapt to everchanging environmental conditions. Intraindividual variability is thought to reflect both endogenous and exogenous modulations of the neural dynamics of the central nervous system. However, how variability is internally regulated and modulated by external cues remains elusive. Here, we address this question by analyzing the statistics of spontaneous exploration of freely swimming zebrafish larvae and by probing how these locomotor patterns are impacted when changing the water temperatures within an ethologically relevant range.
Results
We show that, for this simple animal model, five shortterm kinematic parameters — interbout interval, turn amplitude, travelled distance, turn probability, and orientational flipping rate — together control the longterm exploratory dynamics. We establish that the bath temperature consistently impacts the means of these parameters, but leave their pairwise covariance unchanged. These results indicate that the temperature merely controls the sampling statistics within a welldefined kinematic space delineated by this robust statistical structure. At a given temperature, individual animals explore the behavioral space over a timescale of tens of minutes, suggestive of a slow internal state modulation that could be externally biased through the bath temperature. By combining these various observations into a minimal stochastic model of navigation, we show that this thermal modulation of locomotor kinematics results in a thermophobic behavior, complementing direct gradientsensing mechanisms.
Conclusions
This study establishes the existence of a welldefined locomotor space accessible to zebrafish larvae during spontaneous exploration, and quantifies selfgenerated modulation of locomotor patterns. Intraindividual variability reflects a slow diffusivelike probing of this space by the animal. The bath temperature in turn restricts the sampling statistics to subregions, endowing the animal with basic thermophobicity. This study suggests that in zebrafish, as well as in other ectothermic animals, ambient temperature could be used to efficiently manipulate internal states in a simple and ethological way.
Background
Variability, both inter and intraindividual, is an ubiquitous trait of animal behavior [1]. Intraindividual variability may participate in efficient strategies, as best exemplified by the alternation of exploration and exploitation phases during foraging [2]. It can also endow the animal, or the population, with robustness, i.e., the ability to rapidly and efficiently cope with changing environmental conditions [3, 4]. The idea, known as bethedging, is that a modest loss in fitness associated with phenotypic variability could be balanced by the gain in leniency when facing unexpected and possibly hostile conditions. The origin of interindividual variability may be attributed to genetic, epigenetic, or developmental differences. Intraindividual variability may in turn reflects spontaneous transitions between distinct brain states, i.e., patterns of persistent neural activity [5, 6]. It may also be the signature of endogenous modulations in the production of neuromodulators [7].
Although the functional significance of variability in animal behavior is now largely recognized [8], the way it is regulated and modulated by external cues, as well as its neuronal substrate remain elusive. To address this question, one not only needs to quantify variability, but also manipulate it in a physiologically relevant manner, in an animal that is accessible to both behavioral and neuronal circuit interrogation. Here, we used the zebrafish larva as a model vertebrate as it is uniquely amenable to in vivo whole brain functional imaging [9–11] and to highthroughput behavioral studies [12, 13].
As an ectothermic animal, zebrafish must actively navigate towards regions of its environment that are thermally optimal for its thriving [14], while potentially being exposed to a wide range of temperatures [15]. How fish swim in thermal gradients has been extensively studied [16], and the neuronal circuits underlying this thermotactic process have been identified [17]. Zebrafish larvae integrate thermal signals (change in temperature) over a subsecond time window, and adapt their forthcoming movement accordingly in order to eventually move towards optimal zones.
Here, we focus on the exploratory dynamics at various but spatially uniform temperatures. We use a reductive approach, as previously introduced [18], to quantify its spontaneous locomotion using a finite number of shortterm kinematic parameters. We then quantify how the bath temperature not only impacts the mean of these parameters, but also their statistical distribution (variability) and pairwise covariance. We further assess the timescale over which this behavioral variability unfolds at the level of individual animals. From this detailed analysis, we build a numerical model of zebrafish larvae navigation at all temperatures over the physiologically relevant range. Finally, we use this model to demonstrate how this thermal adaptation of spontaneous swimming pattern may complement the thermotactic mechanism, based on direct gradient sensing, in order for the animal to limit its presence in potentially harmful environments.
Results
A behavioral assay to record spontaneous navigation at different temperatures
Batches of 10 zebrafish larvae aged 57 days were videomonitored at 25 frames/second for periods of 30 min as they freely swam in a rectangular 100 × 45 × 4.5 mm^{3} pool at a constant and uniform temperature (Fig. 1A, see “Methods” section). For each batch, we successively imposed up to 5 values of temperature (18, 22, 26, 30, and 33 ^{∘}C) in a random order. This thermal range spans the nonlethal conditions for larval zebrafish and has been shown to be effectively encountered by the animal in its natural habitat [19]. Each 30minlong recording session was preceded by a 14minlong period of habituation to allow the animals to reach their steadystate exploratory regime. A total of 10 batches per temperature involving 170 different fish were used.
Larval zebrafish swim in discrete bouts lasting for about 100 ms, interspersed with ∼ 1–2 periods of rest. As we aim to probe how the bath temperature impacts the longterm exploratory process, we focus on the characterization of a few kinematic parameters associated with each bout. We thus ignore the fine structure of the swimming events, such as the amplitude of the tail deflection or beating frequency [20, 21] but examine their resulting heading reorientation and linear displacements. The center of mass coordinates and orientation of each larva in every frame were extracted using FastTrack [22] (see “Methods” section). For each identified swim bout, we computed three scalar parameters (Fig. 2A) whose statistics control the fish spatiotemporal exploration [18]: (i) the interbout interval (IBI), δt_{n}, is the idle time following the bout event; (ii) the displacement, d_{n}, is the traveled distance associated with the bout; and (iii) the reorientation angle, δθ_{n}, denotes the change in heading direction.
Zebrafish larvae, as most animals, tend to adopt a distinct navigational behavior at proximity of boundaries [23]. Tracking was thus performed within the innermost region of the arena, at a minimum distance of 5 mm from the walls. As a result, individual fish were not tracked continuously over the entire recording periods, but along trajectories (from one wall to another). In the analysis, we ignored trajectories that last less than 25 s. Example trajectories for three temperatures are shown in Fig. 1C, where each dot indicates the location of a swim bout, while its size reflects the interbout interval. This comparison provides a first qualitative illustration of the effect of temperature on the fish exploration. At low temperatures (18 ^{∘}C), the trajectories are relatively straight, comprising a majority of small discrete forward bouts executed at relatively low frequency. At high temperatures, the trajectories appear much more meandering, with more frequent and ample reorienting maneuvers with longer traveled distances. In the following, we quantify these differences by systematically comparing the statistics of the perbout kinematic parameters at different temperatures.
The bath temperature controls the statistical distributions of the kinematic parameters
For each batch and temperature, a probability density function (pdf) was computed for interbout intervals, displacements and turn angles by pooling all bout events. We then computed an average distribution across batches (Fig. 2B–D, respectively) for the 3 parameters, as well as the temperaturedependence of their mean values (Fig. 2F–H).
A decrease in the bath temperature from 26 to 18 ^{∘}C is associated with an increase of the mean IBI (〈δt〉) from 1 to 1.4 s, while the bout frequency remains essentially unchanged at higher temperatures (2B, F). This increase in the mean values is accompanied by a systematic broadening of the statistical distribution. The perbout displacement exhibits a similar trend (Fig. 2C). This quantity increases in the range 18−−26 ^{∘}C from 1 to 1.5 mm, and remains unchanged at higher temperatures (Fig. 2G).
The turn angle distributions shown in Fig. 2D reveal the existence of two main bout categories [13, 18, 24]. The central narrow peak corresponds to forward bouts while the wide tail is associated with turning events. We adjusted this distribution as a sum of two empirically chosen functional forms in order to extract the fraction of turning bouts p_{turn} (see “Methods” section). This quantity steadily increases with the temperature, from 0.3 to 0.8 (Fig. 2E). This increase in the fraction of turning bouts comes with an increase in their associated reorientation angles δθ_{turn} as shown in Fig. 2H.
The bath temperature controls the persistence time of the orientational state
In a recent study [18], we showed that the orientational dynamics of zebrafish larvae can be described by two independent Markov chains (Fig. 3A). The first one controls the bout type selection, between forward scoots or turn bouts. This process is essentially memoryless, such that the transition rates are simply set by the ratio between either categories, namely p_{turn} and 1−p_{turn}. The second Markov chain controls the orientations of the turning bouts. When a turn bout is executed and if this chain is in the left (right) state, then the animal turns left (right, respectively). This second selection process has been shown to display a persistence over a few bouts: the fish tends to chain turn bouts that are similarly orientated [24–26], a mechanism whose result is to enhance the angular exploration [18].
Here, we examined how this motorpersistence mechanism is impacted by the bath temperature. We estimated the flipping rate p_{flip} — the probability to switch orientation at each bout — by first binning the turning angles into three categories (denoted Δ) and assigning a discrete value to each of them: right turn (Δ=−1), forward bout (Δ=0) and left turn (Δ=+1). We then computed the mean discretized angle value 〈Δ_{n+1}〉 at bout n+1 for the three possible values of the previous bout Δ_{n}, as shown in Fig. 3B. The slope of the linear fit provides a measurement of p_{flip} (see “Methods” section and Eq. 1). This flipping probability increases with temperature from 0.22 at 18 ^{∘}C to 0.45 at 33 ^{∘}C (Fig. 3C), approaching 0.5. Hence, at high temperatures, the orientational persistence essentially vanishes, i.e., the probability to trigger a left vs a right turn becomes independent of the orientation of the previous bout. This increase of the flipping probability counteracts the concurrent increase in turning probability and turn bout amplitude at high temperature by limiting the tortuosity of the trajectories in these conditions.
This approach yields a typical number of bouts 1/p_{flip} over which the turning orientation is maintained. A complementary approach consists in characterizing the actual timepersistence (in seconds) of the orientational state [18]. To do so, we assume that the orientation selection is driven by a hidden twostate continuous signal, of which the turn bouts provide a stochastic sampling. We hypothesize that a forward bout is “transparent”, i.e., it does not interfere with the persistence process, and that the orientational state remains unchanged until a bout in the opposite direction is executed. The procedure for reconstructing the orientational signal is illustrated in Fig. 3D.
For all trajectories, we computed the autocorrelation function (ACF, R_{ΔΔ}) of the reconstructed orientational signals, and averaged them for each temperature (Fig. 3E). The ACF shows a faster decay for higher temperatures, i.e., the time period over which the animal can maintain its orientational state is larger in colder water. The ACFs could be correctly adjusted with an exponential decay, a functional form that is expected for a simple telegraph process [27]. This suggests that the left/right transition over a time interval dt is simply given by k_{flip}dt, where k_{flip} is the transition rate from one state to another. From the exponential fit of the ACFs, we extracted k_{flip}, which we found to increase quasilinearly with the temperature, as shown in Fig. 3F (purple line). The rate k_{flip} is the temporal counterpart of the perbout flipping rate p_{flip}, the two quantities being linked through the interbout interval. Consistently, we found that p_{flip}/〈δt〉 provides a good approximation of k_{flip} for all temperatures (Fig. 3F, red line).
Navigational kinematic parameters are statistically coupled
In the preceding sections, we showed that the bath temperature impacts in a systematic way the statistical distributions of the five kinematic parameters that control the fish spontaneous navigation, namely the interbout interval (IBI), turn amplitude, traveled distance, turn probability, and orientational flipping rate. When examining trajectories recorded at a given temperature, we noticed that they tend to fall in stereotypical categories reminiscent of those most often observed at various temperatures. Some trajectories are tortuous with short IBI, akin to typical hot water trajectories, while other appear to be straighter with less frequent bouts as generally observed in cold water (Figs. 4A and 1C). This is suggestive of the existence of a finite kinematic repertoire accessible to the animals whose relative occurrence may be controlled by the bath temperature.
To test this intuition, we first aimed at establishing the statistical constraints that could set this accessible repertoire. We thus examined the pairwise covariance of the aforementioned kinematic parameters. At short time scale (over one bout), we did not observe any significant correlation between the 3 parameters that can be evaluated on a perbout basis (IBI, reorientation angle, and traveled distance, see Additional file 1: Figure S1A). However, when performing the same analysis on pertrajectory averages, we observed a robust covariance of the parameters. This is illustrated in Fig. 4B which shows the covariance matrices computed for all data and for each temperature. The IBI appears to be strongly anticorrelated with the forward displacement and the flipping rate. In contrast, besides IBI, all pairs of parameter tend to exhibit positive correlations. Importantly, these statistical features are conserved across the entire temperature range.
Temperature controls the distribution probability within a welldefined locomotor repertoire
We sought to evaluate how this intratemperature covariance of the navigational parameters aligned with the intertemperature covariance. To do so, we used the temperatureaveraged parameters to build a 5 temperatures by 5 parameters matrix from which we computed an intertemperature Pearson correlation matrix (Additional file 1: Figure S1B). The latter displays a comparable structure as the mean intratemperature correlation matrix 4B: as we have shown in the previous sections, all parameters increase with temperature, and are thus positively correlated, except for the interbout interval which decreases with the temperature and is therefore anticorrelated with the 4 other parameters.
Hence, intrinsic variability and temperatureinduced behavioral changes both reflect a concerted modulation of the kinematic parameters along a similar axis. To confirm this claim, we performed a principal component analysis on both the intertemperature and intratemperature data. For all temperatures, the first principal component (PC) explains 28 to 45% of the intratemperature variance (Fig. 4C), i.e., significantly more than expected for independent parameters (20%). Due to the small size of the intertemperature matrix (5 samples), the first PC explains more than 90% of the intertemperature variance (Additional file 1: Figure S1C). The first PC is conserved across the temperature range (Fig. 4D, colored bars) and essentially aligned with the intertemperature PC (black squares). The second PC is similarly conserved across temperatures (Fig. 4E) yet less clearly aligned with its intertemperature counterpart.
In order to represent data from various temperatures within the same lowdimensional space, we performed a PCA analysis on the pooled covariance matrix, combining all intratemperature arrays after standardization (Fig. 4C, E, solid gray line). Based on the GuttmanKaiser criterion, we only retained the first two principal components [28] (Additional file 1: Figure S1D). Figure 4F shows the entire dataset projected in this unique 2D PCA space, where the temperature is colorcoded. As the temperature is increased, the accessible locomotor space is shifted towards higher values of both marginal projections, with a concurrent widening of the distribution for the first PC. These observations are thus in line with the view that the trajectories are confined to a manifold defined by the correlation between the various parameters. Each temperature delimits a specific accessible region of this subspace as defined by the PCs projection values.
Singlefish recordings reveal a slow diffusivelike modulation of the locomotor behavior
The experiments on which these analysis were performed are based on simultaneous recordings of 10 fish for each batch. As we cannot track individual fish over the entire session, we cannot evaluate to what extent individual animals’ navigational pattern may vary during the course of the assay. To address this specific question, we performed a second series of experiments in which single animals (N=18) were continuously monitored for 2h at an intermediate bath temperature of 26 ^{∘}C. The same analysis pipeline was implemented. In particular, the recordings were split into successive “trajectories” corresponding to walltowall sequences. We observed that over the course of the assay, the trajectories tended to exhibit strongly distinct features as illustrated in Fig. 5A, reflecting a significant intraindividual behavioral variability.
For each individual, we similarly computed a feature matrix containing, for all successive trajectories, the mean interbout interval, reorientation angle of turn events, displacement, turning probability, and flipping rate. We then performed a PCA on each array. Both the explained variance (Additional file 2: Figure S2A) and the PCA coefficients (Fig. 5B, C) were unchanged with respect to the multifish analysis (5B, C, gray line). This indicates that the covariance structure in the locomotion pattern is similar at the intra and interindividual level.
We thus used the multifish PC space defined in the preceding section to represent the singlefish data. The result for an example fish is shown in Fig. 5D where the successive trajectories are indicated as dots in this twodimensional PC space. This representation reveals a slow diffusivelike exploration of the locomotor space over the course of the experiment, with a progressive transition from one type of trajectory (e.g., long displacements, frequent bouts, frequent turns) to another (e.g., short displacements, longer interbout intervals, and fewer turns).
To quantify the timescale of this itinerant exploration within the locomotor space, we computed the autocorrelation function (ACF) of the projections on the two first PCA components (5E, F, black line). These curves could be captured by an Ornstein–Uhlenbeck (O.U.) process, which describes the dynamics of a random walker within a quadratic energy basin [29, 30], see “Methods” section). The latter allows one to bound the stochastic exploration within a finite region of the locomotor space. From the fit, we extracted the times needed for the dynamical system to reach its stationary regime: τ=2585±58 s for PC1, τ=1980±14 s for PC2 (mean ± s.e.m.). These values clearly demonstrate that the modulation of the exploratory behavior in individual animals takes place over time scales that are orders of magnitude longer than the interbout interval.
This series of experiments allowed us to further assess the relative contribution of the intra and interindividual components in the observed behavioral variability. As the assay is longer (2 h) than the time needed to reach the stationary regime (∼ 2000 s), each recording provides an estimate of the intraindividual variability. The latter was quantified in the PC space as the variance of the PC projections across the entire duration of the assay, averaged over the various individuals. We then separately computed the variance of the PC projections, pooling the data of all animals (Additional file 2: Figure S2D, green). The latter quantity thus encompasses both inter and intraindividual variability. This analysis led to the conclusion that a dominant fraction of the variance (68% on PC1, 53% on PC2) can be explained by the intraindividual variability.
Simulations of spontaneous navigation at various temperatures reveal basic thermophobic behavior without direct gradientsensing mechanism
Having thoroughly characterized the statistical structure of the kinematic parameters and their thermal modulation, we sought to build a minimal stochastic model of the fish navigation in order to generate synthetic trajectories at different temperatures. Each kinematic parameter defines a random variable whose mean is set by the temperature and whose statistical distribution accounts for both the intertrajectory variability and the perbout stochasticity. The dual nature of the variability was mathematically recapitulated by expressing each of the 5 kinematic variables as a product of two stochastic, temperatureindependent variables: one accounting for the trajectorytotrajectory modulation (within a range controlled by the bath temperature, Additional file 3: Figure S3BE, Y column), and the other reflecting the remaining shortterm variability (bouttobout, Additional file 3: Figure S3BE, ε column, see “Methods” section). For the former, we used the copula method to reproduce the observed covariance of the pertrajectory means of the various parameters.
This approach allowed us to generate various trajectories at different temperatures, as illustrated in Fig. 6A. These trajectories are qualitatively similar to those typically observed at the corresponding temperatures (see Fig. 1C for a comparison). To quantify how this stochastic model captures the exploratory behavior, we computed the mean square displacement (MSD, Fig. 6B) and the mean square reorientation (MSR, Fig. 6C) on both the real (dots) and numerical data (solid lines). Overall, the exploratory dynamics appear to be correctly reproduced by the numerical model. Importantly, the intertrajectory variability is also, by construction, correctly reproduced by this minimal model.
This model was used to probe how the temperature dependence of the navigational kinematics may participate in driving the animal along thermal gradients. We first experimentally quantified how zebrafish larvae responded to a linear thermal gradient spanning our temperature range (18−−33 ^{∘}C), by focusing on the steadystate occupation distribution. We found that the larvae favor regions where the temperature is comprised between 23 and 29 ^{∘}C (Additional file 4: Figure S4), i.e., they tend to avoid both extreme (hot and cold) regions. The underlying sensorymotor mechanism is bound to involve both the effect of the temperature on the fish navigation pattern (thermokinesis) and a direct (immediate) response to perceived temperature changes (thermotaxis) [14, 16]. Our model allows us to assess the relative contribution of the kinesis process. In order to do so, we implemented a simulation in which a virtual fish navigates in a rectangular pool (L×45 mm) in which we imposed a linear thermal gradient along the horizontal xaxis spanning the 18−−33 ^{∘}C range. We simulated trajectories of numerical swimmers by continuously updating their exploratory statistics according to the local bath temperature. These changes are entirely controlled by the temperaturedependence of the 5 kinematic parameters, which we linearly interpolated across the thermal gradient. Four gradient strengths were emulated by changing the length L of the pool (L=0.1,0.3,0.5,1 m).
The time evolution of the position distribution along the gradient are shown as heatmaps in Fig. 6D. They reveal a global drift of the population towards the low temperature region for all values of the thermal gradient (Fig. 6E). In all conditions, the distributions were found to converge towards a unique steadystate profile after a finite time. The probability of presence in the steadystate regime displays a quasilinear decay from 18 to 26 ^{∘}C, and remains uniform at higher temperature. The thermokinesis process thus endows the animal with a basic thermophobic behavior, even for minute gradients — orders of magnitude smaller than those imposed in thermotactic assays. In contrast, the avoidance of cold regions seen in experiments (Additional file 4: Figure S4, see “Methods” section) is absent in our simulations and must therefore reflect a direct gradientsensing mechanism.
The dynamic of this thermophobic behavior in the simulations appears to depend on the imposed gradient, as illustrated in Fig. 6F, which shows the mean experienced temperature across the population as a function of time for the three gradients. All the curves display a similar decay associated with a global drift towards the cold region, until a similar plateau value is reached, albeit with different timescales. Due to the diffusive nature of the fish spatial exploration, the settling time is expected to scale with the square of the pool length. Consistently, the four dynamic evolution are found to fall on a unique curve when plotted as a function of t/L^{2} (Fig. 6G). The associated settling times range from 10 min for the largest gradient up to ∼ 14 h for the smallest one.
Discussion
Animal behaviors unfold as trajectories in a high dimensional space of motor actions. To make behavior mathematically tractable, one needs to unveil statistical rules that couple the different components of the behavior and organize them across timescales. This dimensionality reduction approach is a prerequisite to further distinguish between deterministic and stochastic components of the behavior and concurrently discover the underlying neural mechanisms [31, 32]. Leveraging novel techniques for highthroughput behavioral monitoring and automatic classifications has allowed to elucidate the statistical structure organizing selfgenerated behaviors in numerous species, such as C. elegans [33], Drosophila [34, 35], zebrafish [21, 36], or mice [37].
With its boutbased navigation, zebrafish larva offers a relatively simple model for such an endeavor. It has been shown that as few as 13 different swim bout types are sufficient to capture the entirety of its behavioral repertoire [21]. Here we focus on spontaneous exploration in the absence of timevarying sensory cues. For this particular behavior, Marques et al. showed that only 4 different bout types effectively contribute to the navigation. In the present work, we further restricted the discretization to only two categories (forward and turn bouts), and extracted 5 kinematic variables, ignoring fine differences in tail bout execution. We showed that the knowledge of these 5 variables statistics nevertheless suffices to fully characterize the longterm exploratory process. Indeed, synthetic trajectories generated by stochastic sampling from the statistical distributions extracted from the data accurately reproduce the experimentally observed angular and translational dynamics.
Using this reductionist approach, we were able to demonstrate that the variability in the fish exploratory dynamics originates from two separate mechanisms, acting on distinct timescales. Over a few bouts, the fish displacement is akin to a random walk in which multiple stochastic processes set the successive values of two discrete (bout type and turn bout orientation) and three continuous (interboutinterval, linear, and angular displacements) variables that together define its instantaneous inplane velocity. These processes are statistically constrained by mean transition rates and amplitude probability distributions that can be considered invariant at the scale of individual trajectories (i.e., over tens of bouts). These parameters however vary significantly over long time scales: their time modulation takes place over hundreds to thousands of bouts, indicative of a clear timeseparation between the two different processes. Importantly, although we did not observe any significant correlation in the instantaneous locomotor variables, the slow modulation of the kinematic parameters exhibits robust covariance and is thus constrained within a welldefined kinematic manifold.
The present study allowed us to quantify how the water temperature modulates the locomotor statistics of zebrafish larvae. Rather than evoking distinct locomotor patterns, temperature controls the relative occupancy within this subspace: changing the temperature consistently impacts the mean value of the kinematic parameters but leaves their covariance structure unchanged. Temperature thus essentially sets the accessible range of exploratory trajectories within a welldefined continuum of possible locomotor behaviors.
At the circuit level, it is tempting to interpret these observations by considering the brain as a dynamical system exhibiting multiple metastable patterns of activity (brain states) whose relative stability and transition rates define a particular energy landscape [38]. In this view, the shorttime dynamics that select the successive bout properties correspond to a stochastic itinerant exploration of this neuronal landscape. The latter is essentially invariant over minutes but is slowly reshaped via endogenous processes or through temperature changes, leading to a gradual modification of the shortterm statistics.
Slow modulation of locomotor characteristics in zebrafish larvae have been reported in two recent studies [2, 36]. In [2], the authors identified two discrete states, associated with exploration and exploitation during foraging, with typical persistent times of order of minutes. In [36], progressive changes in locomotor statistics were associated with decaying hunger state, as the initially starved animal progressively reached satiety. In contrast with these two studies, the modulation in locomotor kinematics that we observed is continuous and does not reflect spatial heterogeneities in the environment (e.g., local presence of preys) or explicit changes in internal states such as satiety. With respect to hunger state, the use of temperature may offer a practical way to externally drive the internal state to a stationary point in an ethologically relevant way.
The neuronal basis of this internal state modulation process remains to be elucidated. The circuits regulating specific locomotor features, such as the bout frequency [39] or orientational persistence [25, 26] have been identified. However, the fact that the various kinematic parameters display concerted endogenous modulations points towards a global drive. Temperature is known to impact cellular and synaptic mechanisms [40] in such a way that an increase in temperature tends to speed up neuronal oscillatory processes [41, 42]. This may explain the concurrent decrease in the persistent times associated with the orientational persistence and interbout intervals. The thermal modulation of the angular and linear amplitude of the bouts may in turn reflect a temperature dependence of the muscular efficiency rather than neuronal processes [43]. Another possibility is that the temperature drives the activity level of neuromodulatory centers which may also exhibit slow endogeneous modulations. This neuromodulation release would then globally impact the spontaneous dynamics of various premotor centers yielding the observed change in locomotor patterns. The serotonergic neurons of the dorsal raphe constitute an attractive candidate for such a mechanism as their activation has been shown in numerous instances to drive a persistent change in behavior in zebrafish [2, 5, 44], as well as in mice [45].
Our study yields a minimal numerical model of zebrafish locomotion at different temperatures. This model allowed us to probe in silico how the thermal modulation of the exploratory dynamics may contribute to the thermotaxis behavior, thus complementing direct gradientsensing mechanisms [17]. Our simulations indicate that this thermokinesis process endows the animal with the capacity to efficiently avoid hotter regions, but cannot explain the observed avoidance of cold water. As thermal gradient sensing operates within a time window of 400 ms [16], it may be ineffective in conditions where the lengthscale of thermal gradients is much larger than the typical distance traveled per bout. In such conditions, this complementary mechanism may be strategically relevant as it allows the animal to navigate away from potentially noxious regions.
Conclusions
This study establishes the temperature as an effective and practical external parameter to explore behavior variability in vertebrates. Our analysis provides simple latent variables, namely the two first PCA projections, that can be used to efficiently track the animal’s behavioral state. Changes in behavioral states are generally induced through complex protocols, involving a perturbation of a sensorimotor loop, or through abrupt changes in sensory conditions [46]. In such approaches, the change is discrete and generally transient as the animal eventually adapts to the new conditions. In contrast, temperature offers a way to drive a robust, continuous, and chronic shift in behavior that can be easily implemented while performing largescale brain monitoring. Various behavioral states are thought to reflect different levels of attention or arousal, which in turn impact the responses to sensory stimulation. Beyond its utility for studying how a given neuronal circuit may give rise to distinct dynamics, thermal perturbation could also be leveraged to investigate how internal states may enhance or inhibit sensory responses.
Methods
Animals maintenance
Experiments were performed with wild type Danio rerio AB, aged 5 to 7 days postfertilization (dpf). Larvae were reared in Petri dishes containing embryo medium (E3), at 28 ^{∘}C, with a 14/10h cycle of light/dark and were fed with nursery powder GM75 everyday from 6dpf. Experiments were done during daytime, in E3. They were approved by Le Comité d’Éthique pour l’Expérimentation Animale Charles Darwin C2EA05 (02601.01).
Experimental setup
The container consists in a rectangular pool (100 × 45 × 2.5 mm^{3}) made of copper whose surface was protected by a biocompatible heatresistant black layer (RustOleum). It is stuck on two 78W Peltier modules (Adaptive) with thermal tape (3M). A transparent, 2mmthick PMMA cover is placed over the pool with 2 mm spacers to minimize water evaporation, leaving a water thickness of 4.5 mm. The temperature is measured at both ends of the pool with thermocouples type T (Omega). The two left/right error signals (T_{target}−T_{measured}) are used within two independent PID loops implemented on an Arduino Uno board (Arduino) whose coefficients have been optimized manually. Each PID regulates the PWM frequency sent to a Hbridge driving the power sent to the two Peltier modules. A graphical user interface (GUI) written in C++ using the Qt framework is used to monitor the measured temperatures and to impose the target temperatures on both ends. Due to its high thermal diffusivity, the copper piece quickly reaches a uniform temperature and acts as a thermostat for the water. After about 4 min, the temperature of the water in the center of the pool has reached the set temperature (± 0.2 ^{∘}C), which then remains constant over time. The GUI monitors the bath temperatures while grabbing frames from a CMOS camera (FLIR Chameleon3 CM3U313Y3MCS) coupled with a macrolens (Navitar) at 25 frames per second. The whole apparatus is placed in a lighttight box, illuminated with a homogeneous white light emitted by a LED panel (Moritex) placed on the side; a mirror placed at the other side limits significant phototactic bias in the small direction of the pool. All codes mentioned above are available on Github [47] under a GNU GPLv3 license. Blueprints of the box and pool as well as electronic designs are available upon request.
Innocuity of the blackpainted copper pool
Zebrafish larvae are sensitive to minute concentration of chemicals. To check the innocuity of the container, ten zebrafish larvae were left overnight inside the setup. All survived and were swimming actively. We further checked whether the black layer may releases chemical compounds that could impact the animals navigational dynamics. We prepared a stock of E3 heated for 1 h at 45 ^{∘}C in the experimental setup. Three batches of 10 larvae were then sequentially recorded for 30 min in E3 (control), in a Petri dish, then placed for 40 min in the heated E3 and recorded again for 30 min (treated). Both experiments were performed at 28 ^{∘}C. We then computed the mean interbout interval, mean displacement, mean reorientation amplitude for each batch, and the cumulative mean square displacement. Additional file 5: Figure S5A shows that none of these quantities display any significant change for larvae recorded in the heated water with respect to the control. Copper being known to alter lateral line neuromasts [48], we further tested the absence of copper ions by incubating 10 larvae in the heated E3 for 2 h before marking their lateral line neuromasts with DiAsp according to the protocol detailed in [48]. All larvae displayed intact lateral line neuromasts.
Experimental protocols
The pool is filled with E3. A temperature is randomly drawn from 18,22,26,30,33 ^{∘}C and set with the GUI. After 10 min, a batch of 10 zebrafish larvae is introduced in the pool. After 10 min of habituation, the fish kinematics are monitored for 1800 s (half an hour). In order to confirm the full habituation of the animal to the new conditions, we checked that the fish navigation is indeed timeinvariant during the recording period. We evenly split each recording into 3 time windows and computed the statistics of the various kinematic parameters (mean bout frequency, displacement, and reorientation amplitude) for these 3 periods, The distributions within each timewindow are not significantly different (p>0.1, twosample KolmogorovSmirnov test) as shown in Additional file 5: Figure S5B. Fish remain in the pool while we randomly draw a new nontested temperature. After 20 min (temperature regulation and habituation), a new recording of 1800 s is performed. The five temperatures are not systematically tested on all batches, but for each temperature, 10 different batches of 10 fish are used. To check whether the testing order may impact the kinematic of swimming for a given temperature, we separated the assays in two groups depending on whether the previous recording was performed at a lower (ΔT>0) versus higher (ΔT<0) temperature. We then computed the mean interbout intervals, mean displacement and mean reorientation amplitude normalized by their temperature mean. The comparison of the two distributions is shown in Additional file 5: Figure S5C. For all parameters, the statistical difference is nonsignificant, indicating that the testing order has no major impact on the navigational statistics. In total, the experiments involved 17 different batches. The sample size was not statistically determined beforehand.
For singlefish experiments, the same protocol is used except that a single fish was placed in the pool. The recordings last for 2 h and only T=26 ^{∘}C is tested.
For thermal gradient experiments (Additional file 4: Figure S4), 10 larvae are used during 45 min. The first 5 min are recorded with a uniform temperature of 22 ^{∘}C, then a linear gradient is imposed during 40 min, from 18 to 33 ^{∘}C. The gradient direction (i.e., which side is set to either 18 or 33 ^{∘}C) is chosen randomly. 10 different batches are tested. The distribution of presence along the gradient is computed over the last 2 min (5% of the gradient duration) such as to allow enough time for the animals to reach a steadystate.
Tracking and basic analysis
Larvae were tracked offline using the opensource FastTrack software [22]. It generates a text file containing the position of each fish’s center of mass and body angle across frames until they leave the defined ROI. Kinematic analyses were performed using MATLAB (R2020a, Mathworks). Bouts are detected when the instantaneous speed is greater than two times the overall variance of the speed. Putative bouts are then filtered on a distance criterion (bouts with a linear displacement — measured in a time window of ± 0.5 s centered on the bout onset — less than 0.3 mm or greater than 18 mm are rejected) and on a temporal criterion (bouts occurring within 0.4 s after a bout are rejected). Bout timing is defined as 80 ms before the velocity peak. Detection performance was checked manually on randomly selected sequences. From positions, time, and body angles before and after a bout event, we computed displacements, interbout intervals, and turn angles associated with each bout. Data are split into trajectories, from one edge of the ROI (set at 5 mm from the walls) to another. Only trajectories that last at least 25 s, with at least 10 bouts, with 3 bout types (left turn, right turn, and forward scoot) are kept for further analysis. Trajectories last on average 67 s (median 47 s, 95th percentile 178 s) and contain on average 60 bouts (median 44 bouts, 95th percentile 154 bouts). The number of trajectories and the number of bouts retained for further analysis are given for each temperature in Table 1. The effect of the chosen cutoffs (minimum trajectory duration and minimum interbout interval) used for trajectory and bout selection is tested in Additional file 6: Figure S6. This control demonstrates that changing these cutoffs has no significant impact on our results. All MATLAB routines are available on Gitlab [49] under the GNU GPLv3 license.
Bout classification
To discriminate whether a bout falls in the forward or the turning categories, we fitted the onesided (absolute value) reorientation angles distributions with the sum of a zeromean Gaussian distribution and a gamma distribution. The Gaussian corresponds to the part of the distribution close to zero, while the gamma function aims at describing the distribution of high angle reorientations. We manually set the Gaussian width and the scale parameter of the gamma function based on the observed distributions. We fitted the shape parameter for each temperature, ensuring that the slope at high angles in logarithmic scale is well reproduced. Then, we defined a fixed threshold for the angles to be considered as a turn or a forward bout. This threshold is the angle at which the two distributions cross, invariably found around 10^{∘} (10.25±0.23^{∘},mean±sd). This value of 10^{∘} (0.17 rad) was used to classify bouts throughout this work.
Displacement correction
We noticed that the displacement corresponding to a turn event was systematically larger than the displacement associated to a forward event. This is due to the fine structure of a turning bout: first, the fish performs a small reorienting bout, then it scoots forward [20]. Since we do not look at this fine structure, the overall displacement during a turn bout is geometrically overestimated and would bias temperaturetotemperature comparison. We computed the ratio between displacements during turns and the ones during forward swims, and found a factor of 1.6±0.1, regardless of the temperature. Therefore, in all analyses presented in this work, all displacements corresponding to a turn event were corrected by a factor 1/1.6=0.625.
Statistical methods
Probability density functions (pdf) were computed with a kernel density estimation through the builtin Matlab function ksdensity, with a bandwidth of 0.1 for interbout intervals and displacements and 0.5 for turn angles. For the distributions of Fig. 2, a pdf was computed for each batch and the mean and standard error of the mean are computed. For rescaled curves (Additional file 3: Figure S3), data from all experiments were pooled to compute the temperatureaverage quantity \(\overline {X}_{T}\) and rescaled values. Boxplots were made with the builtin Matlab function boxchart, using as input data the means of the respective quantities for trial (one dot corresponds to a batch of 10 fish). For simulations of navigation, averages over temperature were computed by pooling all bout events from all experiments corresponding to this particular temperature. p_{turn} and p_{flip} values were estimated for each trajectory and then averaged. Error bars for those temperature averages and for the pdf shown in Additional file 3: Figure S3 were all computed using bootstrapping with 1000 boots to get the 95% confidence interval through the builtin bootci function. Errors were propagated for the ratio of p_{flip} and 〈δt〉_{T} in Fig. 3F.
Reorientation dynamics
The two Markov chains model has been described in details in a previous study [18]. We first binned the reorientation angles δθ into a ternarized reorientation Δ, with values −1 (right turn R), 0 (forward bout F) and +1 (left turn L). To extract p_{flip}, we analytically derived the mean reorientation Δ_{n+1} given the previous reorientation Δ_{n}. There are 9 combinations of bouts {n;n+1}: {L;L},{L;R},{L;F},{F;L},{F;R},{F;F},{R;L},{R;R},{R;F}. All combinations involving a forward bout yield 0. Remain combinations with two turns in the same direction and two turns in the opposite direction. For a turn in direction L (resp. R), the associated probability corresponds to the case where a flip occurred (i.e., the previous bout was in direction R, resp. L) and the case where no flip occurred (i.e., the previous bout was in direction L, resp. R). Noting \(\Delta _{n}^{R}\) and \(\Delta _{n}^{L}\) the turns in the right and left direction at bout n, the mean reorientation given the direction of the previous bout reads:
These equations can be summed up as:
This is the fit used in Fig. 3B.
A random telegraph signal is a binary stochastic process with constant transition probability per unit of time. In the case where both states are equiprobable, the two transition rates (here noted k_{flip}) are equal. For such processes, the time spent in one or the other state (left or right) is exponentially distributed [27] and the autocorrelation function for a zeromean signal reads :
This is the fit used in Fig. 3E.
Mean square displacement (MSD) <d^{2}> and mean square reorientation (MSR) 〈δθ^{2}〉 were computed using the MATLAB package msdanalyzer [50]. All (x,y) and δθ sequences are pooled by temperature for both data and simulations, the MSD and MSR were computed for each sequence and we show in Fig. 6B and C the ensemble average for each temperature with the standard error of the mean.
Principal components analysis
The “features matrices” were built for each temperature. They include, for each trajectory, mean interbout intervals, turn probability, flip rate (estimated as p_{flip}/〈δt〉,p_{flip} being extracted as explained above, for each trajectory), mean reorientation angle of turning events and mean displacements. Each set was standardized (centered and normalized by its standard deviation) before being processed by the single value decomposition (SVD) algorithm through the builtin pca function. Those 5 intratemperature standardized arrays are then concatenated to form the socalled pooled matrix, that is in turn used to find a common space through PCA. For projection, each set was normalized by the standard deviation of all the pooled data (regardless of temperature) and not centered for comparison purposes. The aforementioned common space was also used to project data from singlefish experiments.
Numerical OrnsteinUhlenbeck process
The singlefish experiments contains an average of 48±16 trajectories (mean ± sd) with a median duration of 42 s and a value of the 90th percentile of 122 s. Each trajectory yields one point in the PC space. We then linearly interpolated the PCA projection values in order to produce a discrete signal defined over the same time vector across the total duration of the assay (7200 s), sampled every second. For each fish, on both PC, we computed the autocorrelation function (Additional file 2: Figure S2BC) before averaging them (Fig. 5E, F, black line is the mean, shade is the s.e.m.). The autocorrelation function of a OU process reads [30]:
where D is the diffusion coefficient and k is the bias term. However, this expression is only valid in the limit where the recording time is much larger than the relaxation time 1/k of the process. When the recording duration is of order of the autocorrelation decay time, the computed autocorrelation function exhibits a negative overshoot, which reflects an incorrect estimate of the longtime mean of the signal. This issue is common to stochastic signals whose mean is unknown. In order to fit the experimental autocorrelation and extract the relaxation time τ=1/k, we used simulated dynamics over similar timewindows. Numerical simulations of the Ornstein–Uhlenbeck (O.U.) process were sequentially implemented using the following equation [51] :
with δt the time interval chosen for the simulation (units s) and \(\mathcal {N}\) is a random number drawn from a normal distribution.
To determine τ, we generated 500 realizations of the O.U. process with D set to 1 and τ set to values in a given range. For each realization, we computed the autocorrelation function (ACF) and averaged them across realizations. We then computed the residual sum of square (RSS) and chose the minimum one to select the best parameter τ. After manually narrowing down the best range for τ (PC1, 2000 to 3000 s, 1000 values; PC2, 1900 to 2100 s, 1000 values), we repeated the previous process 20 times to get 20 “best τ” and we report the mean ± s.e.m. in the text and figure.
It should be noticed that the autocorrelation decay times extracted from this analysis are 2585 and 1980 s, for PC1 and PC2 respectively. These timescales are thus one order of magnitude larger than the typical trajectory duration. This clear timescale separation a posteriori validates the pertrajectory discretization used in our analysis.
Numerical simulations of trajectories
Trajectories were simulated using the framework described in Additional file 3: Figure S3, based on the hypothesis that (1) spatiotemporal dynamics can be reproduced solely from five parameters, (2) perbout values of interbout intervals (δt), displacements (d) and turn angles (δθ) are drawn from a distribution that can be decomposed as \(X=\overline {X}_{T}Y\epsilon \), (3) the pertrajectory values of turning probability (p_{turn}) and flipping probability (p_{flip}) are drawn from a distribution that can be decomposed as \(X=\overline {X}_{T}Y\), and (4) the trajectoryaveraged parameters are correlated. Note that for the simulations we use p_{flip} rather than flipping rate for simplicity in the code implementation.
\(\overline {X}_{T}\), the temperature average
All perbout values of δt, d, reorientation angle of turn events (δθ_{t}) and reorientation angles of forward events (δθ_{f}) are pooled by temperature and the mean is computed. A p_{turn} and a p_{flip} is estimated for each trajectory, pooled by temperature and averaged (Additional file 3: Figure S3BE, left column).
Y, the trajectory means variability
For each trajectory, a mean value is computed for δt, d and δθ_{t/f} while p_{turn} and p_{flip} are extracted. They are then rescaled by the corresponding temperature average value computed above. For each temperature, a cumulative density function (cdf) is computed. They are then averaged across temperatures to get a single Y cdf for each parameters (pdf shown in Additional file 3: Figure S3BE, middle column).
ε, the perbout variability
Similarly, for each trajectory we rescale values of δt, d and δθ_{t/f} by their corresponding trajectory mean. Then, all events are pooled by temperature and a cdf is computed. Finally, we will use the mean cdf, resulting in a single ε cdf for perbout parameters. p_{turn} and p_{flip} are defined for a trajectory, hence they do not have bouttobout variability (pdf shown in Additional file 3: Figure S3BD, right column).
Correlations of means
We compute the Pearson’s correlation matrix of the trajectories’ parameters (trajectory means and probabilities), for each temperature. The coefficients are then averaged to get a single correlations matrix 〈R_{traj}〉_{T}.
Algorithm
After choosing a number n of fish (trajectories), we generate multivariate distributions (copulas) with the MATLAB builtin mvnrnd function, with the mean 〈R_{traj}〉_{T} correlations matrix as input. It produces 5 marginal sets of n Gaussian random numbers, correlated with one another. We then get the corresponding normal cdf, which is in turn used to sample the corresponding Y cdfs, inversing the latter. Finally, those samples are multiplied by the corresponding temperature average \(\overline {X}_{T}\). A bout is generated by sampling a displacement and a turning angle, along with a interbout interval during which the virtual fish stands still, from the generic cdf of ε. Those values are multiplied by the trajectory means drawn earlier, and the new position (x,y) is computed. The next bout is generated, and so on. For the gradient simulations, the same strategy is used, at the notable difference that the temperature averages are determined dynamically given the position of the agent along the temperature gradient. We used reflective boundary conditions. We checked the consistency between parameters distributions from the data and from the simulations, as well as correlations between trajectory means.
Availability of data and materials
All data generated or analyzed during this study are included in the Dryad repository https://doi.org/10.5061/dryad.3r2280ggw.
Abbreviations
 IBI:

Interbout interval
References
 1
Shaw AK. Causes and consequences of individual variation in animal movement. Mov Ecol. 2020; 8(1):12.
 2
Marques JC, Li M, Schaak D, Robson DN, Li JM. Internal state dynamics shape brainwide activity and foraging behaviour. Nature. 2020; 577(7789):239–43.
 3
Gillespie JH. Natural selection for whithingeneration variance in offspring number. Genetics. 1974; 76(3):601–6.
 4
Philippi T, Seger J. Hedging one’s evolutionary bets, revisited. Trends Ecol Evol. 1989; 4(2):41–44.
 5
LovettBarron M, Andalman AS, Allen WE, Vesuna S, Kauvar I, Burns VM, Deisseroth K. Ancestral circuits for the coordinated modulation of brain state. Cell. 2017; 171(6):1411–142317.
 6
McCormick DA, Nestvogel DB, He BJ. Neuromodulation of brain state and behavior. Ann Rev Neurosci. 2020; 43(1):391–415.
 7
MacDonald SWS, Nyberg L, Bäckman L. Intraindividual variability in behavior: links to brain structure, neurotransmission and neuronal activity. Trends Neurosci. 2006; 29(8):474–80.
 8
Hiesinger PR, Hassan BA. The evolution of variability and robustness in neural development. Trends Neurosci. 2018; 41(9):577–86.
 9
Panier T, Romano SA, Olive R, Pietri T, Sumbre G, Candelier R, Debrégeas G. Fast functional imaging of multiple brain regions in intact zebrafish larvae using selective plane illumination microscopy. Front Neural Circ. 2013; 7:65.
 10
Ahrens MB, Orger MB, Robson DN, Li JM, Keller PJ. Wholebrain functional imaging at cellular resolution using lightsheet microscopy. Nat Methods. 2013; 10(5):413–20.
 11
Portugues R, Feierstein CE, Engert F, Orger MB. Wholebrain activity maps reveal stereotyped, distributed networks for visuomotor behavior. Neuron. 2014; 81(6):1328–43.
 12
Burgess HA, Granato M. Sensorimotor gating in larval zebrafish. J Neurosci. 2007; 27(18):4984–94.
 13
Burgess HA, Granato M. Modulation of locomotor activity in larval zebrafish during light adaptation. J Exp Biol. 2007; 210(14):2526–39. Chap. Research Article.
 14
Haesemeyer M. Thermoregulation in fish. Mol Cell Endocrinol. 2020; 518:110986.
 15
Engeszer RE, Patterson LB, Rao AA, Parichy DM. Zebrafish in the wild: a review of natural history and new notes from The Field. Zebrafish. 2007; 4(1):21–40.
 16
Haesemeyer M, Robson DN, Li JM, Schier AF, Engert F. The structure and timescales of heat perception in larval zebrafish. Cell Systems. 2015; 1(5):338–48.
 17
Haesemeyer M, Robson DN, Li JM, Schier AF, Engert F. A brainwide circuit model of heatevoked swimming behavior in larval zebrafish. Neuron. 2018; 98(4):817–8316.
 18
Karpenko S, Wolf S, Lafaye J, Le Goc G, Panier T, Bormuth V, Candelier R, Debrégeas G. From behavior to circuit modeling of lightseeking navigation in zebrafish larvae. eLife. 2020; 9:52882.
 19
Gau P, Poon J, UfretVincenty C, Snelson CD, Gordon SE, Raible DW, Dhaka A. The zebrafish ortholog of TRPV1 is required for heatinduced locomotion. J Neurosci. 2013; 33(12):5249–60.
 20
Budick SA, O’Malley DM. Locomotor repertoire of the larval zebrafish: swimming, turning and prey capture. J Exp Biol. 2000; 203(17):2565–79.
 21
Marques JC, Lackner S, Félix R, Orger MB. Structure of the zebrafish locomotor repertoire revealed with unsupervised behavioral clustering. Curr Biol. 2018; 28(2):181–1955.
 22
Gallois B, Candelier R. FastTrack: an opensource software for tracking varying numbers of deformable objects. PLOS Comput Biol. 2021; 17(2):1008697.
 23
Schnörr SJ, Steenbergen PJ, Richardson MK, Champagne DL. Measuring thigmotaxis in larval zebrafish. Behav Brain Res. 2012; 228(2):367–74.
 24
Chen X, Engert F. Navigational strategies underlying phototaxis in larval zebrafish. Front Syst Neurosci. 2014; 8:39.
 25
Dunn TW, Mu Y, Narayan S, Randlett O, Naumann EA, Yang CT, Schier AF, Freeman J, Engert F, Ahrens MB. Brainwide mapping of neural activity controlling zebrafish exploratory locomotion. eLife. 2016; 5:12741.
 26
Wolf S, Dubreuil AM, Bertoni T, Böhm UL, Bormuth V, Candelier R, Karpenko S, Hildebrand DGC, Bianco IH, Monasson R, Debrégeas G. Sensorimotor computation underlying phototaxis in zebrafish. Nat Commun. 2017; 8(1):651.
 27
Balakrishnan V. Mathematical physics: applications and problems. Cham: Springer International Publishing; 2020.
 28
Kaiser HF. The application of electronic computers to factor analysis. Educ Psychol Meas. 1960; 20(1):141–51.
 29
Uhlenbeck GE, Ornstein LS. On the theory of the Brownian motion. Physical Review. 1930; 36(5):823–41.
 30
Gardiner CW. Handbook of stochastic methods: for physics, chemistry and the natural sciences, Study ed., 2. ed., 6. print edn. Berlin: Springer; 2002.
 31
Krakauer JW, Ghazanfar AA, GomezMarin A, MacIver MA, Poeppel D. Neuroscience needs behavior: correcting a reductionist bias. Neuron. 2017; 93(3):480–90.
 32
Datta SR, Anderson DJ, Branson K, Perona P, Leifer A. Computational neuroethology: a call to action. Neuron. 2019; 104(1):11–24.
 33
Stephens GJ, JohnsonKerner B, Bialek W, Ryu WS. Dimensionality and dynamics in the behavior of C. elegans. PLoS Comput Biol. 2008; 4(4):1000028.
 34
Branson K, Robie AA, Bender J, Perona P, Dickinson MH. Highthroughput ethomics in large groups of drosophila. Nat Methods. 2009; 6(6):451–7.
 35
Mueller JM, Ravbar P, Simpson JH, Carlson JM. Drosophila melanogaster grooming possesses syntax with distinct rules at different temporal scales. PLOS Comput Biol. 2019; 15(6):1007105.
 36
Johnson RE, Linderman S, Panier T, Wee CL, Song E, Herrera KJ, Miller A, Engert F. Probabilistic models of larval zebrafish behavior reveal structure on many scales. Curr Biol. 2020; 30(1):70–824.
 37
Wiltschko AB, Johnson MJ, Iurilli G, Peterson RE, Katon JM, Pashkovski SL, Abraira VE, Adams RP, Datta SR. Mapping subsecond structure in mouse behavior. Neuron. 2015; 88(6):1121–35.
 38
Mazzucato L, Fontanini A, La Camera G. Dynamics of multistable states during ongoing and evoked cortical activity. J Neurosci. 2015; 35(21):8214–31.
 39
Severi KE, Portugues R, Marques JC, O’Malley DM, Orger MB, Engert F. Neural control and modulation of swimming speed in the larval zebrafish. Neuron. 2014; 83(3):692–707.
 40
Hille B. Ion channels of excitable membranes, 3rd edn: Sinauer Associates; 2001.
 41
Tang LS, Taylor AL, Rinberg A, Marder E. Robustness of a rhythmic circuit to short and longterm temperature changes. J Neurosci. 2012; 32(29):10075–85. Chap. Articles.
 42
Rinberg A, Taylor AL, Marder E. The effects of temperature on the stability of a neuronal oscillator. PLOS Comput Biol. 2013; 9(1):1002857.
 43
Batty RS, Blaxter JHS. The effect of temperature on the burst swimming performance of fish larvae. J Exp Biol. 1992; 170(1):187–201. Chap. Journal Articles.
 44
Yokogawa T, Hannan MC, Burgess HA. The dorsal raphe modulates sensory responsiveness during arousal in zebrafish. J Neurosci. 2012; 32(43):15205–15.
 45
Lottem E, Banerjee D, Vertechi P, Sarra D, oude Lohuis M, Mainen ZF. Activation of serotonin neurons promotes active persistence in a probabilistic foraging task. Nat Commun. 2018; 9(1):1000.
 46
Andalman AS, Burns VM, LovettBarron M, Broxton M, Poole B, Yang SJ, Grosenick L, Lerner TN, Chen R, Benster T, Mourrain P, Levoy M, Rajan K, Deisseroth K. Neuronal dynamics regulating brain and behavioral state transitions. Cell. 2019; 177(4):970–98520.
 47
Candelier R, Le Goc G. ThermoMaster. Github. bfb329d. https://github.com/LJPZebra/ThermoMaster.
 48
Hernández PP, Moreno V, Olivari FA, Allende ML. Sublethal concentrations of waterborne copper are toxic to lateral line neuromasts in zebrafish (Danio rerio). Hear Res. 2006; 213(12):1–10.
 49
Le Goc G. ThermoMaster Lab. Gitlab:12407048. https://gitlab.com/GuillaumeLeGoc/thermomasterlab/.
 50
Tarantino N, Tinevez JY, Crowell EF, Boisson B, Henriques R, Mhlanga M, Agou F, Israël A, Laplantine E. TNF and IL1 exhibit distinct ubiquitin requirements for inducing NEMO–IKK supramolecular structures. J Cell Biol. 2014; 204(2):231–45.
 51
Gillespie DT. Exact numerical simulation of the OrnsteinUhlenbeck process and its integral. Phys Rev E. 1996; 54(2):2084–91.
Acknowledgements
We thank the IBPS fish facility staff for the fish maintenance, in particular Stéphane Tronche and Alex Bois. We are grateful to Carounagarane Dore for his contribution to the design of the experimental setup, Geoffrey Migault for engineering expertise, and Raphaël Voituriez for fruitful discussion on the modeling aspects.
Funding
This work was funded by the CNRS, Sorbonne Université, and the Systems Biology Network of Sorbonne Université and supported in part by the Human Frontier Science Program under grant no. RGP0060/2017, by the French Research National Agency under grant no. ANR16CE160017 and the European Research Council under the European Union’s Horizon 2020 research and innovation program grant agreement no. 715980. Funding bodies did not have any role in the design of the study and collection, analysis, and interpretation of data nor in writing the manuscript.
Author information
Affiliations
Contributions
G.L.G., V.B., R.C., and G.D. conceived the project. R.C. designed the setup. G.L.G. and J.L. performed the experiments. G.L.G., S.K., and G.D. analyzed the data. All authors contributed to the manuscript. The authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1
Figure S1: Correlations between parameters. A Pearson’s correlation coefficients between perbout parameters, reorientation angles of turn bouts, interbout interval and displacement. B Pearson’s correlation matrix between temperatureaveraged parameters. C Variance explained by the principal components of the intertemperature matrix. D Eigenvalues of the pooled intratemperature matrix. The red line highlights the KaiserGuttman criterion.
Additional file 2
Figure S2: PCA in singlefish experiments. A Variance explained by the five principal components for each singlefish. BC Autocorrelation function of the projection on PC1 (B) and PC2 (C) from each fish in singlefish experiments. The color code is the same as in A, black line and shaded area is the mean and s.e.m. across fish. D Mean variance of projections across time (intra, purple) and overall variance of projections (green). Error bars for intra is the s.e.m. and error bars for overall is 95 confidence intervals after bootstrapping (n=1000 boots).
Additional file 3
Figure S3: Temperatureindependant rescaling of parameters. A Equation describing parameter X distribution. BE Left to right, temperatureaveraged value, trajectoryaveraged rescaled by temperature averagedvalue and perbout value rescaled by the trajectory average, for B interbout intervals, C displacements, D reorientation angle of turn events, E turning probability.
Additional file 4
Figure S4: Fish position distributions along a linear thermal gradient. Presence probability density function of 10 batches of 10 larvae experiencing a thermal gradient from 18 ^{∘}C to 33 ^{∘}C. Solid line is the mean across batches, shaded area is the s.e.m. Dashed line is the expected value for a uniform distribution.
Additional file 5
Figure S5: Controls. A Example statistics of three batches successively recorded in E3 (control) and in E3 previously heated at 45 ^{∘}C for one hour in the experimental setup. (Left to right, up to bottom) Mean interbout interval, mean displacement, amplitude of turn bouts reorientation angle and mean square reorientation. B Example statistics from all experiments pooled together, dividing time into three windows of 10 min. (Left to right, up to bottom) Mean interbout interval, mean displacement, amplitude of turn bouts reorientation angle. C Example statistics at a given temperature (current T) as a function of the previously tested temperature. One dot is the mean parameter value for one experiment, color encoded current temperatures, lines are the mean for each previous temperature. (Left to right, up to bottom) Mean interbout interval, mean displacement, amplitude of turn bouts reorientation angle.
Additional file 6
Figure S6: Effect of cutoffs in trajectory selection. AB Example statistics when changing cutoffs in trajectory selection. (Left to right, up to bottom) Mean interbout interval, mean displacement, fraction of turns, amplitude of turn bouts reorientation angle, first principal component coefficients, second principal component coefficients. A With a minimum time between two consecutive bout of 200ms, trajectory must last at least 25 s. B With a minimum time between two consecutive bout of 400ms, trajectory must last at least 5 s.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Goc, G.L., Lafaye, J., Karpenko, S. et al. Thermal modulation of Zebrafish exploratory statistics reveals constraints on individual behavioral variability. BMC Biol 19, 208 (2021). https://doi.org/10.1186/s1291502101126w
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1291502101126w
Keywords
 Behavior
 Variability
 Thermokinesis
 Zebrafish
 Navigation
 Locomotion