Ideal free flows of optimal foragers: Vertical migrations in the ocean

We consider the collective motion of animals in time-varying environments, using as a case the diel vertical migration of marine copepods. The animals are distributed in space such that each animal moves optimally, seeking regions which offer high growth rates and low mortalities, subject to costs on excessive movements as well as being in regions with high densities of conspecifics. The model applies to repeated scenarios such as diel or seasonal patterns, where the animals are aware of both current and future environmental conditions. We show that this problem can be viewed as a differential game of mean field type, and that the evolutionary stable solution, i.e., the Nash equilibrium, is characterized by partial differential equations, which govern the distributions and migration velocities of animals. These equations have similarities to equations that appear in the fluid dynamics, specifically the Euler equations for compressible inviscid fluids. If the environment is constant, the ideal free distribution emerges as an equilibrium. We illustrate the theory with a numerical example of vertical migrations of oceanic copepods, where animals are attracted to nutrient-rich surface waters while repulsed from light during daytime due to the presence of visual predators, aiming to reduce both proximity to conspecifics and swimming efforts. For this case, we show that optimal movements are diel vertical migrations in qualitative agreement with observations of copepods.


Introduction
Animals move in space in order to forage, avoid predators, or seek mates, and the resulting movements impact the dynamics of ecosystems (Nathan et al. 2008). The spatial distributions of populations, which emerge from the movements of individuals, govern the rate with which predators meet prey, and therefore how biomass is transported from lower trophic levels to higher trophic levels, as well as the growth rate and mortality of the various populations. In turn, the movements of individual animals are motivated by the spatial structure of their habitat.
A noteworthy example of animal movement is vertical migrations in the ocean, which are ubiquitous (Brierley 2014) and believed to constitute the largest movement of biomass on the planet. Even phytoplankton display vertical migration (Wirtz and Smith 2020), while diel migrations of zooplankton and forage fish give rise to the dynamics of the deep scattering layer (Cisewski et al. 2010). Top predators generally track the movements of their prey, but take into account the limits imposed by physiological constraints such as light, oxygen, and temperature (Thygesen et al. 2016). These migrations that occur across trophic levels have implications not just for the animals that perform them, their prey and predators, but also for ecosystem functioning such as the vertical transport of carbon, a mechanism known as the biological carbon pump (Longhurst and Harrison 1989). Vertical distributions also affect measurements taken from the system and their interpretation, such as trawl surveys or acoustic surveys targeting forage fish or zooplankton (Gauthier and Rose 2005).
It is generally agreed that prey descend at day into the dark deep to avoid visual predation (Bollens and Frost 1989) and that predators follow their prey (Josse et al. 1998). Optimization arguments have been used to explain both the behavior of planktonic prey and planktivorous predators (Mangel and Clark 1988). While these motions seem to have "bottom-up" dynamics, in that the prey initiate the movement and the predators follow, they must be understood in a game-theoretic setting, since they involve the simultaneous decisions of several agents which each pursue their own interests. One phenomenon which adds credibility to the game-theoretic understanding is that the migrations can drive reverse migrations by even lower trophic levels (Ohman et al. 1983): Here, as predatory copepods descend at dawn to avoid visual predation from forage fish, their smaller planktonic prey ascend to spend the day near the surface, where the absence of predatory copepods implies relative safety.
To understand and explain these diel vertical migrations using the notions of game theory, Iwasa (1982) was seminal in modeling the phenomenon as a matrix game. This framework was pursued further by several authors including Sainmont et al. (2013) and Pinti and Visser (2019), with increasing fidelity as well as complexity, both in terms of modeling and computations. Thygesen and Patterson (2018) modeled the phenomenon in continuous space and time, using differential equations. This leads to higher fidelity and in particular resolves the narrow windows of opportunity and risk that arise at dawn and dusk, and also allows a more direct comparison to data. However, the modeling framework is more technically challenging, and to make progress, Thygesen and Patterson (2018) ignored the costs and constraints of locomotion. While this may be reasonable for larger animals, such as tunas, smaller animals such as copepods or larvae are constrained by their maximum swimming speed, and the cost of locomotion can be expected to be a significant part of their energy budgets (Sprung 1984).
Therefore, this study aims to improve on the modeling framework developed by Thygesen and Patterson (2018) by including the cost of locomotion. Thus, we formulate a game, where each player is an animal which chooses a trajectory in space such as to maximize its fitness in a Darwinian sense. We assume that there are effectively infinitely many players, so that model is a differential mean field game in the sense of Lasry and Lions (2007). To achieve a simple and tractable model structure, we focus on the single-species case, where individuals of the same species play against each other. A model with some structural similarity has recently been posed for the collective motion of human pedestrians (Achdou and Lasry 2019).
We apply calculus of variations (Liberzon 2011) to reach a characterization of the Nash equilibrium in terms of partial differential equations, which turn out to have similarities to-but also notable differences from-the governing equations of fluid mechanics. Similar observations have been made in the field of differential mean field games (Lasry and Lions 2007). A seminal mean field game in spatial ecology is that which leads to the ideal free distribution (Fretwell and Lucas 1969), and we show that our model can be viewed as a dynamic extension of this paradigm. Models based on partial differential equations are common in spatial ecology (Okubo and Levin 2001), many of these models being advection-diffusion equations. Most basically, such models assume the movement strategies of animals and determine the resulting distributions. In contrast, we seek the optimal strategies as well as the resulting distributions, across all conceivable movement strategies, but assuming that animals are rational decision-makers that have perfect global information about their current and future environment. Thus, our motivation is similar to that of Cosner (2005) and Cantrell and Cosner (2018), who identified dispersal strategies that resulted in the ideal free distribution: Starting from the axiom that each individual acts in its own interest, but do so in environments shaped by the actions of other animals, we seek equations that describe the dynamics at population level.
To be more specific, we demonstrate the approach using a test case which mimics diel vertical migrations of marine copepods. The animals are attracted to feeding opportunities near the surface, but repulsed from daylight due to the presence of visual predators; they also aim to avoid movements and being too close to conspecifics. This is an archetypal situation where animals have an incentive to migrate in response to persistent fluctuations in their environment. Our ambition is to verify that diel vertical migrations are the optimal behavioral response to this scenario, and to examine the shape of these migrations as we vary key parameters in the model. Relative to the study of Thygesen and Patterson (2018), our method offers the advantage of including cost of motion. Relative to the methodology employed by Pinti and Visser (2019), which was to formulate the problem as a matrix game and solve it using the replicator equation, our method has the advantage of representing time and space and time continuously and in particular resolve dawn and dusk. Moreover, the partial differential equations in our model lend themselves to efficient methods for analysis, including numerical analysis.
Indeed, an important element of our contribution is numerical analysis of the resulting equations. The structural similarity between the dynamics of fluids and those of optimal foragers suggests that inspiration can be found in well-established methods for computational fluid dynamics. However, these methods cannot be used off the shelf, since we pursue periodic solutions which turn out to be unstable, and therefore we establish numerical methods tailored for the specific problem.
Although our study focuses on the case of diel vertical migration in the plankton, we believe that differentiable mean field games are a natural framework in ecology that has wider applicability. The models we pursue here can cover also other cases of predictable migrations, such as seasonal migrations across latitudes. Further, the framework can be modified to cover dynamics of other state variables than position, such as individual body size. Our ambition for the study is to build a case for the general applicability of the approach, by demonstrating it on a fairly specific example. In the discussion we shall return to these perspectives.

