Using neuronal models to capture burst and glide motion and leadership in fish

Linnéa Gyllingberg Department of Mathematics, Uppsala University, Uppsala, Sweden Alex Szorkovszky RITMO Centre for Interdisciplinary Studies in Rhythm, Time and Motion, University of Oslo, Oslo, Norway David J. T. Sumpter Department of Information Technology, Uppsala University, Uppsala, Sweden

Abstract

While mathematical models, in particular self-propelled particle (SPP) models, capture many of the observed properties of large fish schools, they do not always capture the interactions of smaller shoals. Nor do these models tend to account for the observation that, when swimming alone or in smaller groups, many species of fish use intermittent locomotion, often referred to as burst and coast or burst and glide. Recent empirical studies have suggested that burst and glide movement is indeed pivotal to the social interactions of individual fish. In this paper, we propose a model of social burst and glide motion by combining a well-studied model of neuronal dynamics, the FitzHugh-Nagumo model, with a model of fish motion. We begin by showing that the model can capture the motion of a single fish swimming down a channel. By then extending to a two fish model, where visual stimuli of the position of the other fish affect the internal burst or glide state of the fish, we find that our model captures a rich set of swimming dynamics found in many species of fish. These include: leader-follower behaviour; periodic changes in leadership; apparently random (i.e. chaotic) leadership change; and pendulum-like tit-for-tat turn taking. Unlike SPP models, which assume that fish move at a constant speed, the model produces realistic motion of individual fish. Moreover, unlike previous studies where a random component is used for leadership switching to occur, we show that leadership switching, both periodic and chaotic, can be the result of a deterministic interaction. We give several empirically testable predictions on how fish interact and discuss our results in light of recently established correlations between fish locomotion and brain activity.

1 Introduction

A wide range of mathematical and computational models have been proposed to capture the collective motion of fish schools [1, 2, 3, 4, 5, 6, 7]. These models are often referred to collectively as self-propelled particle models (SPPs), with the name coming from a seminal paper by Vicsek et al. in 1995 [8]. The Vicsek model assumes that each fish (or particle) moves with constant speed, while its direction is updated at each time step to be closer to the average direction of individuals within its neighbourhood. A noise term is added to model uncertainty or error in the fish’s direction. Variations of this canonical model extend it to include repulsion, attraction and other social interactions [9, 10, 11, 12].

These models reproduce many of the observed properties of large fish schools [13, 12], but don’t always capture the interactions of smaller shoals. One difference is that the assumption that the fish moves with constant speed is not typical of the swimming behaviour for many species of fish. For example, zebrafish [14], koi carps [15], guppies [16], cod [17], red nose tetra fish [18] and many other species swim by alternating between accelerated motion and powerless gliding [19]. Even other animals, such as spiders, beetles and lizards have his type of intermittent locomotion [20]. For swimming animals, intermittent motion is often referred to as burst and coast or burst and glide. We use the latter term in what follows. Models of burst and glide behaviour have been actively developed in the biomechanics community for decades [19, 17, 21, 22, 23, 24, 25, 26]. It has been shown that there are energetic advantages of such movements, when compared to constant swimming speed. The energetic cost of swimming is minimized during the glide phase, where the body is rigid and the fish decelerate due to water resistance [19]. Once a low velocity is reached, the burst phase starts again and the fish accelerates to a maximum velocity. The cycle is then completed and the fish glides again.

Burst and glide movement is pivotal to detect and quantify social interactions between individual fish [16]. Herbert-Read at al. observed that, when swimming in pairs, guppies respond to the other fish by accelerating their motion when the other fish is nearby. Pairs of fish do not update their speed continuously, but at discrete time points, making Vicsek type of models unsuitable for understanding the social dynamics of small shoals of guppies. Fish living in high predation areas have more pronounced acceleration and deceleration responses [16]. Moreover, Kotrschal et al. found that the high burst speed in response to neighbours evolves when subjected to artificial selection [27].

Another aspect which is not included in the most common self-propelled particle models, but is found in fish schools of several species, is leadership. Schaerf et al. showed that for pairs of freely exploring eastern mosquitofish, it is possible to categorize the fish into leaders and followers — some fish takes the role of leaders and will on average spend more time in front than behind the other fish [28]. In this case, one particular fish is the leader. However, other studies show that fish change leader-follower roles. For example, Nakayama et al. show that leadership in pairs of stickleback changes over time [29] and Milinksi show tit-for-tat like strategies in pairs of swimming sticklebacks when confronting a potential predator [30].

As we outline above, even though studies of interacting burst and glide behaviour are becoming more common, most studies of the burst and glide responses are of isolated fish, with a focus on understanding the energetics. Even though some self-propelled particle models assume variable speed on the part of the fish [31], they do not necessarily model burst and glide behaviour. Cavoli et al. made an interesting start in this direction by incorporating burst and glide behaviour into an SPP model [32]. They have developed a model of two fish dynamics where the speed of individual fish is updated at discrete time points and the heading angle is updated due to attraction and alignment with the other fish. However, this model does not capture the observations from Herbert-Read et al. and Kotrschal et al. that there is a direct response in speed, and not just heading angle, when two fish interact with each other [16, 27]. Thus, developing a simple model for speed response of fish interactions remains an open question.

In order to incorporate social information into burst-and-glide dynamics, it is necessary to consider how burst timings and amplitudes are governed in the nervous system. Continuous rhythmic patterns of movement in fish and other vertebrates are achieved by a neural circuit known as a central pattern generator, with afferent feedback as an important modulating component [33]. Recent work indicates that the control of episodic bursting in larval zebrafish is similarly distributed, yet separable from the system that coordinates the bursts themselves [34]. Although the exact mechanism for this control is still unclear, there is ample evidence for an internal burst/glide state that is subject to modulation by sensory input.

There is also an intriguing similarity between the burst and glide of fish and the burst and recover phases of neurons. The FitzHugh-Nagumo model [35, 36], a well-known model of neuronal firing, is a simplification of the Hodgkin and Huxley’s Nobel prize winning model for action potentials in neurons [37]. The FitzHugh-Nagumo model has two variables, $V$ and $W$ , which do not explicitly relate to specific chemical and electrical properties of neurons. Instead, $V$ models the membrane potential and $W$ is thought of as a recovery variable. The neuron dynamics in the FitzHugh-Nagumo model is given by the following equations:

		$\displaystyle\frac{dV}{dt}=V-\frac{V^{3}}{3}-W+c$		(1)
		$\displaystyle\tau\frac{dW}{dt}=V+a-bW.$		(2)

