Mathematical Modeling of Zebrafish Social Behavior in Response to Acute Caffeine Administration

Zebrafish is a model organism that is receiving considerable attention in preclinical research. Particularly important is the use of zebrafish in behavioral pharmacology, where a number of high-throughput experimental paradigms have been proposed to quantify the effect of psychoactive substances consequences on individual and social behavior. In an effort to assist experimental research and improve animal welfare, we propose a mathematical model for the social behavior of groups of zebrafish swimming in a shallow water tank in response to the administration of psychoactive compounds to select individuals. We specialize the mathematical model to caffeine, a popular anxiogenic compound. Each fish is assigned to a Markov chain that describes transitions between freezing and swimming. When swimming, zebrafish locomotion is modeled as a pair of coupled stochastic differential equations, describing the time evolution of the turn-rate and speed in response to caffeine administration. Comparison with experimental results demonstrates the accuracy of the model and its potential use in the design of in-silico experiments.


INTRODUCTION
Animal experiments are a standard practice for hypothesis testing in preclinical research [1,2]. However, experimental studies involving pharmacological treatment of sentient animals continue to raise ethical concerns regarding the well-being of the animals [3]. Computational methods can enable in-silico experiments that might help in the fulfillment of the 3Rs: Reducing the number of subjects, Refining experimental design and setup, and Replacing the use of live subjects [4][5][6].
Zebrafish (Danio rerio) has emerged as a species of choice in experimental studies in pharmacology where it is used in high throughput drug screening of several psychoactive compounds [7,8]. Its genetic and physiologic similarities with humans have made the zebrafish an attractive species for experimental investigations of human dysfunctional processes [9]. In particular, zebrafish experiments could clarify some of the open questions on anxiety-related behaviors in humans [10]. In these experiments, fish behavior is monitored in an experimental setup to investigate how anxiety-related behavior is modulated by anxiolytic and anxiogenic compounds, such as caffeine, cocaine, and ethanol [11][12][13][14][15]. Experiments on fish treated with such compounds have revealed numerous anxiety-related behaviors, including erratic activity (jump turns and sudden change in direction), thigmotaxis (tendency to stay near the wall), geotaxis (tendency to stay at the bottom of the tank), and freezing [16][17][18].
Previous efforts have leveraged data-driven, mathematical models to accurately describe the locomotion of isolated fish swimming in shallow or deep water tanks [19][20][21][22][23]. Concerning zebrafish, a number of studies have sought to incorporate their unique burst-and-coast swimming style, composed of sudden tail bursts that are followed by coasting phases [24,25]. The general line of approach consists of formulating a stochastic differential equation (SDE) for the turn-rate evolution, in which white noise is superimposed to intermittent excitation in the form of a jump process [26]. The original jump persistent turning walker (JPTW) was later adapted to the study of the effect of psychoactive manipulations in two separate studies [19,21]. Mwaffo and Porfiri [21] investigated the effect of acute ethanol treatment of zebrafish on model parameters of the JPTW, discovering a strong effect of concentration on the parameters of the jump process. Burbano-Lombana and Porfiri [19] expanded on the JPTW to simulate zebrafish response to acute caffeine administration. Not only did the model account for speed modulation during locomotion through an additional SDE, but also did it incorporate a detailed treatment of freezing episodes using discrete-time Markov chain. Overall, these studies suggest that the sensitivity of model parameters to the administration of psychoactive compounds must be considered when performing projective, in-silico experiments.
Other studies have extended individual fish models to groups, thereby including fish social interactions in terms of schooling and shoaling behaviors. In these models, social interaction is introduced as a response function that modulates the speed and turn-rate. Visual stimuli associated with the presence of conspecifics have been often considered in these models [23,[27][28][29][29][30][31][32], where fish tend to align and swim closer to neighboring subjects accommodating to alignment and attraction forces. Related efforts have included hydrodynamic interactions to incorporate lateral line sensing of the flow caused by neighboring subjects [34][35][36][37]. Overall, the mathematical underpinnings of these studies are common to the investigation of the structure of collective behavior of several species, from ants [38] to bats [39].
To the best of our knowledge, models looking at the effect of psychoactive compounds on zebrafish social behavior have never been explained in the literature. Here, we aimed to partially fill this gap by proposing a model that not only captures the effect of caffeine administration on fish locomotory activity but also takes into consideration the influence of the social environment in modulating the pharmacological response. To this end, we modeled fish dynamics in terms of speed and turn-rate, along two time-scales similar to Burbano-Lombana and Porfiri [19]. We defined a slow time-scale that captures the transitions between swimming and freezing states using a discrete-time Markov chain. During the swimming state, we modeled the speed and turn-rate evolution along a fast time-scale as a system of coupled SDEs. In the evolution of the turn-rate, we accounted for social interactions for each subject based on visual cues from neighboring individuals, therein, we utilized different interaction parameters depending on the treatment of the specific subject. To calibrate the model parameters, we relied on the experimental data-set from Neri et al. [40], wherein a group of untreated subjects swam with a caffeine-treated individual.
We investigated the value of the social group in modulating the response of fish to caffeine administration. Specifically, we compared calibrated model parameters for a treated fish swimming with an untreated group with those of a treated fish swimming in isolation from Burbano-Lombana and Porfiri [19]. We further highlighted an asymmetric interaction between the treated individual and untreated subjects, associated with the effect of caffeine on locomotory activity of fish and how it is perceived by untreated subjects. Lastly, we verified the predictive ability of the proposed model in capturing the social behavior of the group by comparing a set of social interaction metrics obtained from in-silico experiments to those from real experiments.

MATERIALS AND EQUIPMENT
Our theoretical endeavor is grounded in experiments from Neri et al. [40] (approved by the Animal Welfare Committee of New York University: protocol number 13-1424) on the effect of acute caffeine treatment on social behavior. Below, we summarize the main components of the experimental framework and data analysis from Neri et al. [40].

Experiment Setup and Procedure
The setup consisted of a circular tank of diameter d 90 cm filled with water at depth h 10 cm. Cameras were used to record fish behavior at 40 frames/s for 5 minutes (T exp 300 s). Videos were processed by an in-house multitarget tracking system developed in MATLAB [41,42].
Experiments were performed on groups of five adult subjects, including four untreated individuals and one treated individual, at four different caffeine concentrations: 0 (vehicle), 25, 50, and 70 mg/L. For each trial, five fish were randomly chosen from the holding tank. 50 fish were chosen at random to conduct ten experimental trials for each caffeine concentration (200 fish in total). One of the fish was kept in a 0.5 L beaker of a caffeine solution for 1 hour. Four untreated fish were introduced to the circular arena at the same time the beaker with the treated fish was placed in the arena. After 10 minutes of habituation, the treated fish was hand-netted from the beaker and released into the arena. The average fish body length (BL) was approximately 3 cm.

Data Post-processing
Fish trajectories were obtained by tracking the centroid of each fish. Figure 1 illustrates representative trajectories from each concentration. The trajectory of the i-th fish is denoted by (x i (kΔ), y i (kΔ)), where Δ 0.025 s is the sampling time, and k ∈ 1, . . . , K T exp Δ . Position increments between consecutive readings were used to obtain the velocity v i kΔ [v i,x kΔ , v i,y kΔ T and the speed v i kΔ v 2 i,x kΔ + v 2 i,y kΔ . To calculate the turn-rate, ω i (kΔ), we estimated the fish heading, θ i (kΔ), by fitting three consecutive positions, (x i ((k − 1)Δ), y i ((k − 1)Δ)), (x i (kΔ), y i (kΔ)), and (x i ((k + 1)Δ), y i ((k + 1)Δ)), along a circle [20]. The turn-rate was then inferred from the heading increment, δθ i (kΔ), between the two lines connecting the center of the circle with the (k − 1)-th and (k + 1)-th centroid position on the circle as ω i (kΔ) δθi(kΔ) 2Δ . Without loss of generality, we take i 1 as the treated fish throughout this paper.
Fish trajectories were also used to score the time spent freezing, an anxiety-related behavior in zebrafish [18]. According to Kopman et al. [43], a fish was considered to be in a freezing episode if it stayed within 2 cm radius for at least T F 2 s. From experimental data, we defined a binary Boolean variable Γ i (nT F ), with n ∈ 1, . . . , T exp T F that recorded instances of swimming (Γ i (nT F ) 1) and freezing (Γ i (nT F ) 0).
Four experimental trials were discarded due to recording issues (two from 0 mg/L, and two from 50 mg/L). We omitted four additional experimental trials due to insufficient data points for experimental analysis and parameter calibration (two from 25 mg/L, and two from 70 mg/L), whereby the fish spent less than 10 s in the swimming state and more than two BL away from the wall. For this reason, the experimental results presented in this paper may differ from those presented in Neri et al. [40] that relies on the same data-set.

METHODS
Here, we introduce the proposed data-driven framework to study the effect of caffeine treatment on individual and social behavior. Concerning our previous work [19], this study contributes a detailed model of social behavior, including attraction and alignment between subjects. Most importantly, these parameters are functions of the caffeine concentration and vary between treated and untreated subjects.
With respect to the state of the art on social behavior, the proposed model brings forward the critical role of the freezing response, by developing a two-time-scale modeling dichotomy where freezing evolves a slow time-scale that dictates when the animal is swimming or motionless. During locomotion, we used two coupled stochastic differential equations (SDEs) to describe the evolution of the turn-rate and the speed. The variables and notation used in the manuscript are included in Supplementary Table S1 in the supplemental material.

Zebrafish Kinematics
The fish were swimming in a shallow water tank, such that we could consider a two-dimensional (2D) model to describe their motion. Each fish was modeled as a rigid body, moving in a global reference frame [X, Y] with origin O. The position of the centroid of fish i at time t was denoted as [x i (t), y i (t)] T . We also measured the heading θ i (t) ∈ [ − π, π] as the angle between the swimming velocity and the global reference frame. Hence, the pose of fish i, that is, its position and orientation [44], was described as a three-dimensional vector [x i (t), y i (t), θ i (t)] T , as shown in Figure 2A. The evolution of zebrafish pose was modeled as a first-order kinematic model with initial conditions x i (0) x 0,i , y (0) y 0,i , and θ i (0) θ 0,i . Here, v i (t) and ω i (t) were the speed and turn-rate of the fish, respectively. We developed a mathematical model for the timeevolution of v i (t) and ω i (t) to predict the individual and social response of zebrafish.

Freezing Model
We adopted a discrete-time Markov chain to capture the transitions between freezing and swimming states. Building on the work of Burbano-Lombana and Porfiri [19] on isolated animals, for the i-th fish, we introduced a binary process Γ i (nT F ) that takes values 0 (freezing, F) and 1 (swimming, S), where n ∈ {1, . . . , Y}, Y Tsim TF , and T sim was the total simulation time. The Markov chain was determined by probabilities of persistence in swimming and freezing states, p S,i and p F,i , respectively, and probabilities of state transition, given by p SF,i 1 − p S,i and p FS,i 1 − p F,i , respectively.
The speed and turn-rate of the i-th fish were Such that during a freezing episode both the speed and turn-rate are zero and during swimming they evolve based on the SDEs described below.

Locomotion and Interaction Models
Speed and turn-rate in the swimming state were modeled as a system of two coupled SDEs. In the model, we included social interaction terms that modulate fish locomotion based on the visual cues from neighboring conspecifics. As illustrated in Figure 2B, we described fish schooling between the focal fish, i, and the neighboring fish, j, in terms of the relative orientation, Further, we examined fish shoaling in terms of the relative position of the neighboring fish with respect to the focal fish expressed in terms of the distance between the pair of fish, s ij (t), and relative angle, θ ij (t).
To model the evolution of the speed, we adopted the following logistic model, similar to Burbano-Lombana and Porfiri [19] for a single fish [45]: where η i [s −1 ] and σ v,i [s − 1 2 ] were the linear expansion rate and the strength of the added noise, respectively; W v,i (t) was a standard Wiener process; and g(ω S,i (t)) [m −1 ] encapsulated the effect of the turn-rate. Specifically, the speed response function was where std ω,i was the standard deviation of the absolute instantaneous value of the turn-rate [19]. This function captured the need of fish to slow down when turning, while attaining larger speeds during straight swimming. This model offers a first approximation of speed modulation during social behavior. For each fish, the model required the calibration of two parameters, assuming that the body length is common to the entire group: η i , and σ v,i . In this basic incarnation, the model does not incorporate speed-based social interaction, which has been proposed by several authors to play some role in the social response of social fish [46][47][48][49][50]. The choice of neglecting social interactions mediated by the speed was due to the need of reducing the number of model parameters, magnified by the presence of individual differences in the treatment of the group.
The turn-rate dynamics were captured by the JPTW [26,33], where ω * S,i (t) [rad s −1 ] was the turn-rate interaction response function; f w (ϕ w,i (t), d w,i (t)) was the wall interaction function where ϕ w,i (t) was the projected angle to collision and d w,i (t) was the distance from the wall; α i [s −1 ] was a positive parameter quantifying the relaxation rate; σ ω,i [rad s − 3 2 ] was the strength of the added noise; W ω,i (t) was a standard Wiener process; and J i (t) was the jump noise term encapsulating sudden changes in the turn-rate.
Due to the presence of the caffeine treatment, the social interaction gains should vary in the group. We expected untreated fish to respond differently to a treated subject compared to an untreated one. Further, we anticipated the influence of treated on untreated to be different from the influence of untreated on treated, such that their interaction should be asymmetric. These claims are grounded in two propositions from the literature. First, the anxiogenic value of caffeine has been shown to influence the tendency of the caffeinetreated fish to interact with untreated conspecifics [15,51]. Second, the psychostimulatory nature of caffeine is known to influence the locomotory response of the animals [52], which may underlie differences in the appraisal of treated fish by untreated individuals. Accordingly, the turn-rate response function was written as where k p,ij [rad m −1 s −1 ] and k v,ij [rad m −1 ] were the attraction and alignment gains of fish i toward fish j, respectively. For each trial, the model required calibrating 2N − 1 pairs of gains.
T , swimming at speed v i (t) and turn-rate ω i (t) (B) A close-up look at the interaction between a pair of fish within a group of five fish. Alignment and attraction between the i-and j-th fish are functions of the relative orientation, ϕ ij (t), and relative position, in terms of the distance between fish, s ij (t), and relative angle, θ ij (t).
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org October 2021 | Volume 7 | Article 751351 We categorized these parameters based on the direction of interaction as follows: where TU, UT, and UU identified the response of the treated to untreated fish, the untreated to the treated fish, and the interaction between untreated subjects, respectively. The presence of Γ j (nT F ) in Eq. 6 was used to selectively limit the social response of fish to the group members that were actively swimming. Fish that were freezing were excluded from the social interaction model, based on calibration of the model on real data as well as biological observations that suggested zebrafish are more responsive to dynamic, rather than static stimuli [53]. The wall interaction function was written as follows [19,20]: where the intensity of wall interactions, a w [rad s −1 ], and the sensitivity of the fish to visual stimulus to the wall, b w [cm −1 ], were two positive parameters. We hypothesized that all fish interact in the same way with the environment, such that the two parameters a w and b w were the same for the entire group and for every trial. The selection of the form in Eq. 8 encapsulated the wall avoidance behavior of the fish and ensured that fish remain within the boundary of the tank; this selection did not capture wall-following behavior.
We finally modeled the jump noise for the i-th fish as a compounded Poisson process, Here, A k,i (t)'s were independent and identically distributed Gaussian random variables with zero mean and variance , and the total number of jumps at time t, m i (t), was such that its increments are Poisson random variables with parameter λ i (t′′ − t′) for time t′, t′′ and t′′ > t′, with λ i [s −1 ] being frequency of jumps.

Model Calibration
For each fish in the group, i ∈ {1, . . . ,N}, we calibrated the set of locomotion and social interaction model parameters. The transition probabilities for the discrete-time Markov chain model were obtained by simply counting instances of freezing and transitions to swimming in the experimental time-series. On the other hand, maximum likelihood estimation was applied to calibrate the locomotion model parameters.
In summary, we calibrated the following parameters: transition probabilities, p FS,i and p SF,i ; linear expansion rate, η i ; strength of added noise on speed, σ v,i ; relaxation rate, α i ; strength of added noise on turn-rate, σ ω,i ; intensity of jump turns, c i ; frequency of jump turns, λ i ; alignment gains of treated to untreated fish, k v TU , untreated to treated fish, k v UT,i , and between untreated fish, k v UU,i ; attraction gains of treated to untreated fish, k p TU , untreated to treated fish, k p UT,i , and between untreated fish, k pUU,i . Given that five fish comprised each of the groups, a total of 58 parameters were calibrated per trial.

Calibration of the Discrete-Time Markov Model for Freezing
We obtained the binary sequences {Γ i (nT F )} Y n 1 } from the experimental time-series for each fish in the group. Similar to Burbano-Lombana and Porfiri [19], we estimated the transition probabilities as follows: where N SF,i and N FS,i were the number of transitions by the i-th fish from swimming to freezing and from freezing to swimming, respectively. N SS,i and N FF,i were the number of instances in which the fish maintained the swimming or freezing state, respectively. Estimated transition probabilities for the treated fish in the group are shown in Supplementary Table S2. For completeness, in Supplementary Table S3, we also report a summary of the transition probabilities for the discrete-time Markov chain of the untreated fish in terms of mean and standard deviation calculated across all trials.

Calibration of the Locomotion and Interaction Models Through Maximum-likelihood Estimation
Using the experimental sampling time Δ as the time-step for discretization, we approximated Eqs 3, 5 via the Euler-Maruyama method as follows [54]: where ϵ vi (k) was a standard Gaussian random variable, utilized to approximate the added noise. We followed the same discretization approach to approximate the JPTW in Eq. 5, leading to.
where ϵ 1 ωi (k) and ϵ 2 ωi (k) were standard Gaussian random variables and ζ i (k) was a Bernoulli process with a probability Δλ i . Wall interaction was not included in the approximation of the JPTW in Eq. 12 since we performed calibration only when the fish were more than 2 BL away from the wall. For each individual, we consolidated unknown parameters in two vectors, φ v,i and φ ω,i , one for the speed and the other for the turn-rate dynamics, in Eqs 11, 12, respectively. These vectors were determined by solving two independent optimization problems for the speed and turn-rate. The parameters were estimated for each fish in the group independently for every trial.
For the approximated logistic equation in Eq. 11, the vector of unknown parameters for each fish was where we used a scaling factor, κ, to avoid singularities at near zero swimming speed [19]. The search was conducted within a set of admissible values χ v selected from previous work [22]. The optimization problem was solved by using as input the K * i samples of the speed obtained by excluding instances of freezing or swimming in proximity of the wall.
The maximum-likelihood estimation problem was formulated aŝ The likelihood function, l v,i (φ v,i , v S,i (kΔ), ω S,i (kΔ)), was derived from the model approximation in Eq. 11 as where H (x, σ) was the Gaussian distribution at x with zero mean and variance σ 2 . Further, q i (kΔ) was given by Heuristically, we found that κ 5 guarantees convergence of the optimization problem. A similar approach was adopted to calibrate the discrete JPTW in Eq. 12. For each fish, we solved the optimization problem for the vector of unknown parameters for each fish, φ ω,i [α i , σ ω,i , c i , λ i , k pij , k vij ] T , with j ∈ {1, . . . , N}, j ≠ i, where the interaction gains were categorized in accordance with Eq. 7. We used an input of K * i samples of the turn-rate obtained by excluding instances of freezing or swimming close to the wall. In addition, the search was conducted within a set of admissible values χ ω selected from Butail et al. [27] and Mwaffo et al. [22]. The maximum-likelihood estimation problem was formulated aŝ where χ ω was in R 6 for the treated fish (i 1) and χ ω was in R 8 for the untreated fish (i ≠ 1). The likelihood function l ω,i (φ ω,i , v S,i (kΔ), ω S,i (kΔ)) was defined as and z i (kΔ) was The locomotion parameters of each treated fish for all trials are displayed in Supplementary Table S4. A summary of the parameters of the untreated fish is shown in Supplementary  Table S5, in terms of mean and standard deviation calculated across all trials. Supplementary Table S6 displays the attraction gains of the treated fish k p TU , and the attraction gains of the untreated subject towards treated neighborsk p UT and untreated neighborsk p UU where a hat denotes the average of untreated individuals in each trial. Similarly, Supplementary Table S7 contains the alignment gains of the treated fish k v TU , and the alignment gains of untreated subjects towards treated neighborsk v UT and untreated neighborsk v UU . We discarded two additional trials from 25 mg/L and one additional trial from 50 mg/L due to divergence of the estimator, where interaction gains converged to their upper bounds.

Calibration of the Wall Function
We relied on the work of Burbano-Lombana and Porfiri [19] to obtain the wall function parameters in Eq. 8. The wall interaction function was calibrated for a fish swimming alone, from the dataset of Neri et al. [40], using a wall-corrected turn-rate from the real time-series of the turn-rate of fish swimming alone, where ω a (kΔ) was the turn-rate of the fish swimming alone and ω c (kΔ) was the corrected turn-rate. Next, ω c (kΔ) was plotted against the distance from the wall d w (kΔ) where only the positive values of the corrected turn-rate were considered to capture wall avoidance. A robust non-parametric locally weighted least squares (RLOESS) function in MATLAB was used to fit the signal to a parametric exponential function. As such, the wall interaction parameters were obtained by calculating the average across all trials as a w 11.68 rad s −2 and b w 0.19 cm −1 .

RESULTS
We began our analysis of the model by examining the influence of caffeine concentration on fish locomotion in terms of the variations of relevant model parameters. With respect to parameters on freezing response and locomotion, we compared with model parameters obtained in Burbano-Lombana and Porfiri [19] to assess the effect of the social environment on fish response to caffeine administration. Finally, we conducted in-silico experiments to demonstrate the predictive power of the model in anticipating experimental results on schooling and shoaling.

Analysis of Model Parameters
First, we investigated the effect of caffeine concentration and social environment on the locomotion parameters of the treated fish, utilizing two-way ANOVA with caffeine concentration and social environment (single or group) as independent variables. Second, we conducted ANOVA comparisons with caffeine concentration as a single independent variable to compare the interaction parameters across concentrations. Post-hoc comparisons were conducted using Tukey's HSD (honestly significant difference). The significance level was set to 0.050 throughout. We found that caffeine concentration did not influence the Markov chain transition probabilities p FS (F 3,50 0.424, p 0.738) and p SF (F 3,50 0.125, p 0.944), neither in isolation nor in group (shown in Figures 3A,B, respectively). No difference was detected across social environment with respect to p FS (F 1,50 0.630, p 0.443). Although we registered a dependence on the social environment with respect to p SF (F 1,50 5.416, p 0.027), we did not detect any variation in post-hoc analysis. The interaction between the two independent variables was found to be not significant for both p FS (F 3,50 1.733, p 0.181) and p SF (F 3,50 0.812, p 0.497).
Likewise, the linear expansion rate, η, was not influenced by either caffeine concentration (F 3,50 1.264, p 0.297) or social environment (F 1,50 0.698, p 0.407), shown in Figure 4A. Further we did not detect differences in the interaction of the independent variables on η (F 3,50 0.048, p 0.986). In terms of the strength of added noise on the speed evolution, σ v , we found a dependence on caffeine concentration (F 3,50 3.039, p 0.038; Figure 4B), which, however was not accompanied by variations in post-hoc analysis. We found that the presence of the social environment had an effect on σ v (F 1,50 33.21, p < 0.001), and post-hoc analysis indicated a decrease in the strength of added noise in the presence of untreated subjects for 0 mg/L. We did not detect a significant interaction between the independent variables on σ v (F 3,50 1.088, p 0.363).
With respect to the turn-rate model parameters, we did not detect an effect of caffeine concentration on the mean reversion rate, α (F 3,50 1.368, p 0.263). Although we found α to be affected by the social environment (F 3, 50 15.49, p < 0.001; Figure 5A), post-hoc analysis did not reveal significant differences between concentrations. Likewise, we did not discover a significant interaction between caffeine concentration and social environment on α (F 3,50 0.519, p 0.672). While caffeine concentration was found to have an influence on the strength of added noise in the turn-rate evolution, σ ω (F 3, 50 2.926, p 0.043; Figure 5B), no variations were identified in post-hoc analysis. We determined a modulatory role of the social environment (F 3,50 24.83, p < 0.001), where σ ω increased in the presence of a social group for 50 mg/L in post-hoc analysis. No significant interaction was detected between the independent variables with respect to σ ω (F 3,50 0.866, p 0.465). With respect to intensity of jumps, c, we found caffeine concentration to play a modulatory role (F 3,50 5.760, p 0.002; Figure 5C), with post-hoc analysis revealing a decrease in the intensity of jumps for treated fish swimming in isolation from 50 to 70 mg/L. In addition, we found the social environment to influence c (F 1, 50 15.90, p < 0.001), where we detected an increase in the jump intensity in the presence of FIGURE 3 | Comparisons of discrete-time Markov chain parameters of the treated fish across caffeine concentrations and social environment (single or group). The bars represent the mean value of the probability of transition from freezing to swimming (A), and the mean value of the probability of transition from swimming to freezing (B). The striped bars correspond to the calibrated parameters for the case of a single treated fish from Burbano-Lombana and Porfiri [19]. The solid bars correspond to the calibrated parameters for the case of a treated fish swimming in a social group. The vertical red error bars represent standard errors of the means. untreated subjects for 0 mg/L in post-hoc analysis. We did not identify a significant interaction between caffeine concentration and social environment with respect to c (F 3,50 0.747, p 0.529). Finally, the frequency of jumps, λ, was not affected by caffeine concentration (F 3,50 2.166, p 0.104). In contrast, we detected significant differences across social environment (F 1,50 13.65, p < 0.001; Figure 5D). Post-hoc analysis revealed that fish swimming in isolation had higher values of λ than those swimming in group for the 25 mg/L treatment. We registered a significant interaction of the independent variables on λ (F 3,50 2.924, p 0.048). Next, we investigated the effect of caffeine concentration on the interaction gains in the turn-rate model, as shown in Figure 6. We identified an effect of caffeine concentration on the attraction gain of the treated fish towards untreated fish, k p TU (F 3,22 3.323, p 0.038), but post-hoc analysis did not detect differences between concentrations. The average attraction gain,k p UT , of the untreated fish towards treated fish was not found to vary with caffeine concentration (F 3,22 0.588, p 0.629). We determined that the average attraction gain of the untreated fish towards other untreated subjects,k pUU , varied with caffeine concentration (F 3,22 3.679, p 0.028), and post-hoc analysis brought to light a decrease from 0 to 25 mg/L. Finally, the alignment gains were indistinguishable with respect to caffeine concentration: k v TU (F 3,22 1.252, p 0.315),k v UT (F 3,22 0.756, p 0.531), andk v UU (F 3,22 0.596, p 0.459).
In summary, among all the freezing and locomotion parameters, we only found the intensity of jumps to depend on caffeine concentration, yet, without differences with respect to vehicle-treated individuals. Comparisons across social environment revealed variations in the strength of added noise on both speed and turn-rate and in the jump parameters. Swimming in group reduced the strength of the added noise on the speed evolution of vehicle-treated subjects, and it increased the strength of the added noise on the turn-rate evolution at the intermediate concentration. Further, while the presence of a social group increased the intensity of jumps of vehicle-treated subjects, it reduced the frequency of jumps of individuals treated at a low concentration. Parameters pertaining to social response were generally robust with respect to caffeine concentration, except for the attraction of untreated fish towards other untreated subjects, with low caffeine concentration causing a reduction in alignment.

In-silico Experiments
We conducted in-silico experiments to validate the developed model and investigate its ability to predict the social behavior of fish detected from experimental time-series [40], for a range of interaction metrics that quantify schooling, and shoaling.
Schooling is a measure of fish tendency to align their bodies during swimming [55,56]. The degree of alignment among the four untreated fish was scored in terms of the instantaneous polarization [57], where N 5 was the number of fish in the experiment. Polarization varies between 0 and 1, where 1 identified the case in which untreated fish are perfectly aligned in the same direction. The alignment between the treated fish and the untreated group of fish was scored in terms of the relative instantaneous polarization, R (kΔ), Relative polarization ranges between −1 and 1, where 1 corresponded to the group of untreated fish pointing in the same direction of the treated fish, and −1 indicated that the treated fish is pointing in the opposite direction to the group of untreated fish. These quantities were averaged in time to compute the average polarization and the average relative polarization.
To quantify fish shoaling, the tendency of fish to swim nearby, we computed the inter-individual distance, d ij (kΔ), between each pair in the group. We scored the average distance between the treated and untreated subjects, and the average distance among untreated individuals.
We conducted in-silico experiments using the model parameters for the case of in-group swimming, shown in solid bars in Figures 3-6. Ten simulations were performed for each of the four caffeine concentrations. For each fish in all 40 trials, the interaction gains were sampled from a Gaussian distribution with mean and standard deviation of the corresponding parameter at that concentration. On the other hand, since we did not find any effect of caffeine concentration on the transition probabilities and locomotion parameters for in-group swimming, those parameters were taken as the average of all fish across all experimental trials based on treatment. The initial conditions x i (0), y i (0), θ i (0), Γ i (0), v i (0), and ω i (0) were chosen uniformly at random in their respective intervals. Time-series of four trajectories for each caffeine concentration are shown in Figure 7; videos are presented at the link provided in the data availability statement. Note that the wall function adopted in this study does not take into consideration the wall following behavior of zebrafish, thus explaining the differences in wall interactions between in-silico trajectories in Figure 7 and experimental ones in Figure 1.
We performed statistical analysis to compare the social interaction metrics across caffeine concentrations, and validate the in-silico results against those obtained from real experiments. For this purpose, we conducted two-way ANOVA with caffeine concentration and data-type (experiments or in-silico) as independent variables. Post-hoc comparisons were conducted using Tukey's HSD. The significance level was set to 0.050 throughout.
We detected influence of caffeine concentration on the average polarization, P (F 3,61 7.781, p < 0.001), shown in Figure 8A. Post-hoc analysis revealed differences in experimental results, where the average polarization was found to increase from 0 to 50 mg/L. Comparisons across data-types did not indicate differences between real and in-silico experiments (F 3,61 1.354, p 0.249). Likewise, no interaction between the independent variables was identified on P (F 3,61 0.675, p 0.571). With respect to the average relative polarization, R ( Figure 8B), we did not find an effect on either caffeine concentration (F 3,61 1.354, p 0.071) or data-type (F 1,61 0.229, p 0.634), although we identified significant interaction (F 3,61 3.855, p 0.014).
Next, we determined that the shoaling tendency between the treated fish and untreated subjects, in terms of the average distance d T−U , was consistent across caffeine concentrations (F 1,61 1.849, p 0.179) and data-types (F 1,61 1.461, p 0.234), as shown in Figure 9A. No interaction was detected between the independent variables on d T−U (F 3,61 0.262, p 0.853). In contrast, we detected an effect of caffeine concentration on the average distance between the untreated fish, d U−U (F 3,61 12.16, p < 0.001; Figure 9B). Post-hoc analysis revealed a decrease in d U−U from 25 to 50 mg/L in the experimental data-set. Similar differences were found in the in-silico data-set where d U−U was larger for 25 mg/L than 0 and 50 mg/L. While comparisons between data-types revealed a significant difference (F 1,61 11.29, p 0.0.001), the results were indistinguishable between real and in-silico experiments in post-hoc analysis. Finally, we did not identify a significant interaction between the independent variables on d U−U (F 3,61 2.354, p 0.081).

DISCUSSION
In this work, we developed a modeling framework to study the effect of acute caffeine treatment on the social behavior of zebrafish. We contributed two key advances to previous work on modeling collective behavior of zebrafish. First, similar to the analysis with respect to zebrafish swimming alone in Burbano-Lombana and Porfiri [19], we included the freezing response of each individual within the group, which is necessary to capture anxiety-related behavior in response to caffeine [58]. For this purpose, we developed a two-time-scale modeling dichotomy. Along a slow time-scale, we used a discrete-time Markov chain to describe the transition between swimming and freezing states. At a fast time-scale, we modeled the evolution of the speed and turnrate during swimming as a system of coupled SDEs: a logistic equation to represent the speed and a JPTW to describe the turnrate. Second, we granularly tracked the directional interaction between each pair of fish based on the treatment of each fish within the pair. This approach takes into consideration previous experimental work highlighting the effect of caffeine on the behavioral response of treated fish and its appraisal by untreated conspecifics [51,52]. We calibrated the model on real experimental data from previous work [40], where we studied groups of caffeinetreated subject and untreated individuals swimming in a shallow water tank. For each group of five individuals, we estimated 20 parameters, entering the Markov chain and the SDEs. Calibration employed a combination of maximum likelihood estimation and classical plug-in estimation. We displayed our results on two fronts. First, we compared the model parameters obtained for a treated fish swimming with untreated subjects with those obtained by Burbano-Lombana and Porfiri [19] for the case of an isolated fish. Second, we compared the social interaction metrics, in terms of average polarization, average relative polarization, and average inter-individual distance, between the experimental and in-silico data-sets.
In contrast with our expectations, we did not observe a modulatory role of caffeine on freezing and locomotion parameters. Our expectation was based on a number of previous studies documenting a robust dependence of zebrafish behavioral response to acute caffeine administration [18,59]. This was particularly evident for individuals swimming with untreated subjects, for which we failed to detect any effect of caffeine treatment. Likely, the explanation for the abolishment of a dose-dependent response should be sought in the presence of the social environment, which, indeed was responsible for a few, salient variations in locomotion parameters associated with the speed and turn-rate evolution. It is tenable that the presence of social cues had a leveling role on the anxiogenic effect of caffeine, which is indirectly evidenced by the tendency to enhance white noise with respect to the jump noise in the turn-rate evolution. Jumps have been associated with erratic activity of the animal, in the form of C-and U-turns, so that their reduction in favor of steady swimming offers an indication of an axiolytic value of the social environment, also discussed by Neri et al. [40]. With respect to the effect on untreated subjects, we recorded a decrease in their tendency to shoal with each other, which highlights an interesting, albeit indirect, effect of caffeine treatment. Caffeine treatment of one selected individuals might bear an anxiolytic effect on the rest of the group that reduce their tendency to stay close [15,51]; understanding this counter-intuitive finding should be the object of future research.
The calibrated model is in good agreement with experimental observations on social metrics, related to shoaling and schooling. While this agreement should be desired in any calibrated model, it is not obvious to attain. In fact, in-silico experiments do not contain the fine-grain variations that are unique to the experimental subjects, whereby we excluded from the simulations any statistical variation in the locomotion and freezing parameters. Accounting for variations in the social gains due to caffeine administration through a simple normal distribution seems sufficient to capture the emergent response of the groups, as well as the role of the treated individual.
The proposed model is not free of limitations. First, we assumed that the interaction between fish is solely based on visual stimuli. Incorporating other mechanisms of social interactions, such as hydro interactions [35,37], may help refine the mathematical model, especially in terms of short-range interactions related to the perturbations they create in the fluid environment [37]. Second, the current model does not incorporate wall following behavior observed in real experiments, whereby interaction with the wall is limited to a simple repulsion [20]. Third, the model is purely twodimensional, thereby failing to capture salient anxiety-related responses that have been documented in zebrafish, such geotaxis [18,60]. Fourth, the entire modeling framework is based on a single psychoactive compound, which bears limitations in the generalizations of the predictions to other substances that impinge on anxiety [11,14]. Along this line, the most fundamental limitations of the model is the lack of a direct link between the molecular composition of the substance or the brain mechanisms it affects and the parameters of the model. In its present incarnation, the model requires knowledge of all its parameters to perform in-silico experiments, without allowing for exploring different substances or even untested concentrations on caffeine.
Despite these limitations, the proposed model offers a first step in the design of in-silico experiments that can aid the 3R's with respect to zebrafish experimentation. The proposed model can be FIGURE 9 | Comparisons of the shoaling tendency of the fish, measured in terms of the average distance between treated fish and untreated fish (A), and average distance between untreated fish (B), across caffeine concentrations and data-types (experiment or in-silico). Different letters on top of the bars indicate a significant difference (p < 0.050) in Tukey's HSD post-hoc analysis across caffeine concentrations, comparing interaction metrics in experiment (standard font) or in-silico (italic font) data-type. The vertical red error bars represent standard errors of the means.
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org October 2021 | Volume 7 | Article 751351 used to reduce the number of experiments, by affording statistical insight into the sample size. Likewise, the model can be used to refine existing data-sets, by informing model-based analysis of the data and, potentially, assist in verification and tracking. Finally, pilot studies could be conducted on a computer, thereby reducing the number of subjects utilized in experimental research.

DATA AVAILABILITY STATEMENT
The experimental data-set and representative videos of the insilico experiments can be found in the repository of our research laboratory at https://github.com/dynamicalsystemslaboratory/ CaffeineSocialBehavior.

AUTHOR CONTRIBUTIONS
MP designed and supervised the research. MT developed the computer codes, conducted model calibration, performed statistical analysis, and wrote a first, preliminary draft of the