Migrations as a dynamic game between individuals of the same species
Consider a population of animals distributed in a spatial domain ⊂ d ; the dimension d of the space may be one, two or three. We are interested in how the spatial density of animals varies over time, i.e., the function N(t, x), t ∈ , x ∈ . Animals at position (t, x) move vertically with velocity V(t, x). Our current aim is to pose equations that govern N and V. We are interested in phenomena which take place on a shorter timescale than the life span of an individual, so we ignore population dynamics and require that the continuity equation Our interest is in a bounded domain with no-flux boundary conditions, which express that no animals leave or enter the domain across the boundary: where n(x) is normal to the boundary of at x. Together, the continuity Eq. (1) and the no-flux boundary Eq. (2) express that a constant number of animals redistribute themselves in space by moving according to the velocity field V. This field V arises from decisions by the animals and thus derives from fitness optimization, as detailed in the following.
An animal at time t, position x and moving with velocity v is able to harvest energy at a rate g(t, x) but spends energy on locomotion at a rate 1 2 |v| 2 . At the same time, it is subjected to a density-dependent mortality which is large enough to affect the decisions of the animals yet small enough that the resulting population dynamics can be ignored. For simplicity we let the density-dependent mortality be governed by a coefficient 1 > 0 which is constant in space and time. Our interest is the situation of diel migrations where g and 0 are periodic in time with period T = 1 day.
The animal aims to maximize its net energy gain while minimizing its mortality and therefore has an incentive to seek out regions in space where the harvest rate g(x, t) is high and the mortality (t, x) is low, but it must do so without excessive movements. According to unified foraging theory (Mangel and Clark 1986), these conflicting objectives are traded off against each other by combining them into a net gain rate which includes harvest, cost of locomotion, and mortality. Here the parameter F converts mortality to energy and is the expected harvest of energy until the animal dies, so that probability of dying dt in a small time interval dt is energetically equivalent to a loss of energy F dt . We describe F as the "fitness" of the animal, since the energy harvested could be channeled into reproduction; this term should not be understood too literally.
For complete consistency, the fitness F would depend on time and space and would be obtained as the expected future energy harvest, i.e., where (X s , V s ) is the position and velocity of an animal which is alive at time t and at position X t = x , and which behaves optimally thereafter. S s is the probability that this animal is alive at time s > t , i.e., S s = exp(− ∫ s t (u, X u ) du) . However, we once again assume that our model covers a small time span relative to the lifespan of the individual, so that temporal fluctuations in fitness can be ignored, and therefore we can consider F constant in time. Moreover, we aim to establish a game between the animals and pursue the Nash equilibrium of this game; at this equilibrium, all animals have the same fitness, so that F does not vary with space. Finally, by focusing on a single day and realizing that the environmental conditions may vary during the remaining lifetime of the animal, F cannot be computed solely from present-day conditions. In summary, we consider F constant and a free parameter.
If an animal chooses a trajectory {X t ∶ t ∈ [0, T]} with corresponding velocity {U t = dX t ∕dt = V(t, X t )} , then its change in fitness over one period [0, T] is Our fundamental assumption is that each animal behaves optimally in the sense of choosing a periodic trajectory {(X t , U t ) ∶ 0 ≤ t ≤ T} such as to maximize this integral, for a given density function N ∶ (t, x) ↦ N(t, x) , i.e., the solution (N, V) must have the property that if all animals are distributed in space according to N(t, x), then no single animal can benefit by deviating from the strategy U t = V(t, X t ) . From a game-theoretic perspective, this is the definition of a Nash equilibrium in the mean field game, where each individual plays against the continuum of the remaining population (Lasry and Lions 2007). Mathematically, it can also be interpreted as a symmetric Nash equilibrium in mixed strategies in a two-player game, where N(t, x) represents the probability density function of the opponent. From the point of view of evolution and adaptive dynamics (Geritz et al. 1998), the model is related to the strategy U t = V(t, X t ) being evolutionary singular: If the resident follows this strategy, then no rare mutant can obtain a higher invasion fitness. Note, however, that our model does not include population dynamics and therefore we do not follow the entire program of adaptive dynamics, where the hypothetical invader meets a resident population in ecological equilibrium.
From this fundamental assumption, we can derive the variational principle that the change in fitness (6) over one period must be insensitive to perturbations X t of the trajectory, i.e., the first variation vanishes. As is standard in calculus of variations (Liberzon 2011), we reach the Euler-Lagrange equations, which are d coupled ordinary differential equations governing X t and U t : Inserting the specific forms for r, we find Here, we omit the arguments (t, X t ) for notational clarity.