The cubic-shaped dependence in the first equation captures the essential bistability of the membrane potential $V$ and, for certain parameter values, the system displays intermittent switching between firing and recovery, at a rate governed partially by the slow exponential decay of the recovery variable $W$ . The constant $c$ represents a tonic (i.e. constant) input to the neuron.

An interesting property of the FitzHugh-Nagumo model, is that when two FitzHugh-Nagumo model neurons are coupled, the system displays chaotic dynamics [38]. The evidence for chaos in neural mechanisms is extensive [39, 40, 41]. Chaotic dynamics occur on both macroscopic and microscopic scales in brain dynamics [42] and in both humans and animals [43, 44]. Neuroscientific studies show that chaos is also essential for many neural mechanisms. For example, Ohgi et al. show that chaotic dynamics are fundamental in the learning and control of the dynamical interactions between brain, body and environment in human infants [45].

There are currently no existing self-propelled particle models with internal neural states explicitly defined, despite sensory information being increasingly a focus for collective motion studies [46, 47, 48]. Instead, models tend to be based on heuristic “rules of interaction” where the ability to turn sensory information into actions, such as averaging neighbour headings to obtain a turning angle, is assumed to be instantaneous and left out of the model. More complete biological models, on the other hand, are difficult to scale up to multiple agents, and the high complexity makes it difficult to connect specific mechanisms to collective outcomes. However, there are intriguing starts in this direction made by O’Keeffe et al., where an SPP model is combined with the Kuramoto model [49].

Less direct analogies have previously been made between collective motion and activities in the brain, for example Marshall et al. draws parallels between decision-making in ants and decision-making in primate brains [50], Passino et al. shows that some main characteristics of cognition in brains of vertebrates are also present in swarms of honey bees [51], and on a more general level, Riberio et al. discusses the connections between animal group dynamics and networks in brains [52]. And there are existing models connecting sensory systems to neuronal dynamics. For example, Kuniyoshi et al. have presented a mathematical model coupling the sensory system in muscle spindles to neuronal dynamics in human infants [53]. Furthermore, there are indications of correlations between different swimming activities and activities in fish brains [54, 55], although as yet there is no proven link between burst and glide and neuronal bursting.

In this article, we consider the possibility that even simple deterministic neural models are rich enough in dynamics to capture seemingly stochastic burst-and-glide motion patterns in fish, and that this can be used as a basis for self-propelled particle modelling. We begin with a two variable dynamical model, inspired by the FitzHugh-Nagumo model, one representing the internal state and the other representing speed, as a model for a single fish. We then introduce social interaction, and characterize the behaviour of a pair of fish as a function of key parameters. The aim of our study is primarily to investigate the mechanisms at work [56, 57, 58, 59, 60], rather than to compare the model in detail to data. Nonetheless, we show that the patterns of interactions created in the model mimic those observed in single and pairs of guppies. We discuss several testable hypotheses arising from our model.

2 Single fish model

We modify the Fitzhugh-Nagumo model so that recovery is governed by sensory feedback of the fish speed rather than internal dynamics. Our model has two variables: the speed of the fish, $v$ , and the burst potential, $b$ , which indicates whether the fish is bursting or gliding. We model motion along a line, so that $v(t)$ defines the motion entirely. The velocity itself is subject to propulsive force and drag:

\displaystyle\frac{dv}{dt}=g(b)-kv\,.

(3)

Here $g(b)\geq 0$ is the propulsive force, which is a function of the internal state. To approximate the bursting impulse, we assume this function increases in a step-like manner from zero towards a positive constant at a critical value of $b$ . We begin with linear drag with coefficient $k$ , which is a suitable approximation for fish with low Reynolds number ( $\lesssim 1000$ ) such as Trinidadian guppies or zebrafish (a quadratic drag is often included in models of larger fish [17, 21, 23, 24, 25] and we will consider this in Section 3.5). The fish therefore alternates between gliding (i.e. decelerating) when $\frac{dv}{dt}=g(b)-kv<0$ , i.e., when $g(b)<kv$ , and bursting (i.e. accelerating) when $\frac{dv}{dt}=g(b)-kv>0$ , i.e., when $g(b)>kv$ . The full model is now:

	$\displaystyle\frac{db}{dt}=b-\frac{b^{3}}{3}+c-v$		(4)
	$\displaystyle\frac{dv}{dt}=g(b)-kv,$		(5)

where for $g(b)$ we choose the sigmoidal inverse tangent function:

g(b)=a\left(\frac{\arctan(z_{1}(b-b_{0}))}{\pi}+\frac{1}{2}\right).

(6)

This model again resembles the FitzHugh-Nagumo model, but with a sigmoidal rather than linear function for $g(b)$ . In fact, although the functional forms are relatively simple, the shapes of these functions (cubic and sigmoidal, respectively, in the fast variable, and linear in the slow variable) are very common in fast-slow type neural models [61].

The equation for the internal state $b$ has one main parameter, $c$ , which controls the bursting rate and amplitude. This parameter represents the tonic control of frequency and amplitude of motion that is known to come from the brain stem [62, 33].

In Figure 1, we see the time series for the one fish model for two different values of $c$ . Both $c=0.95$ and $c=1$ result in burst and glide swimming, but $c=1$ gives a higher frequency and amplitude of the speed. In the modelled scenario of a fish swimming down a channel, the fish accelerates about 4 cm down and then glide about 4.3 cm (for $c=0.95$ ). For $c=1$ the fish accelerates about 2.3 cm, and then glides about 4.8 cm. Figure 2 (a) shows phase plane trajectories for five different values of $c$ . For $c=0.9$ , there is a stable equilibrium at $(b^{*},v^{*})\approx(-1,2.3)$ , meaning that the fish do not move with intermittent burst and glide locomotion, but with a constant speed of $\approx 2.3$ cm/s. As $c$ is increased to $c=0.95$ , there is a limit cycle, meaning there is burst and glide behaviour. As $c$ increases further, the limit cycle increases leading to increasing amplitude of the speed. In Figure 2 (b) we see period length of the limit cycle as a function of $c$ . We see that there is a Hopf bifurcation at $c\approx 0.927$ , below which there is a stable equilibrium. On the other side of the bifurcation, the limit cycle period $T$ rapidly decreases from infinity before approaching a value on the order of one second for $c>1$ . This means that the higher value of $c$ , the more often the fish bursts.