If an animal has velocity
where ∇V is a matrix field with elements V i ∕ x j . We obtain: Importantly, this equation must only hold where animals are found, i.e., at points (t, x) such that N(t, x) > 0 . In places (t, x) where N(t, x) = 0 , we do not define the velocity. Any trajectory that spends time in a region with N = 0 must lead to an increase in fitness over a period which is no greater than what is obtained by the animals that follow optimal trajectories. In this study we disregard this side condition, as our numerical case has N(t, x) > 0 everywhere, but we shall return to it in the discussion.
It is instructive to notice the structural similarity between this Eq. (7) and the one describing conservation of momentum in the flow of a compressible inviscid fluids, i.e., the Euler equation (Batchelor 1967, p. 164). See also the discussion in (Lasry and Lions 2007). With our notation, conservation of momentum in a fluid would imply where (N, V) would be fluid density and velocity, u would be an external mass-specific potential acting on the fluid, and p would be pressure. We see that in this analogy, the habitat quality (g − 0 F)∕ corresponds to the potential u driving the motion. Notice that the resulting "force" on the individual animal is in the direction of decreased quality g − 0 F . This sign may seem counter-intuitive, but is consistent with the single-animal problem (Thygesen, in prep). Apart from the sign, it is an expected analogy that the external environment-specifically, spatial fluctuations in growth and mortality-drive the motion of animals. The term 1 F∇N in (7) indicate how animals should respond to the density of conspecifics. This corresponds to the pressure term in fluid dynamics, if we define the "pressure" as p(t, x) = −N 2 (t, x) 1 F∕ . Note the sign and that the force on the individual animal, arising from the pressure gradient, is toward increased densities, which also may seem counter-intuitive but follows the same logic as the external forces. In summary, optimal foragers satisfy equations that are reminiscent of those governing fluid flows, although the signs and forces and the functional form of pressure could probably not have been guessed. The underlying reason for this similarity is the principle of least action (Goldstein et al. 2002). This principle, which permeates physics, states that molecules in a fluid follow trajectories which render the socalled action integral stationary, in the same way optimal foragers follow trajectories which render the fitness integral stationary.

Steady state and the ideal free distribution
Our principal interest is on the dynamics, but it is useful to also consider steady states. This, of course, assumes that the harvest rate g and the mortality do not vary with time.
In one spatial dimension ( d = 1 ), we seek solutions where N∕ t = 0 and V ≡ 0 , i.e., each animal has chosen a fixed spatial location. Then, (7) simplifies to By integration, we find g − F 0 − 1 FN = c where c is an integration constant. This is the celebrated ideal free distribution (Fretwell and Lucas 1969): In that terminology, is the suitability of the habitat at x, which in our case decreases linearly with the density N. In the ideal free distribution, the animals distribute themselves according to N(x) such that no animal can benefit from relocating; a consequence is that all animals experience the same suitability c. The side condition then states that all void habitats ( N(x) = 0 ) have a smaller basic suitability g(x) − F 0 (x) , implying that no animal can benefit from moving into a void habitat.
We can therefore view our model as an extension of the ideal free distribution to dynamic situations and taking into account cost of locomotion. This explains that we characterize the general solution as an ideal free flow of optimal foragers. In a true steady state, consistency requires that the fitness F can be found from the integral (5) and thus satisfies for all x. Thus the integration constant vanishes ( c = 0 ), and we reach or, isolating the density N(x): If the total abundance, i.e., the integral ∫ N(x) dx , is known, then this allows us to determine F and subsequently the distribution N(x). Alternatively, if we assume steadystate also for population dynamics, then the fitness F can in principle be determined as the energetic cost to produce an offspring, and this allows us to determine the distribution N(x) and the total abundance ∫ N(x) dx.
This establishes an equilibrium solution to the governing Eqs.
(1) and (7) for the case where the parameters g and are constant in time. This equilibrium is typically hyperbolic, i.e., there exist stable and unstable manifolds consisting of solutions that converge to, or diverge from, the equilibrium. These stable solutions come into play when the environment is stationary, but the population is initially out of steady state and therefore has to redistribute to reach the ideal free distribution. Similarly, a terminal reward can be added to the model so that the population will eventually depart from the ideal free distribution along the unstable manifold to pursue the terminal reward. This hyperbolic structure is a general feature of dynamic optimization problems (Liberzon 2011). An important consequence of this, for the case of time-varying parameters, is that one should not attempt to solve the governing equations as initial value problems: Such solutions will quickly diverge along the unstable manifolds, which renders the solutions useless.

A numerical example
In this section we discuss a particular case of diel vertical migration in the ocean. The model is chosen to illustrate the mathematical framework as simply as possible, rather than to mimic a specific system. We envision a species of zooplankton, e.g., copepods, that are small enough that their movements are constrained by viscous drag and hence the power required is quadratic in the speed, in agreement with (7). We consider a relatively shallow habitat such as a continental shelf sea (alternatively, a deep lake), so that it is plausible that the simplifying condition N(t, x) > 0 holds everywhere. The copepods are subject to predation by planktivorous fish, which rely on visual detection, so that the instantaneous mortality depends on local light levels.
The copepods themselves, in turn, rely on mechanosensing, so that their feeding rates do not depend on light levels but only on the abundance of prey, which varies with depth but not with time.

Model specification
We now detail the model, i.e., the functional forms of energy harvest rates and mortalities. Specific parameters in the model are given in Table 1 and argued for at the end of this section.

Mortality
The density-independent mortality 0 stems from visual predators. We take the resulting mortality to derive directly from the local light intensity, which governs the distance at which a predator can detect its prey. Our starting point is therefore the surface irradiance, which varies with the time of day according to the periodic function The resulting irradiance is seen in Fig. 1 (left panel). Note that Thygesen and Patterson (2018) take into account astronomy and atmosphere optics; such an effort would hardly be worthwhile in the present more idealized context. Below the water surface, light decays exponentially with absorption (8) I s (t) = I 0 1 + exp (A cos(2 t∕T)) . We next introduce the local detection radius r(t, x), which is defined as the distance at which a predator can detect a prey at depth x and time t. We follow the reasoning of Aksnes and Utne (1997) and take as detection criterion that enough scattered photons from a prey arrive on the predator's retina and trigger a neural response in the predator. Aksnes and Utne (1997) take into account the spherical spread of light between prey and predator, the exponential decay of light along this line due to absorption, and saturation in the neural response, to arrive at the transcendental equation Here I is the light level from (9) and K is a saturation parameter reflecting the transformation of the light energy to neural activity. is a compound parameter which combines the optical properties of the prey with the sensitivity of the predator, and which can be characterized as the detection distance in clear water ( k 1 = 0 ) when light is abundant ( I → ∞ ). Once the detection distance r has been found, we posit that the predation mortality scales quadratically with the detection distance: Here, m is a free parameter which combines the speed, abundance and attacking success of predators. The justification is that predators cruise through the water with a speed that exceeds those of the zooplankton, effectively clearing cylindrical volumes (Kiørboe 2008a) in which the radius is r. This implicitly assumes that the density of predators is constant; a simplification that could be bypassed by explicitly modeling the distribution and movements of predators. The resulting mortality is depicted in Fig. 1, right panel.
The density-dependent mortality, in turn, is considered constant and fixed (Table 1). This is a coarse simplification (10) r 2 exp(k 1 r) = 2 I(t, x) K + I(t, x) .
of how densities of conspecifics affect the vital rates of an organism, and its justification is solely mathematical simplicity. We return to this in the discussion.