3 Two fish model

3.1 Coupling between fish

With our single fish model as a starting point, we now build a two fish model. As for the single fish model, the two fish are assumed to move along a one dimensional line, e.g., swimming up a channel, and are at the points $x_{1}(t)$ and $x_{2}(t)$ respectively at time $t$ . We let $v_{1}(t)$ and $v_{2}(t)$ be the velocity of the two fish, at time $t$ , and $b_{1}(t)$ and $b_{2}(t)$ be the internal state of the fish as described in Section 2.

We couple the equations for the two fish using a response function $f(x)$ , where $x$ is the distance (negative when the other fish is behind, positive when in front) between the fish. As we saw in the single fish model, the burst frequency and the maximum speed increases as $c$ increases. As we want the fish to burst when it is behind another fish, the fish are coupled in the way that their internal state is affected by the distance between the fish, making $f$ act as an extra boost on $c$ :

	$\displaystyle\frac{db_{1}}{dt}=b_{1}-\frac{b_{1}^{3}}{3}+c+f(x_{2}-x_{1})-v_{1}$		(7)
	$\displaystyle\frac{db_{2}}{dt}=b_{2}-\frac{b_{2}^{3}}{3}+c+f(x_{1}-x_{2})-v_{2}$		(8)
	$\displaystyle\frac{dv_{1}}{dt}=g(b_{1})-kv_{1}$		(9)
	$\displaystyle\frac{dv_{2}}{dt}=g(b_{2})-kv_{2}$		(10)
	$\displaystyle\frac{dx_{1}}{dt}=v_{1}$		(11)
	$\displaystyle\frac{dx_{2}}{dt}=v_{2}.$		(12)

We want the follower fish to try to catch up with the leader. Thus, for positive $x$ , $f(x)$ is relatively large. Since most fish have wide peripheral vision, the fish will also react if the other fish is just behind, meaning that $f(x)$ is also positive (but smaller) for $x$ less than zero. It approaches zero for negative values of $x$ . To model this, we use the following sigmoid function:

f(x_{2}-x_{1})=d\frac{1}{1+\exp(-z_{2}(x_{2}-x_{1}-x_{0}))}.

(13)

This is a logistic function which has $f(x)=\frac{d}{2}$ at the threshold $x_{0}$ and approaches $d$ as $x$ becomes large, and 0 as $x$ goes to $-\infty$ . The parameter $z_{2}$ determines how steep the threshold response is: a high value of $z_{2}$ gives a steeper, more sudden, threshold response, whereas a low value gives a shallower threshold. Since most fish have a wide peripheral vision in our model, we want a relatively less steep threshold response – a fast switch would mean that $f(x)\approx 0$ for $x<0$ , and thus entail no peripheral vision. In our simulations, we set $x_{0}=15$ cm and $z_{2}=0.15$ cm^-1. This means the fish have responsiveness up to around 5 cm behind, but will react more strongly when the fish is swimming in front. From the single fish model we know that the burst and glide dynamics are highly sensitive to small fluctuations of the value of c. Thus $d$ takes values of size $\sim 0.01$ . Since $f(x)\to d$ as $x\to\infty$ , the sigmoid function assumes that the fish would interact with the other fish infinitely far in front. This is, of course, an unrealistic assumption over large distances – a fish would not interact with a fish 50 meters in front – but since we are modelling short range interactions, this assumption can be used as a simplification for now.

3.2 Temporal dynamics

We now examine the temporal dynamics of the two fish model. The model gives rise to five different behaviours: (1) constant swimming speed with no distance between the fish, (2) aperiodic leadership switching, (3) periodic leadership switching, (4) one leader-one follower swimming, (5) no leader swimming. These are shown in Figure 3, where we see time series of the difference in distance ( $x_{1}-x_{2}$ ) and the speed of the fish ( $v_{1}$ and $v_{2}$ ) for five different values of $d$ . In Figure 3 (a) and (b), where $d=0.005$ , we see that the fish undergo one burst each, and then reach an equilibrium speed of around 2.8 cm/s. At this equilibrium the distance between the fish quickly also reaches an equilibrium at $x_{1}-x_{2}=0$ . This corresponds to a steady, non-bursting swimming on the part of the fish.

When $d$ is increased to $0.01$ (Figure 3 (c)-(d)) both fish exhibit burst and glide motion. The fish undergo aperiodic leadership switching: the leadership of the fish alternates without any regular periodicity. The speed of the fish (Figure 3 (d)) is displayed from $t=130$ to $t=160$ . In the beginning of this time series, fish 1 is the leader. Fish 1 bursts about a second before fish 2, and fish 2 has a higher maximum speed in order to try to catch up. At $t\approx 143$ s there is switch in leadership and fish 2 bursts 1-2 seconds before fish 1, while fish 1 has a higher maximum speed in order to try to catch up. After 3 bursts each, there is yet another leadership switching, and fish 1 becomes the leader again.

Increasing $d$ even further to $0.02$ (Figure 3 (e)-(f)), the leadership switching happens periodically: after around 7 bursts, the fish switch leader-follower roles. There is a leadership switching at $t\approx 144$ s, where fish 1 becomes leader instead of fish 2. In the time series of the speed of the fish (Figure 3 (f)), we see that before the switching, fish 1 nearly catches up with fish 2, bursting about a second after, with higher maximum speed than fish 2. The burst by fish 1 happens closer and closer to the bursts of fish 2 on each cycle of burst and glide, and at $t\approx 143$ s, fish 1 and fish 2 burst almost simultaneously, with the same maximum speed. After this, leadership switches so that fish 2 behind and catching up to fish 1. Now fish 2 has a higher maximum speed and burst about one second after fish 1.