Energy harvest rate
We assume that the energy harvest rate is greatest near the surface and decreases with depth. Since we do not wish to include an explicit model of the distribution of the phytoplankton or microzooplankton that are the food source of our focal organism, we simply assume that harvest rate g(t, x) depends logistically on depth: Note that we take the harvest rate to be constant in time. This is consistent with copepods detecting prey using mechanical cues and assumes that behavior and prey density is constant during the day, which appears to be a reasonable approximation.

Model parametrization
The numerical values for parameters, as given in Table 1, are fixed mainly for illustrative purposes, but mimicking a case of an oceanic copepod of length 1 mm.
The maximum light level I 0 = 1 Wm −2 mm −1 just below the surface corresponds to blue light on a typical clear day (Wozniak and Dera 2007, p. 4), and the light amplitude A = 13 implies that the light levels at midnight are reduced with a factor exp(13) ≈ 4.4 ⋅ 10 5 corresponding to a full moon (Bond 1861). The day is roughly as long as the night, meaning that the scene is low latitudes and/ or around equinox. The absorption coefficient k 1 = 0.02 per meter corresponds to blue light in fairly clear ocean water (Wozniak and Dera 2007, p. 6). A saturation constant of K = 1 Wm −2 mm −1 in the vision of the predatory fish implies that they are generally not limited by light saturation, except partly at noon and at the surface. For the upper limit on detection distance we take = 2 cm in (12) g(t, x) = g 0 1 + exp(k 2 (x − x clin )) . The density-independent mortality rate 0 (t, x) (in per hour) as a function of depth and time of day coarse agreement with (Munk and Kiørboe 1985). With these parameters, it is the spherical spread of light between prey and predator, rather than absorption, that determines the detection probability, as detection distances are short compared to the length scale of absorption. The drag coefficient of 2 ⋅ 10 −6 Jh∕m 2 corresponds to a spherical body of diameter 1 mm and a viscosity of 1 mPa s (Sharqawy et al. 2010, Fig. 8). Here we have used Stokes' law F = 6 r v for the viscous drag, the well-known relationship between force, velocity, and power, a muscular efficiency of 26 % and a hydrodynamic efficiency of 1 % (Kiørboe 2008b), and we have converted time from seconds to hours.
For the maximum rate of energy uptake, Kiørboe et al. (2014) reports a daily egg production under ideal conditions of 65 % of the body mass, with a 36 % efficiency from intake to eggs, and a body dry weight of 13 g. We assume that the dry weight of the eggs is all protein and therefore has energy density 16 kJ/g. This gives an energy uptake of 15 mJ/hour. To reflect that uptake in the field is likely lower than under ideal laboratory conditions, we round this down and use a maximum uptake of g 0 = 10 mJ/hour. For the vertical profile, we take a fairly gradual transition with a cline at x clin = 100 m and a 1∕k 2 = 50 m transition zone, noting that the vertical structure of the oceans display considerable variation (de Boyer et al. 2004).
The mortality parameter m is fixed so that the maximum predation mortality is m 2 = 0.01 pr. hour, and the fitness F corresponds to the energy harvested over a time span of 200 hours. These mortalities and life spans are within the range reported by Hirst and Kiørboe (2002). Note that we study a single day during the yearly cycle, and that constant uptakes and mortalities can be added to the model without changing the optimal migrations, so that fitness does not need to balance growth and mortality exactly.
Realizing that many of these parameters are very coarsely estimated, we use these parameters as a baseline and investigate also alternative scenarios, where selected key parameters are modified. We return to the issue of parameters in the discussion.