For $d=0.05$ , there is always one leader: whenever the following fish starts to catch up with the leader fish, the leader burst again, which we see in Figure 3 (g). In Figure 3 (h), we see that the follower fish (fish 1) has a higher maximum speed than fish 2, in order to try to catch up. Even though the follower fish has a higher maximum speed, the distance swum after one burst and glide period is approximately the same for both the follower and the leader fish.

For $d=0.07$ , there is no leader, instead the fish burst and glide with the exact same maximum speed, but take turns in who is bursting and who is gliding (Figure 3 (i) and (j)). This results in a distance between the fish that oscillates between $\approx-16$ cm and $\approx 16$ cm, with the mean distance between the fish equal to 0.

Another way to look at these interactions is to consider the speed-distance phase plane. In Figure 4 (a) we see the response function $f(x_{2}-x_{1})$ as a function of $x_{2}-x_{1}$ and in Figure 4 (b) the speed-distance phase plane is plotted. We notice that while the response function is asymmetrical, a maximum speed is reached both when the fish is in front of and behind the other fish. However, the maximum speed is higher when the fish is behind than in front, which is also what we see in Figure 3 (f).

3.3 Qualitative comparison to data

The aim of this paper is to investigate whether the model we propose qualitatively reproduces the type of motion we see between real fish. The behaviour in all of the phases, are found in experiments of pairs of fish. For example, the intermittent leadership switching in our model corresponds to stickleback swapping leadership roles [29]; the constant leader-follower dynamics for higher values of the coupling parameter corresponds to the leader-follower relationship found in, for example, pairs of eastern mosquitofish [28]; and the tit-for-tat like swimming is similar to the behaviour demonstrated by Milinksi [30].

To make a further qualitative comparison between the model and pairs of fish swimming together, we compared properties of the motion of pairs of guppies (Poecilia reticulata) to those same properties in the model. Empirical data originated from a subset of the data used in [16] where pairs of fish were allowed to explore an open rectangular arena (1000 x 900 mm) with a depth of 45 mm for $\sim$ 16 minutes. Trials were filmed with a camera recording at 1920 x 1080 pixels at 24 frames per second. Fish in the videos were tracked using the software CTrax and tracking errors corrected with the associated Fixerrors GUI in MATLAB, providing coordinates of the centre of mass and heading (in radians) of each fish on each frame. Trajectories were analysed when fish were more than 100 mm of the corners of the arena and when the difference in heading between the pair was less than 45 degrees. Full details of the empirical methods can be found in [16].

In Figure 5 we see a phase plane plot based on data from [16], compared to the model, together with a short time series from both the model and data from the same data set. We see that the model produces a symmetrical phase plane just as the data from the experiments. However, the minimum speed for the model is higher and the maximum speed is lower than the speed in experiments. In the short time series comparison, we see that the data from the experiments are noisier, but have approximately the same amount of burst and glides within 15 seconds.

3.4 Role of parameters in dynamics

We now further investigate the dynamics for the different types of swimming behaviour. Since the first swimming behaviour, constant swimming speed with no distance between the fish, consists of an equilibrium point for both speed of the fish and distance between the fish, we will not investigate this swimming behaviour in our further analyses. For the other dynamics, we will use a similar technique as first described by Lorenz [63]. When studying what later became known as the Lorenz equations, he plotted the maximum value of his $z$ variable on consecutive journeys round the continuous time attractor. In doing so, he found that the maximum obeyed something similar to a tent map [63], which is known to produce chaotic dynamics.

Inspired by this approach, we find the peaks (maximum) of the difference in distance between the fish, $x$ , and then plot the value of the $n$ th maximum, $M_{n}$ , against the next maximum $M_{n+1}$ . These plots are shown in Figure 6. For our system, this approach does not give a tent map for any of the parameter values, but it does give a clearer impression of the dynamics. In Figure 6 (a) we see that the Lorenz map for $d=0.01$ has a complicated structure, where as in Figure 6 (b) when $d=0.02$ , the Lorenz map shows a periodic pattern. In Figure 6 (c) and (d), i.e. when $d=0.05$ and $d=0.07$ respectively, the maximum value have reached an equilibrium point at 4.53 and 16.1 cm, respectively.

Another approach for investigating the differences between the various behaviours is to look at the autocorrelation function, ACF, for the peaks of the distance difference, $M_{n}$ . In Figure 7 we plot the autocorrelation function of the peaks of distance, $M_{n}$ , for $d=0.01$ and $d=0.02$ . When $d=0.01$ (Figure 7 (a)), there is a weak periodic pattern of length 2-3 in the ACF. This is most likely due to the fact that many leadership switchings happens after 2-3 bursts, but apart from that, there is no clear periodic structure. In Figure 7 (b), where $d=0.02$ , we see a clearly periodic pattern in the ACF. The periodicity is of length $\approx 14$ , which corresponds to leadership switching after $\approx 7$ bursts. Since the peaks of the distance difference for both $d=0.05$ and $d=0.07$ consists of one equilibrium point, we do not analyse the ACF for these values of $d$ .

To get further understanding of how the model is affected both by the coupling parameter $d$ , and the external input parameter $c$ , we also looked at ‘heat map’ bifurcation plots with $d$ as a bifurcation parameter, for different values of $c$ (figure 8). Each subplot is produced by, for each value of $d$ , simulating the dynamical system 10 times, with random initial position for fish 1. For each iteration, find the maximum points of the difference in distance, $M_{n}$ and the minimum points of the difference in distance, $N_{n}$ . These are stored in a histogram with values ranging from -35 to 35 cm, from which we make a heat plot of this two dimensional histogram, where colour indicates how frequent a value is.

In Figure 8 (a), we see that when $c=0.92$ , the fish move with constant speed and with no separation between them until $d\approx 0.058$ , after which the system displays aperiodic leadership switching until $d\approx 0.07$ , where there is no leader or follower, and instead the fish take turns in bursting and gliding. Increasing $c$ to 0.9255 (Figure 8 (b)), the fish move with constant speed and with no separation in distance between them until $d\approx 0.0085$ , when they instead display aperiodic leadership switching. When $d$ exceeds 0.017, the system shifts to periodic leadership switching, which is retained until $d\approx 0.0325$ . For $0.0325<d<0.058$ , there is always one leader and one follower, where the leading fish is dependent on the initial conditions. After that, the there is no leader or follower, instead the fish has swimming dynamics as shown in Figure 3 (i)-(j). Increasing $c$ to 0.9265, the system goes immediately from no separation in distance to periodic leadership switching when $d\approx 0.0035$ . After that, there is one leader-one follower dynamics, where the leading fish is dependent on the initial conditions, until $d=0.53$ , after which there is no leader or follower. When $c=0.93$ , the one fish model has passed the bifurcation point. Thus, there is always burst and glide swimming, no matter the value of $d$ . For $d\in[0,0.005]$ there is periodic leadership switching and after that the system exhibits one leader-one follower dynamics until $d=0.0445$ , where there is no leader or follower.

From all these bifurcation diagrams we see that there is a quite narrow range of values for the coupling parameter that give rise to interesting leader-follower dynamics – when $c=0.9255$ all four behaviours are found within the interval $[0.008,0.065]$ of the coupling parameter $d$ (Figure 8 (b)). The range of values for the external input parameter $c$ is even narrower: when $c$ is decreased by just 0.0055 from $c=0.9255$ to $c=0.92$ , the bifurcation diagram in Figure 8 (a), shows that the the model does not display one leader-one follower dynamics for any values of $d$ , and when $c$ is increased by just $0.0045$ from $c=0.9255$ to $c=0.93$ , the bifurcation diagram in Figure 8 (d) shows that the fish do not exhibit chaotic leadership switching.

3.5 Quadratic drag

Many design choices were made when developing the model, but the main alternative assumption that requires further investigation is the effect of having velocity decrease in proportion to $v^{2}$ . This is the natural choice from a hydrodynamical perspective, in particular for large fish. We now examine the model where both terms in the equation for $dv/dt$ are squared. This provides the same solution for the nullcline ( $dv/dt=0$ ), and we have the following set of equations:

		$\displaystyle\frac{db}{dt}=b-\frac{b^{3}}{3}+c-v$		(14)
		$\displaystyle\frac{dv}{dt}=a^{2}\left(\frac{\arctan(z_{1}(b-b_{0}))}{\pi}+\frac{1}{2}\right)^{2}-(kv)^{2}.$		(15)

In order to obtain a similar range of speeds, the speed decay was doubled to $k=0.8$ . While the nullclines are now identical, the phase-space trajectories are still of a different shape to the linear case, as can be seen in Figure 9(a). There is a threshold for bursting for a single fish at $c\approx 0.79$ , after which the bursting period decreases monotonically, until at $c\approx 0.98$ where a canard explosion takes place [64].

For interacting pairs, the same sequence of behaviours as the linear drag case is observed as $c$ increases, as shown in Figure 10. However, the speed bursts are clearly sharper than for the linear case.

Another design choice we made, which is worth mentioning here, is the choice of arctan function for $g(b)$ . This appears to be a necessary choice to get the dynamics we report here. When we tried different sigmoid functions for $g(b)$ , e.g. a logistic function, the two fish model did not display the rich dynamics described in Section 3.

4 Discussion

The primary contribution of this paper is to show that burst and glide behaviour in fish can be reproduced by a burst and recover model of neurons. The fish velocity is viewed as being analogous to the recovery variable in the FitzHugh-Nagumo model for neural dynamics, and the internal state is analogous to the membrane potential variable. The model reproduces, for different parameter values, many aspects of the intermittent locomotion found in many species of fish [14, 15, 16, 17, 18]. In particular, for low coupling the fish move with a constant speed. For slightly higher values the coupling parameter, $d$ , the fish move by burst and glide swimming and switch leadership irregularly. When increasing the strength of the coupling further, the fish switch leaders at regular intervals. For even stronger coupling there is always one leader, with the leader depending on the initial distance between the fish, similar to that seen in pairs of eastern mosquitofish [28]. These results can be replicated with quadratic drag, but the burst periods are much steeper, which makes the model less biologically realistic for guppies, but might capture the interactions of larger fish species.

In three of the four distinct forms of leader-follower dynamics, the follower fish bursts around a second after the leader (with the exact timing depending on parameter values), and has a higher maximum speed than the leader. However, the distance swum within one burst-and-glide period by the follower is always (slightly) shorter or the same distance as the distance swum by the leader within one burst-and-glide period. These qualitative model predictions would be relatively easy to test and may provide support for our model. For even higher values of the coupling parameter the fish alternate positions and have the same maximum speed and distance travelled within one burst and glide period. When one fish reaches its minimum speed, the other reaches its maximum speed.

Another testable prediction arising for our model lies in the shape of the response function $f$ . Previous experimental studies suggest that the response in speed is symmetric in the position of the other fish: the fish respond with a high speed both when the other fish is in front and behind [27]. As seen in Figure 4 (b), our model reproduces the same results but without that assumption. This emphasizes the contributions of a mathematical model like ours: without a model, the natural assumption from experimental data is that the bursting response is as strong when the other fish is behind than in front [47, 65, 66]. However, we find that a highly asymmetric response function $f$ , for which the fish responds less strongly to a conspecific behind it than in front, also gives rise to the sort of dynamics observed in experiments, when coupled with the delay induced by the ‘recovery’ of the neuronal firing. This finding is reminiscent of the result of Perna et al., who showed that (in models) two very different response functions can give rise to the same behaviour [67]. In our case, the recovery process gives the impression that the fish in front is responding to the fish behind, while in fact it has simply completed a burst and glide cycle, and is ready to burst again.

Existing models of leadership switching in fish include a random component to allow the switching to happen. For example, Markov chain models are often used to match experimental data of leadership switching [68, 29, 69]. Our model demonstrates that leadership switching can occur in a purely deterministic interaction. The irregular leadership switching in our model apparently results from a form of deterministic chaos. The Lorenz map (Figure 6 (a)) shows a clear deterministic form, while the autocorrelation indicates rapid decay in the correlation of movements. Evidence for chaos in animal behaviour has previously been presented in, for example, activity of single Leptothorax allardycei ants [70] and colonies of foraging ants [71].