Computational scheme
We determine the periodic-in-time solution to the Eqs. (1) and (7) with no-flux boundary conditions (2) numerically by discretizing time and space and solving the resulting system of algebraic equations using a Newton method. To ensure convergence of the Newton method, we apply a homotopy perturbation method (Alexander and Yorke 1978;He 1999). Specifically, we include a homotopy perturbation parameter which modifies the original Eqs.
(1) and (7) to the following system Following the homotopy perturbation principle, the case = 1 recovers the original system of equations, while for = 0 we know an analytical solution: N(t, x) is constant in space and time, and V(t, x) ≡ 0 . Increasing the parameter from 0 to 1, we track the solution, until the solution to the original system is obtained with = 1.
For a given value of ∈ (0, 1] , we find the solution using an iterative Newton method. At each iteration, we linearize the equations around the previous guess and solve the resulting linearized PDE system. Our experience with the numerics is that these systems are all well posed so that the good initial guess from the previous value of ensures that the Newton method converges swiftly. We have not pursued the theoretical question if these systems of equations are well posed. We discretize both time and space using the spectral collocation method. As functional basis we choose the Fourier basis with equidistant quadrature nodes for time, since we pursue periodic-in-time solutions, while for space we use a Legendre basis with Gauss-Lobatto quadrature nodes, respectively (Kopriva 2009). The computations are performed in Python. Figure 2 displays the density N(t, x) corresponding to the baseline parameters in Table 1 (top left panel). During the daytime, the animals move to deeper waters to avoid visual predation, while at dusk they return to upper levels to benefit from better feeding conditions during the night. The panel also displays the trajectories of 3 individual animals, chosen at the 25, 50, and 75 percentile of the distribution, respectively. The figure also displays the results of modifying the underlying parameters. In the top right panel, the cost of motion has been increased by multiplying the parameter with a factor 10. As expected, the effect of this is that the animals perform migrations with narrower range. In the bottom left panel, the parameter 1 describing densitydependent mortality has been increased with a factor 3, resulting in a less concentrated distribution of animals. In the bottom right panel, the fitness F has been increased with a factor 1.5, corresponding to a situation where conditions are expected to improve, or where an extra constant source of food is available. The effect of this is that the trade-off between energy harvest and mortality is altered, so that the animals pursue more cautious strategies, i.e., follow deeper trajectories. Figure 3 shows the environment that a single animal experiences during the course of a day. The quantile in these plots refer to the vertical distribution of animals, so that an animal with quantile 0 is always the closest animal to the surface while an animal with quantile 0.5 always has 50% of the animal population above it. Note that the model assumes that animals maintain their quantiles during the day, i.e., they do not pass each other, as this would be suboptimal. In the left panel, we see the light experienced by the animals. We see that animals, despite migration downward at dawn, experience much higher light levels during the day, although the light levels gradually decreases around noon, as the animals continue their migration. In the right panel, we see the instantaneous rate of change of fitness, i.e., the integrand r in (6) given by (4). We see that animals closer to the surface experience greater daily fluctuations, where the nighttime provide positive change due to feeding opportunities under relative safety, while the daytime poses threats of mortality that outweigh feeding opportunities. For the deeper animals, the pattern is qualitatively identical but less pronounced, as these animals experience less fluctuating environments. Figure 4 shows the contributions to the change in fitness, integrated over the single day, for different animals given in terms of their quantiles. The left panel shows the net growth (uptake minus cost of motion) as well as loss of fitness due to mortality, while the right panel focuses on cost of transport. Note the different scales on the two horizontal axes; the cost of motion is orders of magnitude lower than  Table 1 Fig . 4 The components in the daily change in fitness. In both panels, the horizontal axis gives the change in fitness (in joules) while vertical axis gives the quantile of the animal in question, as in Fig. 3. Left panel: Net growth (uptake minus cost of motion) (blue), loss of fitness due to mortality (red), and total gain of fitness (black) . Right panel: Total cost of motion (green) the other fitness components. In the left panel, we see that all animals have the same loss of fitness of approximately 0.07 J during the day, which follows from the construction of a Nash equilibrium between the animals. However, the different animals reach this total in very different ways: The deepest animals have hardly any growth and low mortality, corresponding to an expected life time of approximately 25 days, while the animals closest to the surface only have an expected lifetime of about 7 days, but compensate with significant growth. In the right panel, we see that the cost of motion is generally low in comparison to the other components, but also that the animals in the middle spend the most energy on migrating.