Chaotic dynamics is a crucial part of neural mechanisms [39, 40, 41, 42, 43, 44] and there are existing models connecting sensory systems to neuronal dynamics [53]. Our model appears to be the first in which chaos comes explicitly from a social interaction in a perception-action loop, rather than internally from the neurons. Chaos arises from the coupling between the intake of visual stimuli and the internal dynamics of the fish response. There are increasing experimental studies linking neuronal activities in the brain to locomotion in fish [55, 54, 72]. Hopefully, our model can be used as a starting point for modelling this linkage.

While leadership switching is reasonably widespread in Nature, in our model it occurs when the coupling $c$ is close to the bifurcation point of the one fish model. Indeed, there is a very narrow range of parameters where the fish movements are chaotic. If leadership switching is ubiquitous in Nature, then one explanation for this might be an evolution of the neural parameter toward the Hopf bifurcation point at which burst and glide behaviour emerges. Similar mechanisms have been suggested in mathematical models describing foraging of ants and decision making slime mould [73, 74]. The argument in this case is that being near to a bifurcation point offers flexibility in foraging. Similarly, it has been argued that biological chaos is generally favoured by selection [75]. In our case, such an argument would suggest that evolutionary pressure takes fish responses towards the chaotic switching regime, where the fish regularly switch leadership.

While burst and glide swimming is pivotal in single fish and pairs of fish, intermittent locomotion becomes less evident in schools of fish. Instead, leadership becomes less apparent and fish move more smoothly and are synchronized when schooling. A future direction of our work would thus be to understand how this transition from turn-taking to coordinated movement occurs as the number of fish in a shoal increases.

Acknowledgements

Thank you to Torbjörn Lund, Ernest (Liu Yu) and Maria-Rita D’Orsogna for a series of interesting discussions about modelling fish burst and glide and to the Mittag Leffler Institute for hosting those discussions. Thank you to James Herbert-Read for sharing the experimental data and explaining how the experiments were conducted. Thank you to Niclas Kolm for input on the biology and to Annika Boussard for showing Linnéa the guppies in the Kolm lab.

References