Discussion
We have considered a population distributed in space, where each individual moves to maximize its fitness, and shown that the governing equations for the bulk movement of the population coincides with those that describe an inviscid compressible fluid.
Our case study concerns diel vertical migrations of a single species in the aquatic regime, but the problem and the approach has wider applicability in ecology. Most directly, it could be adapted to seasonal migrations, taking into account both costs of migrations and density dependence. In the aquatic regime, seasonal migrations driven by fluctuations in food availability and predation risks are common among the fishes and have been investigated with optimization models, for example roaches migrating between lakes and streams (Brönmark et al. 2008). We speculate that the approach could be applied to seasonal migrations of birds, insects, and terrestrial herbivores. Each of these applications possesses features that would require adaptation of the framework, such as intermittent migration of herbivores (Owen-Smith et al. 2010) and the effects of wind on longrange migration of insects (Chapman et al. 2015). Since the analysis provides both spatial distributions and trajectories of individual animals, data both at individual and population level can be used to inform, parametrize, and validate the models. The outcome of such modeling studies could be to predict changes in distributions and migratory behavior as results of changing environmental conditions. At a higher level of abstraction, the approach can be generalized from movements in physical space to describe also changes in other state variables of the animal, such as body size. This would make the approach applicable to sizestructured communities (Andersen 2019), where organisms face a growth/mortality trade-off which is influenced by the size composition in their community. In a broader sense, our study was motivated by the quest for mathematical techniques to support the trait-based approach to ecology (Kiørboe et al. 2018). Here, we believe that a natural and promising starting point is mean-field games, where individuals aim to maximize fitness in environments that are shaped by the behavior of other animals.
Many of these potential applications involve multiple interacting species, so a topic of future studies is to include predator-prey interactions or competition between multiple species. While this is straightforward, in principle, the modeling efforts will grow with the complexity of the system, as will the computational challenges, so case studies must be designed carefully. A first step could be to apply the present modeling framework to the system of Thygesen and Patterson (2018).
When parameters are constant in time, the ideal free distribution emerges as a special solution, and thus our model is a dynamic generalization of the ideal free distribution. There are several alternative approaches to such dynamic ideal free distributions: Cosner (2005) posed a partial differential equation, the equilibrium solution of which is the ideal free distribution. Cantrell et al. (2021) used the framework of adaptive dynamics to derive an evolutionarily stable dispersal strategy, and recently, Cantrell et al. (2010) extended the framework to path-dependent fitness in time-periodic environments. A key distinguishing element between models is the assumed information available to the animals. In standard chemotaxis models (Okubo and Levin 2001), animals infer on local conditions from sensory inputs, and Cantrell and Cosner (2018) showed that local information could yield the ideal free distribution in a stationary environment, but not in a time-periodic one. Here, we have made the idealization that animals are perfectly aware of both current and future conditions in the entire domain. Thus, our model can explain motion in response to periodic fluctuations in conditions, which have remained stable over evolutionary timescales. Diel and seasonal migrations are the most immediate examples of such motion.
We have not addressed if and how actual animals can follow the optimal trajectories that we have identified. In our case study, the solutions resemble qualitatively the widely documented diel vertical migrations of zooplankton (Klevjer et al. 2016), but we have refrained from a more quantitative comparison at a specific time, place, and for a given species. It is plausible that similar motions can be implemented with simple behavioral rules, for example, pursuing a constant light intensity while avoiding high densities of conspecifics as and excessive movements. The rationale behind the current study is that evolution has established some mechanism which yields patterns that resemble optimal ones. An interesting avenue of future research would be to extend the sensitivity study lying under Fig. 2 to examine also to which degrees such different trajectories can be obtained by simple behavioral rules. Such a study would arguably be more worthwhile in a more specific setting, parameterized for a particular species and system, as opposed to the idealized and generic settings of the current study.
In the periodic solutions we establish, all animals obtain by definition the same increase in fitness over the day, even if this increase may be distributed over the period differently for different animals. However, the animals closer to the surface obtain this increase by aggressively pursuing high growth rates, thus also enduring high mortalities. In contrast, the deeper animals have slower growth and lower mortality; also, different animals incur very different costs of locomotion. Thus, the population spans a niche in which different zones call for different specialized adaptations, and it is plausible that there are limits to how wide this niche can be while still remain in the scope of the species. These considerations are beyond the present modeling study, but would be interesting to pursue.
At a technical level, we established the governing equations using calculus of variations. The most common approach in studies of mean field games is to include diffusion and pose a Hamilton-Jacobi-Bellman equation, which would govern the optimization problem of the individual animal, in which the density of conspecifics enter as data.
Here, without diffusive motion in the model, calculus of variations seems more direct and also yields the similarity with the Euler equations of fluid dynamics. However, the model can yield regions that are completely void of animals due to disadvantageous conditions at a given point in time, in the same way the ideal free distribution may predict void regions. Voids are numerically challenging, so we chose parameters for our numerical case study to avoid such regions. It would be interesting to investigate if such void regions can be resolved numerically using techniques from free surface flows, but pragmatically they can also be avoided by including a small diffusivity as a regularization parameter. With diffusion, the Hamilton-Jacobi-Bellman approach appears stronger. We intend to pursue more efficient numerical methods in future studies, so that wider regions in parameter space can be explored, and so that regions void or almost-void of animals can be investigated.
Computational practicalities aside, ecological reasoning would justify the inclusion of diffusion in the model of motion (Okubo and Levin 2001): The small aquatic animals we have in mind generally move unpredictably, which can be modeled as biased random walks and represented with diffusion terms in the governing equations, and are subject to turbulence which can also be represented by diffusion. In the interest of sparsity we have omitted diffusion from the model, noting that the main effect of diffusivity-to make animals spread out-is found in the model thanks to the density dependence of fitness, which gives the animals an incentive to disperse.
Our numerical analysis has employed spectral methods, while the common choice in numerical analysis of mean field games is finite difference methods (Achdou and Capuzzo-Dolcetta 2010; Achdou and Lasry 2019). The time periodicity of the domain, and the expected smoothness of the solution, motivated our choice. While spectral methods in this case give relatively high fidelity with relatively few basis functions, it is more difficult to establish properties of the discretized operators which guarantee convergence of the iterative solvers. Although the homotopy principle performed well, we believe that this issue underlies the numerical problems we encountered when pushing the limits of parameters space, in particular, as regions in space become almost-void of animals. We aim to investigate this further in future studies.
Our model is built using simplified functional forms, which could be expanded in future studies. First, our cost of movement is quadratic in the speed. This describes the energy spent on propulsion in creeping flow, which is justifiable for small animals. Other forms could include inviscid flows, which would lead to a cubic relationship between speed and cost, lost opportunity as described, e.g., by Thygesen et al. (2016), or increased conspicuousness to predators . A second simplification is that the cost depends linearly on the density of conspecifics. A mechanistically based density dependence is not trivial to include in the model, as it could involve both interference in the foraging of conspecifics and attraction of predators to patches of conspecifics. Both of these extensions point in the direction of two-species models, or more generally multi-species models, which we aim to address in future studies. In the present study, the role of density dependent is to avoid singular distributions, i.e., all animals following the same trajectory, whereas the specifics of the density-dependent term is rather arbitrary. Another simplification is the form of the nutrient uptake, i.e., the term g in (4). It simple to extend this uptake to reflect a specific habitat with higher fidelity, but a more substantial extension, which would require quite different techniques for analysis, would be to include the dynamics and constraints of stomach fullness.
In our case study, the parameters have been chosen with a view toward oceanic copepods, but without being explicit about the system and the species. The justification for this is that the case is primarily motivated to illustrate the theoretical framework, and it must be acknowledged that several key parameters are not well known; in particular, the density dependence and the fitness. Therefore, sensitivity studies remain an important element of any study.
Our model is based on a time scale separation, where we argue that population dynamics can be ignored when studying a single day, even if mortality is high enough that animals take it into account when choosing diel strategies. The copepods in our case have a life time of weeks, which is arguably around the lower limit for the time scale separation to be justifiable. It would be easy to include mortality in the conservation equation, and correspondingly make the fitness a function of time, but we would then pursue quasi-periodic solutions rather than strictly periodic ones, and the benefit of including mortality explicitly would be questionable. Reproduction would be more difficult to include, as such a model would probably also have to include seasonal fluctuations in conditions, the stage structure of the populations in question, and the life history of the individual. Extending the model in these directions would perhaps be simple from a conceptual perspective, but tedious in terms of modeling and computations, and we believe that such an effort would not result in increased clarity compared to our approach of investigating a fast periodic pattern during which life history can be neglected.
We have emphasized the analogy between our governing equations and the Euler equations of fluid dynamics which describe compressible inviscid fluids. In this analogy, the density-independent terms in growth and mortality serve as an exogenous potential that varies in time and space and which drives the motion, while the density-dependent terms correspond to the thermodynamic concept of pressure. This analogy gives conceptual insight which may appeal in particular to ocean ecologists, as fluid flows already play a fundamental role in shaping oceanic habitats as well as the vital rates of animals that move and forage in these habitats. The analogy also gives access to a rich toolbox for analysis, including established numerical methods.
Nevertheless, there are important differences between the motion of animals and those of fluid particles, even within the model. First, while the potential in a physical model corresponds to net growth in a behavioral optimization model, we expect a compressible fluid to concentrate where the potential is lowest, while we expect animals to aggregate where the net growth rate is highest. Similarly, it may seem disturbing that the "pressure" of animals is negative and decreases with the density. These observations are related to the solutions of interest: While focus of fluid mechanics is, generally, on (Lyapunov) stable solutions or attractors, the solutions of interest for the motion of animals are intrinsically unstable when viewed as initial value problems.
This also contrasts our model to the typical advectiondiffusion equations in spatial ecology: In our case, the solutions are stabilized by the animals looking ahead in time when taking decisions and responding to anticipated events such as sunrise, so information travels backward while densities evolve forwards. This is the situation for optimal control problems (Liberzon 2011) as well as for mean-field differential games (Lasry and Lions 2007; Achdou and Capuzzo-Dolcetta 2010). For numerical analysis, one can therefore not rely on time stepping (neither forward nor backward), which is the reason we in this paper have focused on timeperiodic solutions, as was also done by Achdou and Capuzzo-Dolcetta (2010). Alternatively, one may consider mixed initial/terminal value problems, as was done by Achdou and Lasry (2019) for human pedestrians.
In conclusion, we have established a framework consisting of partial differential equations governing the optimal motion of individual members of a population in response to periodically fluctuating environmental conditions. The equations are similar to those governing fluid flows, so that their solutions describe ideal free flows of optimal foragers. We have investigated a particular case, in an idealized setting, of planktonic copepods performing diel vertical migrations in the water column in response to sunrise and sunset, and shown that the model predicts patterns which are similar to those empirically observed. In general, we believe that the framework of differential mean field games has wide applications in ecology. More specifically, we believe that viewing migrating animals as elements in ideal free flows can illuminate animal migration and offers rich opportunities for further development, both in terms of theoretical analysis and specific case studies.
Author contributions UHT conceived the study and the analysis. MM detailed the analysis and the numerical implementation. UHT drafted the initial manuscript and both authors completed the manuscript.
Funding MM is partially funded by The Centre for Ocean Life, a VKR Centre of excellence supported by the Villum Foundation.
Data availability Not applicable.

Code availability
The code used to perform the calculations and generate the figures will be placed in a publicly accessible GitHub repository after acceptance of the manuscript.

Declarations
Ethics approval Not applicable.

Conflicts of interest
The authors have no relevant financial or nonfinancial interests to disclose.