[1] B. L. Partridge, “The structure and function of fish schools,” Scientific american, vol. 246, no. 6, pp. 114–123, 1982.
[2] I. Aoki, “A simulation study on the schooling mechanism in fish,” Bulletin of the Japanese Society of Scientific Fisheries (Japan), 1982.
[3] C. W. Reynolds, “Flocks, herds and schools: A distributed behavioral model,” in Proceedings of the 14th annual conference on Computer graphics and interactive techniques, pp. 25–34, 1987.
[4] A. Huth and C. Wissel, “The simulation of the movement of fish schools,” Journal of theoretical biology, vol. 156, no. 3, pp. 365–385, 1992.
[5] A. Czirók, H. E. Stanley, and T. Vicsek, “Spontaneously ordered motion of self-propelled particles,” Journal of Physics A: Mathematical and General, vol. 30, no. 5, p. 1375, 1997.
[6] A. Czirók, M. Vicsek, and T. Vicsek, “Collective motion of organisms in three dimensions,” Physica A: Statistical Mechanics and its Applications, vol. 264, no. 1-2, pp. 299–304, 1999.
[7] G. Grégoire, H. Chaté, and Y. Tu, “Moving and staying together without a leader,” Physica D: Nonlinear Phenomena, vol. 181, no. 3-4, pp. 157–170, 2003.
[8] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel type of phase transition in a system of self-driven particles,” Physical review letters, vol. 75, no. 6, p. 1226, 1995.
[9] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks, “Collective memory and spatial sorting in animal groups,” Journal of theoretical biology, vol. 218, no. 1, pp. 1–11, 2002.
[10] G. Grégoire and H. Chaté, “Onset of collective and cohesive motion,” Physical review letters, vol. 92, no. 2, p. 025702, 2004.
[11] D. Strömbom, “Collective motion from local attraction,” Journal of theoretical biology, vol. 283, no. 1, pp. 145–151, 2011.
[12] M. Romenskyy, J. E. Herbert-Read, A. J. Ward, and D. J. Sumpter, “Body size affects the strength of social interactions and spatial organization of a schooling fish (pseudomugil signifer),” Royal Society open science, vol. 4, no. 4, p. 161056, 2017.
[13] K. Tunstrøm, Y. Katz, C. C. Ioannou, C. Huepe, M. J. Lutz, and I. D. Couzin, “Collective states, multistability and transitional behavior in schooling fish,” PLoS Comput Biol, vol. 9, no. 2, p. e1002915, 2013.
[14] A. V. Kalueff, M. Gebhardt, A. M. Stewart, J. M. Cachat, M. Brimmer, J. S. Chawla, C. Craddock, E. J. Kyzar, A. Roth, S. Landsman, et al., “Towards a comprehensive catalog of zebrafish behavior 1.0 and beyond,” Zebrafish, vol. 10, no. 1, pp. 70–86, 2013.
[15] G. Wu, Y. Yang, and L. Zeng, “Kinematics, hydrodynamics and energetic advantages of burst-and-coast swimming of koi carps (cyprinus carpio koi),” Journal of Experimental Biology, vol. 210, no. 12, pp. 2181–2191, 2007.
[16] J. E. Herbert-Read, E. Rosén, A. Szorkovszky, C. C. Ioannou, B. Rogell, A. Perna, I. W. Ramnarine, A. Kotrschal, N. Kolm, J. Krause, et al., “How predation shapes the social interaction rules of shoaling fish,” Proceedings of the Royal Society B: Biological Sciences, vol. 284, no. 1861, p. 20171126, 2017.
[17] J. J. Videler and D. Weihs, “Energetic advantages of burst-and-coast swimming of fish at high speeds,” Journal of Experimental Biology, vol. 97, no. 1, pp. 169–178, 1982.
[18] G. Li, I. Ashraf, B. François, D. Kolomenskiy, F. Lechenault, R. Godoy-Diana, and B. Thiria, “Burst-and-coast swimmers optimize gait by adapting unique intrinsic cycle,” Communications biology, vol. 4, no. 1, pp. 1–7, 2021.
[19] D. Weihs, “Energetic advantages of burst swimming of fish,” Journal of Theoretical Biology, vol. 48, no. 1, pp. 215–229, 1974.
[20] D. L. Kramer and R. L. McLaughlin, “The behavioral ecology of intermittent locomotion,” American Zoologist, vol. 41, no. 2, pp. 137–153, 2001.
[21] R. Blake, “Functional design and burst-and-coast swimming in fishes,” Canadian Journal of Zoology, vol. 61, no. 11, pp. 2491–2494, 1983.
[22] F. E. Fish, J. F. Fegely, and C. J. Xanthopoulos, “Burst-anc-coast swimming in schooling fish (notemigonus cr ysoleucas) with implications for energy economy,” Comparative Biochemistry and Physiology Part A: Physiology, 1991.
[23] E. G. Drucker, “The use of gait transition speed in comparative studies of fish locomotion,” American Zoologist, vol. 36, no. 6, pp. 555–566, 1996.
[24] E. Akoz and K. W. Moored, “Unsteady propulsion by an intermittent swimming gait,” arXiv preprint arXiv:1703.06185, 2017.
[25] D. Floryan, T. Van Buren, and A. J. Smits, “Forces and energetics of intermittent swimming,” Acta Mechanica Sinica, vol. 33, no. 4, pp. 725–732, 2017.
[26] P. Paoletti and L. Mahadevan, “Intermittent locomotion as an optimal control strategy,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 470, no. 2164, p. 20130535, 2014.
[27] A. Kotrschal, A. Szorkovszky, J. Herbert-Read, N. I. Bloch, M. Romenskyy, S. D. Buechel, A. F. Eslava, L. S. Alòs, H. Zeng, A. Le Foll, et al., “Rapid evolution of coordinated and collective movement in response to artificial selection,” Science advances, vol. 6, no. 49, p. eaba3148, 2020.
[28] T. Schaerf, J. Herbert-Read, and A. Ward, “A statistical method for identifying different rules of interaction between individuals in moving animal groups,” Journal of the Royal Society Interface, vol. 18, no. 176, p. 20200925, 2021.
[29] S. Nakayama, J. L. Harcourt, R. A. Johnstone, and A. Manica, “Initiative, personality and leadership in pairs of foraging fish,” PLoS One, vol. 7, no. 5, p. e36606, 2012.
[30] M. Milinski, “Tit for tat in sticklebacks and the evolution of cooperation,” Nature, vol. 325, no. 6103, pp. 433–435, 1987.
[31] S. Mishra, K. Tunstrøm, I. D. Couzin, and C. Huepe, “Collective dynamics of self-propelled particles with variable speed,” Physical Review E, vol. 86, no. 1, p. 011901, 2012.
[32] D. S. Calovi, A. Litchinko, V. Lecheval, U. Lopez, A. P. Escudero, H. Chaté, C. Sire, and G. Theraulaz, “Disentangling and modeling interactions in fish with burst-and-coast swimming reveal distinct alignment and attraction behaviors,” PLoS computational biology, vol. 14, no. 1, p. e1005933, 2018.
[33] Ö. Ekeberg, “A combined neuronal and mechanical model of fish swimming,” Biological cybernetics, vol. 69, no. 5, pp. 363–374, 1993.
[34] T. D. Wiggin, T. M. Anderson, J. Eian, J. H. Peck, and M. A. Masino, “Episodic swimming in the larval zebrafish is generated by a spatially distributed spinal network with modular functional organization,” Journal of neurophysiology, vol. 108, no. 3, pp. 925–934, 2012.
[35] R. FitzHugh, “Impulses and physiological states in theoretical models of nerve membrane,” Biophysical journal, vol. 1, no. 6, pp. 445–466, 1961.
[36] J. Nagumo, S. Arimoto, and S. Yoshizawa, “An active pulse transmission line simulating nerve axon,” Proceedings of the IRE, vol. 50, no. 10, pp. 2061–2070, 1962.
[37] A. L. Hodgkin and A. F. Huxley, “A quantitative description of membrane current and its application to conduction and excitation in nerve,” The Journal of physiology, vol. 117, no. 4, pp. 500–544, 1952.
[38] Y. Shim and P. Husbands, “The chaotic dynamics and multistability of two coupled fitzhugh–nagumo model neurons,” Adaptive Behavior, vol. 26, no. 4, pp. 165–176, 2018.
[39] H. Korn and P. Faure, “Is there chaos in the brain? ii. experimental evidence and related models,” Comptes rendus biologies, vol. 326, no. 9, pp. 787–840, 2003.
[40] M. R. Guevara, L. Glass, M. C. Mackey, and A. Shrier, “Chaos in neurobiology,” IEEE Transactions on Systems, Man, and Cybernetics, no. 5, pp. 790–798, 1983.
[41] C. A. Skarda and W. J. Freeman, “How brains make chaos in order to make sense of the world,” Behavioral and brain sciences, vol. 10, no. 2, pp. 161–173, 1987.
[42] J. Wright and D. Liley, “Dynamics of the brain at global and microscopic scales: Neural networks and the eeg,” Behavioral and Brain Sciences, vol. 19, no. 2, pp. 285–295, 1996.
[43] W. Freeman and G. V. Di Prisco, “Eeg spatial pattern differences with discriminated odors manifest chaotic and limit cycle attractors in olfactory bulb of rabbits,” in Brain theory, pp. 97–119, Springer, 1986.
[44] P. Rapp, I. Zimmerman, A. Albano, G. Deguzman, and N. Greenbaun, “Dynamics of spontaneous neural activity in the simian motor cortex: The dimension of chaotic neurons,” Physics Letters A, vol. 110, no. 6, pp. 335–338, 1985.
[45] S. Ohgi, S. Morita, K. K. Loo, and C. Mizuike, “Time series analysis of spontaneous upper-extremity movements of premature infants with brain injuries,” Physical Therapy, vol. 88, no. 9, pp. 1022–1033, 2008.
[46] A. Strandburg-Peshkin, C. R. Twomey, N. W. Bode, A. B. Kao, Y. Katz, C. C. Ioannou, S. B. Rosenthal, C. J. Torney, H. S. Wu, S. A. Levin, et al., “Visual sensory networks and effective information transfer in animal groups,” Current Biology, vol. 23, no. 17, pp. R709–R711, 2013.
[47] J. E. Herbert-Read, “Understanding how animal groups achieve coordinated movement,” Journal of Experimental Biology, vol. 219, no. 19, pp. 2971–2983, 2016.
[48] J. Larsch and H. Baier, “Biological motion as an innate perceptual mechanism driving social affiliation,” Current Biology, vol. 28, no. 22, pp. 3523–3532, 2018.
[49] K. P. O’Keeffe, H. Hong, and S. H. Strogatz, “Oscillators that sync and swarm,” Nature communications, vol. 8, no. 1, p. 1504, 2017.
[50] J. A. Marshall, R. Bogacz, A. Dornhaus, R. Planqué, T. Kovacs, and N. R. Franks, “On optimal decision-making in brains and social insect colonies,” Journal of the Royal Society Interface, vol. 6, no. 40, pp. 1065–1074, 2009.
[51] K. M. Passino, T. D. Seeley, and P. K. Visscher, “Swarm cognition in honey bees,” Behavioral Ecology and Sociobiology, vol. 62, no. 3, pp. 401–414, 2008.
[52] T. L. Ribeiro, D. R. Chialvo, and D. Plenz, “Scale-free dynamics in animal groups and brain networks,” bioRxiv, 2020.
[53] Y. Kuniyoshi and S. Sangawa, “Early motor development from partially ordered neural-body dynamics: experiments with a cortico-spinal-musculo-skeletal model,” Biological cybernetics, vol. 95, no. 6, pp. 589–605, 2006.
[54] T. W. Dunn, Y. Mu, S. Narayan, O. Randlett, E. A. Naumann, C.-T. Yang, A. F. Schier, J. Freeman, F. Engert, and M. B. Ahrens, “Brain-wide mapping of neural activity controlling zebrafish exploratory locomotion,” Elife, vol. 5, p. e12741, 2016.
[55] E. A. Naumann, J. E. Fitzgerald, T. W. Dunn, J. Rihel, H. Sompolinsky, and F. Engert, “From whole-brain data to functional circuit models: the zebrafish optomotor response,” Cell, vol. 167, no. 4, pp. 947–960, 2016.
[56] J. M. Epstein, “Why model?,” Journal of artificial societies and social simulation, vol. 11, no. 4, p. 12, 2008.
[57] M. A. Bedau, “Can unrealistic computer models illuminate theoretical biology,” in Proceedings of the 1999 genetic and evolutionary computation conference workshop program, pp. 20–23, Citeseer, 1999.
[58] L. Gyllingberg, A. Birhane, and D. J. Sumpter, “The lost art of mathematical modelling,” arXiv preprint arXiv:2301.08559, 2023.
[59] P. E. Smaldino, “Models are stupid, and we need more of them,” in Computational social psychology, pp. 311–331, Routledge, 2017.
[60] P.-A. Braillard, C. Malaterre, T. Breidenmoser, and O. Wolkenhauer, “Explanation and organizing principles in systems biology,” Explanation in Biology: An Enquiry into the Diversity of Explanatory Patterns in the Life Sciences, pp. 249–264, 2015.
[61] F. K. Skinner, N. Kopell, and E. Marder, “Mechanisms for oscillation and frequency control in reciprocally inhibitory model neural networks,” Journal of computational neuroscience, vol. 1, no. 1, pp. 69–87, 1994.
[62] S. Grillner, “Neurobiological bases of rhythmic motor acts in vertebrates,” Science, vol. 228, no. 4696, pp. 143–149, 1985.
[63] E. N. Lorenz, “Deterministic nonperiodic flow,” Journal of atmospheric sciences, vol. 20, no. 2, pp. 130–141, 1963.
[64] M. Krupa and P. Szmolyan, “Relaxation oscillation and canard explosion,” Journal of Differential Equations, vol. 174, no. 2, pp. 312–368, 2001.
[65] J. Gautrais, F. Ginelli, R. Fournier, S. Blanco, M. Soria, H. Chaté, and G. Theraulaz, “Deciphering interactions in moving animal groups,” PLOS Computational Biology, 2012.
[66] Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, and I. D. Couzin, “Inferring the structure and dynamics of interactions in schooling fish,” Proceedings of the National Academy of Sciences, vol. 108, no. 46, pp. 18720–18725, 2011.
[67] A. Perna, G. Grégoire, and R. P. Mann, “On the duality between interaction responses and mutual positions in flocking and schooling,” Movement ecology, vol. 2, no. 1, pp. 1–11, 2014.
[68] J. L. Harcourt, T. Z. Ang, G. Sweetman, R. A. Johnstone, and A. Manica, “Social feedback and the emergence of leaders and followers,” Current Biology, vol. 19, no. 3, pp. 248–252, 2009.
[69] S. Nakayama, R. A. Johnstone, and A. Manica, “Temperament and hunger interact to determine the emergence of leaders in pairs of foraging fish,” PloS one, 2012.
[70] B. J. Cole, “Is animal behaviour chaotic? evidence from the activity of ants,” Proceedings of the Royal Society of London. Series B: Biological Sciences, vol. 244, no. 1311, pp. 253–259, 1991.
[71] L. Li, H. Peng, J. Kurths, Y. Yang, and H. J. Schellnhuber, “Chaos–order transition in foraging behavior of ants,” Proceedings of the National Academy of Sciences, vol. 111, no. 23, pp. 8392–8397, 2014.
[72] X. Chen, Y. Mu, Y. Hu, A. T. Kuan, M. Nikitchenko, O. Randlett, A. B. Chen, J. P. Gavornik, H. Sompolinsky, F. Engert, et al., “Brain-wide organization of neuronal activity and convergent sensorimotor transformations in larval zebrafish,” Neuron, vol. 100, no. 4, pp. 876–890, 2018.
[73] S. C. Nicolis, J. Fernández, C. Pérez-Penichet, C. Noda, F. Tejera, O. Ramos, D. J. Sumpter, and E. Altshuler, “Foraging at the edge of chaos: Internal clock versus external forcing,” Physical review letters, vol. 110, no. 26, p. 268104, 2013.
[74] N. Zabzina, A. Dussutour, R. P. Mann, D. J. Sumpter, and S. C. Nicolis, “Symmetry restoring bifurcation in collective decision-making,” PLoS computational biology, vol. 10, no. 12, p. e1003960, 2014.
[75] R. Ferriere and G. A. Fox, “Chaos and evolution,” Trends in ecology & evolution, vol. 10, no. 12, pp. 480–485, 1995.