Modelling optimal behavioural strategies in structured populations using a novel theoretical framework

Understanding complex behavioural patterns of organisms observed in nature can be facilitated using mathematical modelling. The conventional paradigm in animal behavior modelling consists of maximisation of some evolutionary fitness function. However, the definition of fitness of an organism or population is generally subjective, and using different criteria can lead us to contradictory model predictions regarding optimal behaviour. Moreover, structuring of natural populations in terms of individual size or developmental stage creates an extra challenge for theoretical modelling. Here we revisit and formalise the definition of evolutionary fitness to describe long-term selection of strategies in deterministic self-replicating systems for generic modelling settings which involve an arbitrary function space of inherited strategies. Then we show how optimal behavioural strategies can be obtained for different developmental stages in a generic von-Foerster stage-structured population model with an arbitrary mortality term. We implement our theoretical framework to explore patterns of optimal diel vertical migration (DVM) of two dominant zooplankton species in the north-eastern Black Sea. We parameterise the model using 7 years of empirical data from 2007-2014 and show that the observed DVM can be explained as the result of a trade-off between depth-dependent metabolic costs for grazers, anoxia zones, available food, and visual predation.

The complex behavioral responses and sophisticated life traits of organisms which are observed in the natural world are often considered to be outcomes of long-term evolutionary processes, and their quantitative description via mathematical modelling is usually challenging. Moreover, it is intuitively clear that optimal behaviour and/or life traits of an organism should gradually alter with maturation and progression through different developmental stages: a successful behavioral strategy for juveniles may not be effective for adults, which often experience a different environment. Structuring of populations in terms of size or developmental stage usually greatly enhances the complexity of our evolutionary models, and consequently existing modelling frameworks dealing with evolution in structured populations are somewhat less developed compared to those for unstructured populations 1,2 . The aim of this paper is to propose a framework to reveal optimal strategies for different developmental stages in stage-structured population models with continuous growth. Our study is based on a revisited concept of evolutionary fitness related to the long-term selection process and can be applied to deal with both scalar and function valued life traits or behaviours. As an important ecological study case, we explore patterns of regular diel vertical migration of zooplankton in aquatic ecosystems under variable environmental constraints.
The conventional wisdom of evolutionary modelling used to predict optimal life traits or behaviour is often based on the generic idea that a certain quantity known as the evolutionary fitness should be maximised [3][4][5][6] . However, the definition of fitness for an organism or a subpopulation is generally subjective and may depend on the personal preference of the researcher. For example, the existing models of diel vertical migration of zooplankton use different definitions of population fitness such as the individual reproductive value 7,8 , the inverse mortality 9,10 , the ratio between the food intake and the mortality [11][12][13] , and the 'venturous revenue' 14 among others. As an unfortunate result, models of optimal behaviour using different definitions of fitness can predict clearly distinct patterns [14][15][16] . In recent works, however, a new mathematically straightforward approach for identifying the evolutionary fitness has been proposed [16][17][18] . This approach considers long-term dynamics of competing subpopulations which are described by different inherited units (e.g. behavioural strategies, life traits, genotypes). The evolutionary fitness can be defined based on the comparative ranking order of such subpopulations. As a result, the long-term evolutionary outcome only depends on the choice of the underlying population model rather than a speculative metric of fitness 16 .
Here, we extend the ideas of the approach mentioned above 16,17 to explore optimal behavioral patterns in a generic age-or size-structured population model with continuous age or size settings. In particular, we demonstrate how the evolutionary fitness can be derived for a generic von Foerster model with an arbitrary mortality term. Then we show that the optimal patterns of behaviour which are observed in the model can be derived by applying the variational principle of natural selection. Using this principle we obtain equations providing the optimal trajectories of regular vertical migration of zooplankton.
Regular diel vertical migration (DVM) of marine and freshwater zooplankton is often regarded as the largest synchronised movement of biomass on Earth 19,20 . Understanding DVM is so important because it plays a key role in the carbon exchange between the deep and surface waters, the oceans biological pump 21,22 . Typically, DVM consists of planktonic grazers ascending to plankton-rich surface waters to feed at night, then descending to deeper waters and remaining there during the day 19,23 . There exist several explanations of what causes organisms to migrate, but the most accepted hypothesis is currently that zooplankton perform DVM to avoid visual predation, mostly by planktivorous fish, by staying in the deeper and darker areas during daylight hours and ascending at night when visual predators are unable to see them [23][24][25][26] . The DVM phenomenon has been extensively studied both empirically [23][24][25][26] and theoretically using a number of mathematical models 8,[26][27][28] . However, we should stress that most of the existing models of DVM were based on the principle of maximization of some initially prescribed fitness criterion which was a personal choice of researcher, thus the generality of the modelling results remains questionable. Furthermore, the role of key factors affecting DVM is still poorly understood. This particularly concerns the difference in migration behaviour for different developmental stages and the influence of environmental factors such as food, temperature, oxygen, and predators on the amplitude and timing of migration for each stage. We argue that mathematical modelling based on our revisited concept of fitness and also backed up by long-term empirical observation can provide us with a better understanding of zooplankton DVM 7,13,27,29,30 .
Here we implement the novel theoretical framework to explore the DMV patterns of dominant zooplankton herbivores in the north-eastern Black Sea ecosystem. We differentiate between strategies of different zooplankton developmental stages and consider how key environmental parameters such as food availability, predation, and habitat size affect the optimal DVM. Using the model, we address the long-standing question of why the migration depth of zooplankton is two-three times greater than the size of the euphotic zone [31][32][33] . Our modelling research into optimal DVM is backed up by 7 years of empirical observation (2007-2014) of the migration of two abundant copepod species, Calanus euxinus and Pseudocalanus elongates, across seasons and under different biotic and abiotic conditions. Our study demonstrates that it is the depth-dependent variability of the metabolic costs of the grazers -related to the specific oxygen regime in the Black Sea -rather than trophic pressure by visual predators that determines the choice of the lowest migration depth. The model also predicts that the absence of diel migration of earlier stages of zooplankton might be mostly due to their high mortality rather than due to a high cost of DVM, contrary to what was believed previously.
The paper is organized as follows. Section 2.1 introduces the generalized variational principle of modelling natural selection and introduces the generalised fitness function. In Section 2.2, a generic size-or age-structured model is introduced and the evolutionary fitness for this model is derived. In Section 2.3, the generic model is implemented to explore patterns of optimal zooplankton DVM. In particular, Subsection 2.3.1 provides empirical evidence for zooplankton DVM in the Black Sea; subsections 2.3.2 explain parameterisation of the model coefficients; Subsection 2.3.3 shows the modelling results on the optimal DVM. The discussion in Section 3 summarises the main results and provides further ideas on revealing optimal behavioral strategies using the new framework.

Results
Establishing general variational principles of natural selection. Here we provide a generic mathematical framework of modelling natural selection in a self-replicating system with inheritance and introduce a mathematically rigorous definition of evolutionary fitness.
Consider some population, where organisms are described by inherited elements v, strategies or life traits, for instance. Mathematically, an element v can be a scalar, a vector, or a function, so we consider that v belong to a certain function space V. In particular, for a stage-structured population, v can be a vector of functions describing each developmental stage. From now on we will refer to v as a strategy for the sake of simplicity. For simplicity, we consider that strategy v is passed unchanged to each offspring from its parent, i.e. as in the case of a clonal reproduction (note that our methodology can be extended to a more complicated case allowing for mutations 16 ). We assume that v belongs to a compact domain in a metric space V equipped with a Borel measure μ*. We also assume that the measure of any open set (except ∅) is greater than zero.The need for introducing a metric and using a measure comes from the following.
We assume that elements of V can be somehow compared in terms of similarity (or closeness) to each other. Mathematically, this means that we need to introduce a distance between elements v; thus the space V is required to be a metric space. The existence of a metric allows us to generate neighbourhoods and open sets in V. Our framework should also be able to quantitatively describe and compare the sizes of different sets in V. For this purpose, we use the concept of σ -algebra of ∑ subsets in V and then introduce a measure of those subsets. Note that almost for any choice of subsets ∑, there will be always non-measurable sets. Thus, it is logical to use only such σ -algebra which will contain all subsets which are of interest for further applications (e.g. all open sets). The minimal σ -algebra containing open sets is known as a Borel algebra, and a measure introduced on it is called a Borel measure. Note that in practical applications, almost all considered measures are Borel measures. For examples, in the finite dimension space, a Borel measure will be extension of the 'classical' concept of volume. In infinite dimension spaces (function spaces), a Gaussian measure is typically used 34 .
The presence of the subpopulation with strategy v at time t is indicated by a non-negative numeric value η v t ( , ) which can be understood as a generalized density. This may be the density of the subpopulation with strategy v, but we can also describe the presence of strategies in the population via the logarithmic scale of the biomass, for example. Alternatively, we can characterize the presence of strategies in the population via any positive power of population size which can be strategy-dependent. We postulate that η v t ( , ) should satisfy the following requirements: ) approaches zero this signifies extinction of v. • It is a continuous function over the space V and it is integrable over V with regard to μ*.
• It is a smooth function of time.
for any t. This takes into account the natural assumption that limitation of resources for the population restricts its population growth.
We suggest that we know the equation which govern the temporal dynamics of η v t ( , ); for example, it can be derived from the underlying model of population dynamics.
We can now introduce the following definition of ranking of strategies.
Definition 1 (Ranking order of strategies). We state that element v′ is better (or fitter) than element w′ ( ) of point w′ such that the ratio of generalised densities (1) tends to zero uniformly in the neighborhoods uniformly Note that the above definition should be considered as a partial ordering and may also depend on the initial conditions in the system. Using this definition of the ranking order of two strategies one can easily conclude that the subpopulation with a strategy of lower ranking should eventually go extinct. This is formally given by the following theorem.

Theorem 2.
If v′ is a better strategy than w′ ( ′ ′  v w ) then there exists a neighborhood ′ O w ( ) such that the generalised density η w t ( , ) tends to zero with → + ∞ t uniformly in ′ O w ( ). The proof of the above theorem is given in Supplementary Material SM1(i). Using the above definition of ranking and the theorem, one can now define the evolutionary fitness in the system as follows.

Definition 3 (Evolutionary fitness).
In the case where there exists a functional J v ( ) such that it preserves the ranking order of strategies given by, i.e. J v J w v w ( ) ( )  > ⇒ , this functional is referred to as an evolutionary fitness.
Maximizing a given fitness function provides the variational principle of modelling natural selection: only the strategies v* realizing the global maximum of the fitness J will remain in the population, the others will go extinct. Thus, the strategies v* will be evolutionary optimal. We should stress again that, since the ranking order for a particular population model may depend on the initial conditions, the evolutionary fitness can also depend on the initial conditions. Also, the fitness functional is not unique: any increasing function of J will also be an evolutionary fitness.
According to the formulated variational principle of selection, we need to find an evolutionary fitness based on model equations and then find the strategy or strategies which maximise its value. Note that in general, finding J can be challenging. In the next section we will show how the evolutionary fitness can be derived for an age-or size-structured population model described by the von Foerster equation.
Deriving evolutionary fitness in an age (stage)-structured model. We consider a single population model with structuring described by a von-Foerster-type equation 35,36 . The population is characterised by the density z v W t ( , , ) at the moment of time t with body weight W and behavioral strategy v. Here by v we understand the overall set of strategies across all developmental stages, where each strategy can be a function with known mathematical formulation with a fixed number of evolving parameters or can be an unspecified function-valued trait, essentially an infinite number of evolving parameters. The evolution is governed by: where r v W ( , ) describes the increase in the body weight of individuals due to growth; A v W ( , ) is the linear mortality rate. The second mortality term in the right hand side describes effects of competition between adults and juveniles across all possible strategies as well as possible effects of predation, harvesting, or any external forcing. Note that in this term, we incorporate the effects of strategy into R, and consider that y t ( ) equally affects the mortality of each stage. www.nature.com/scientificreports www.nature.com/scientificreports/ The production of offspring which initially have a minimal body weight W min is due to the reproduction of the whole cohort of adults which is given by is the reproduction coefficient. Biologically, W A is the weight at which organisms become mature and can reproduce, after which > b v W ( , ) 0. The increase of the body weight W is described by the following growth equation: min where r v W ( , ) is the body growth rate. Model (2) with structuring in terms of body weight can be transformed 37 to the equivalent model with structuring in terms of age τ by introducing the following variable ρ τ and the renewal equation becomes , ) should be understood as a re-scaled population density; τ ⁎ v ( ) is the time at which organisms mature and can begin to reproduce.
The weight W and the age τ are related by ( , ) 0 which should guarantee a one to one correspondence between τ and W), in particular, the maturation age τ* is given by A . For simplicity, we will further consider model (5) -(6) for discrete stages = ... i n 0, 1, , . For each stage i and for a given strategy v we assume the model coefficients to be constant. Thus, for the mortality and production rates we have Here we formally assume the reproduction to be possible from stage = i 1. However, the model still allows us to consider the case where some reproduction rates = b 0 i for > i 0. Note also that the reproduction of the final stage ends at the age τ + n 1 , whereas the organisms can live for longer. The evolutionary fitness is in the above structured model given by the following theorem. Theorem 4. In the age-structured model (5)-(6) with discrete stages = ... i n 0, 1, , , the strategy which maximises the evolutionary fitness will outcompete the other strategies. The fitness is defined by where λ i is the solution (the eigenvalue) of the characteristic equation provided in SM1(ii); R denotes the real part of this eigenvalue. The proof of the above theorem is given in supplemental Material SM1(ii).
In this paper, we will consider the case where the number of developmental stages is 3 (i.e. = i 0, 1, 2) and the reproduction starts from the oldest stage only (i.e.
, it is shown that in this case the characteristic equation determining the fitness J is given by Thus, to find the optimal strategies for all three stages = v v v v ( , , ) 1 2 3 , one needs to maximise fitness given by (7). Note that in the given model the fitness does not depend on initial conditions. Study case: exploring optimal DVM of zooplankton. Empirical observation of DVM. We now use the theoretical framework introduced above to explore patterns of diel vertical migration (DVM) of herbivorous zooplankton. The modelling study is motivated by our long-term empirical observation of DVM of zooplankton (2019) 9:15020 | https://doi.org/10.1038/s41598-019-51310-w www.nature.com/scientificreports www.nature.com/scientificreports/ in the north-eastern Black Sea. Samples were collected throughout all seasons in the years 2007-2014 for two of the most abundant herbivorous mesozooplankton species: Calanus euxinus and Pseudocalanus elongatus. Details on samples collection and data analysis are provided in the supplementary material (SM2) as well as in a satellite paper 38 . Figure 1 show a typical pattern of DVM for the investigated species; the graphs represent the variation in the mean depths of the copepod distribution across the 24 h period for all developmental stages. Note that the zooplankton population is scattered around a certain depth in each instant, and the figures show the spatially average depth of the zooplankton distribution in the column (calculated using the computational algorithm from SM2). One can see that for both species only the copepods from the older developmental stages (starting from CIV) exhibit a pronounced migration: they feed on phytoplankon at depths of h = 35-45 m at night and stay in deeper waters (h = 110-120 m) during the daytime.
Interestingly, the amplitude of DVM for both species varies across seasons, as is demonstrated in Fig. 1C,D, which show the lower and upper depths for migrating females through the daily cycle (denoted by the open and filled squares, respectively) as well as the daily average depth of non-migrating stages of the same species (denoted by semi-filled circles). One can see that the deepest depth of DVM increases in summer and decreases in winter. In the literature, seasonal variation of migration depths is often connected to seasonal change of the width of (i) the oxygen zone in the north-eastern Black Sea and (ii) the zone with suitable water temperature ranges for the species considered 32 . However, the role of the oxygen zone and temperature conditions in DVM are not yet well understood and this model study is intended to shed some light on this long-standing question [31][32][33] .
The data confirms the link between the amplitude of DVM and the variation in the boundary of the oxygen and the temperate zone. Figure 2A shows the seasonal variation of the sigma-theta profile (its value of 15.7 approximately corresponds to the minimal comfortable oxygen concentration of 0.4 mg/l) as well as the depth of www.nature.com/scientificreports www.nature.com/scientificreports/ the temperature level of T = 12 °C which specifies the critical boundary for comfortable living for both species preferring colder temperatures. The influence of the oxygen level and the warm waters depths on DVM can be seen more clearly in Fig. 2B,C. Figure 2B demonstrates a positive correlation between the night time depths and the profile of the temperature level of T = 12 °C for both species (Pearson correlation coefficients are = .
r 0 814 for C. euxinus and P. elongatus, respectively). Figure 2 shows a stronger correlation between the maximal depth of migration and the density sigma-theta of 15.7 (Pearson correlation coefficients are = .
r 0 841, for C. euxinus and P. elongatus, respectively). We also investigated the vertical profiles of Chl-a concentration as the index of autotrophic phytoplanktonthe primary food for zooplankton grazers-in the water column. We found, however, that phytoplankton distribution is rather variable both seasonally and from year to year. For modelling purposes we combine the data on chlorophyll and construct the annual average density of chlorophyll using our observational data for 7 years (see Fig. 2D), and do not show here the chlorophyll profile for each year and season. We fit the data with the curve which is further used in modelling; h is the depth (for details see the next section). This curve exhibits a poor fit in the surface waters, denoted by the dashed vertical line, but we can still implement the above fit for modelling purposes since individuals rarely enter the upper surface zone. www.nature.com/scientificreports www.nature.com/scientificreports/ Parametrisation of the model. We further need to parameterise the model coefficients and relate them to the daily trajectory of the individual grazer in the water column, i.e. the vertical coordinate h(t) which describes the strategy v in our original model settings. We consider that h = 0 at the water surface and that the positive direction is downwards. We re-scale time so that a 24 h period corresponds to the interval ∈ t [0, 1], and assume that due to periodicity of this function we have (1). Times t = 0 and t = 0.5 will correspond to midnight and midday, respectively. The life cycle of a zooplankton grazer consists of 6 distinct copepodite stages (see the previous section), however, here we combine some stages and only consider three age groups consisting of the youngest juvenile stages (CI-III), older juvenile stages (CIV-V) and reproductive adults (CVI). Their DVM patterns are mathematically described by the and h t ( ) A , respectively. Note that one can easily extend the approach to model all 6 developmental stages. For simplicity, we neglect the hatching time and consider = R v const ( ) and consider that the reproduction time of adults is fixed τ τ − = = T const 3 2 0 . The reproduction coefficient b(v) expressing the number of eggs produced by a female per day is given by Here W 0 is the egg carbon weight. The first term in the integrand stands for the energy gain (measured in carbon units) that zooplankton obtain from feeding on phytoplankton with the density P. This feeding rate is multiplied by an indicator function S t ( ) A , which is equal to 1, when zooplankton actively feed on phytoplankton and is zero otherwise (e.g. during migration or when staying in deep waters). We assume that grazing obeys a Holling type II law with α A and α being coefficients with the well-known ecological meaning 39 . The coefficient ε A is the food consumption efficiency.
We use the following parametrisation of the vertical distribution of phytoplankton P in the water column based on a logistic curve p 0 where P 0 gives the maximum of the phytoplankton density (we divide P 0 by 2 for normalisation); the parameter h p is the depth at which P is half its maximum. By increasing or decreasing h p we can explore the influence of food distribution on the strength of vertical migrations. Figure 2E shows that the approximation (9) can be used up to the warm surface waters, i.e. until the depths of 10-15 m. We consider that grazing by zooplankton has only a small impact on the phytoplankton and does not change its profile.
The term m h t ( ( )) Ab A describes losses due to basal metabolism (in carbon units). The metabolic rates for the considered calanoids grazers largely depend on the oxygen concentration which decreases with h and it sharply drops as an individual approaches the level of sigma-theta close to 15.7 31 . We use the following parametrisation Here h m is a characteristic depth where m Ab the depth at which the metabolic rate is half its maximum, σ m describes the sharpness of the transitional layer, and m A denotes the maximal basal metabolic cost. M h t S t ( ( )) ( )

A A A1
denotes the metabolic cost spent on active movements in the water when feeding and moving upwards during vertical migration (ascending phase). The indicator S t ( )

A1
is equal to 1 during the feeding and ascending phases of migration, otherwise it is zero. For the descending phase of DVM, we suggest that the organism only spends energy on basal metabolism. The dependence of M h t ( ( )) A A on depth is similar to that of m Ab with a different maximal value denoted by M A0 .
For each stage, the mortality rate within the day is parameterised by Here the first term is the mortality due to visual predators and it varies with the light intensity throughout the day described by t cos(2 ) 1 π − + . The coefficient γ i is the product of the constant predator density and its attack rate which is calculated at the highest light intensity (we divide γ i by 4 for normalisation). Following previous studies, we consider that the density of visual predators is fixed and the risk of predation shows a sigmoidal decrease with depth due to light attenuation 15 . However this assumption is not crucial for the our key findings. The second term stands for extra mortality near the surface and the bottom. It is parameterised as  www.nature.com/scientificreports www.nature.com/scientificreports/ maximal possible weights. The dependence of α i , m ib and M i on W is given the allometric expression (12) discussed below; ε i is the food consumption efficiency for stages = i Y J , . The indicator functions S i and S i1 for stage i has the same meaning as those in (8).
The change in the body weight is given by Eq. (4) and maturation times τ v ( ) i ( = i Y J , ) can be computed by integrating this equation. As such, the maturation time for the young juveniles τ v ( ) Y will correspond to the minimal carbon weight of juveniles W J , whereas the overall maturation time τ v ( ) J (i.e. reaching the adult stage) will correspond to the weight of an adult W A .
Note that we also need to formally set the functions S t ( ) i and S t ( ) i1 . Using the available empirical evidence we can consider that where θ is the Heaviside step function, c 0 is a critical vertical velocity of swimming. In other words, vertically moving grazers at a speed higher than c 0 do not consume food.
(here we consider that θ = (0) 1); thus this function is zero at the ascending phase and it is switched on at the other phases of DVM.
One can formally search for the exact solution h t ( ) i the optimal control problem. However, it is technically much easier to construct the optimal solution using a piecewise approximation of DVM. Interestingly, this simplification has its important biological rationale since it conforms to empirically observed cases of DVM, indicating that both the descending and the ascending velocities are fairly constant throughout DVM 40,41 . Moreover, our empirical data indicate that zooplankton remain nearly at the same depth when grazing at nighttime and while staying in deep waters. The DVM trajectories can be parameterised as follows Here H i0 and H i1 are the shallowest and the deepest depths, respectively; c i0 and c i1 the ascending and descending speeds; times t ij indicate the end of each phase of movement. By due continuity we have become algebraic ones and the equations for the optimal parameters can be obtained analytically by differentiation of fitness (7) (see SM4 for detail). One can prove that in the case where = c c i i 0 1 , the optimal DVM will be symmetrical with respect to = .
t 0 5. This allows us to reduce the number of parameters: each stage is now characterised by only 3 independent parameters.
The parameters for b v r v a v ( ), ( ), ( ) i i are taken from literature or estimated using our empirical data. They are summarized in Table 1.
The velocity of migration c i0 and c i1 is estimated as 30-50 m/h for all stages using acoustic scattering across the water column 40,41 . From the existing data, we can assume that the magnitudes of upwards and downwards migration velocities are close to each other.
The maximal mortality rate γ i due to visual predation during day is highly variable since it is a product between the predator attack rate and the predator density. In this paper, we consider it to vary as 0.6-1.0 day 1/ 14,15 . The natural mortality rate is assumed to be high for early stages (CI-CIII) (γ = . day 0 11/ Y 0 ) and negligibly small for later stages (CIV-CVI), γ ≈ 0 i0 . The thickness of transition layers σ u,d in the profiles of the temperature and the sigma-theta curve can be estimated from our data which gives, respectively, σ u = 0.05-0.21/m and σ d = 0.05-0.41/m (see SM3 for detail). The maximal values of mortality rates (given by 2δ d,u ) are hard to assess. For example, in the absence of oxygen, a copepod individual dies within less than 1 hour, which makes δ d extremely high. Here, for the sake of efficiency of numerical procedure, we consider some smaller values as δ = day 21/ d : increasing δ d does not affect the modelling results since the optimal trajectory of DVM does not enter the dangerous hypoxic zone. We also set day for the upper unfavorable zone with high temperature. This value is within the reported ranges 42 . Using our observation, we explore the following ranges h d = 50-150 m and h u = 0-30 m.
The parameters describing phytoplankton distribution (9) in the water column are estimated from our data averaged across all years (see Fig. 2E). This gives = . ± . . . To model DVM of the other zooplankton species P. elongatus, we used the parameters as in Table 1 and apply the same allometric relation for the following body and egg weights: Optimal trajectories of DVM. Numerical optimization of fitness J was done using the MATLAB function 'fminsearch' which implements the Nelder-Mead simplex algorithm of optimisation. The eigenvalues in (7) were computed numerically using the MATLAB function 'fzero' (we also found that the dominant eigenvalue is always real). To make sure that we find the global maximum, we consider different starting points for optimisation. We also checked the correctness of the optimisation results by numerically solving the equations expressing the first derivatives of fitness which should vanish at a point of maximum (see SM4 for detail). Figure 3A represents the model environment which zooplankton grazers would potentially experience in a typical summer (with = h m 140 d and = h m 20 u ): this includes profiles of their food, chlorophyll a, their daytime predation risk, the temperature and sigma-theta curves as well as variation of basal metabolic costs. Patterns of optimal DVM for the environment in Fig. 3A are shown in Fig. 3B  ). Unlike previous explanation in the literature, which suggested that zooplankton of early stages do not have enough energy to perform vertical migration, we found that they do not undergo DVM due to a high predator-independent mortality for the given stages. Our model predicts that it becomes more optimal for very early developmental stages to stay and feed in shallow waters to maximise their growth rate to be able to reach the next stages as soon as possible: stages CIV-VI generally have a lower natural mortality rate. On the contrary, the older developmental stages (CIV-VI, = i J A , ) perform night feeding in upper layers and stay in deeper waters during the day time. The obtained pattern is close to the one observed empirically (cf. Fig. 1A,B). Interestingly, stages CIV-CV and adult females (CVI) show very close DVM trajectories both in terms of migration times and migration depths. The DVM predicted by the model are triggered by daily variation of the predation load: they disappear when zooplankton mortality rate becomes constant. The choice of the optimal depths of migration is a more complicated issue and is addressed in detail below. www.nature.com/scientificreports www.nature.com/scientificreports/ Our simulation shows that the maximal depth of optimal DVM for the older stages is determined by the spatial gradient of metabolic costs rather than predation. Indeed, the depths at which zooplankton stay during the daytime are characterised by low basal metabolic rates due to low oxygen concentration ( ≈ h m 120 ), (cf. Fig. 3A). On the other hand, the consumption rate of zooplankton by visual predators becomes extremely small starting at much shallower depths ( ≈ h m 55 ) due to light attenuation and has no further effect on the maximal DVM depth. To explore the relative roles of predation and metabolism on DVM, we completed simulations using the same model parameters but assuming a depth-independent metabolic cost. We found that DVM is still observed in this case, but the migration depths during the day time become shallower. The DVM for Fig. 3C demonstrating much shallower depths of migration as compared to depth-dependent metabolic cost in Fig. 3A.
We further investigated the influence of the locations of the upper and lower unfavorable zones characterised by h d and h u on the depths of DVM denoted by H H , The results are presented in Fig. 4. In this figure, the dependence of h u on H A0 is denoted by squares 1, whereas the influence of h d on H A1 is denoted either by circles 2 (constructed for spatially variable metabolic costs) or by triangles 2′ (constructed for a hypothetic regime with spatially constant metabolic costs). As one can see, for spatially variable metabolic costs, an increase of h d starting from ( ≈ h m 60 ) will result in a pronounced increase of the migration depth. However, such an increase does not occur for the scenario with spatially constant metabolic costs. On the other hand, when the depth h u of the unfavorable zone approaches the surface, this only has a partial effect on the upper depth H A0 of DVM. This is explained by the fact that a high trophic pressure by visual predators near the surface does not allow grazers to feed in shallow waters even if the temperature regime becomes comfortable. Interestingly, these model predictions correlate well with our empirical observation shown in Fig. 2B,C.
We also investigated the dependence of DVM on other key model parameters. Varying the predation pressure on zooplankton, we found that an increase in γ i results in larger time spent by grazers in deep waters (Fig. 5A) without significantly affecting the depths of migration. At low predator pressure (e.g. γ = . day 0 1 1/ i ) migration ceases completely. An increase in the available food as described by P 0 reduces the time spent by grazers in deep waters which is shown in Fig. 5B. The model predicts that migration should cease at high phytoplankton density ( μ > P gC l 45 / 0 ). The model predicts that variation of the reproductive period T 0 within the realistic parameter range does not largely affect DVM.
Finally, we found that for both considered zooplankton species their daily food consumption and their maturation times τ i correspond well to estimates from the literature [31][32][33]48 . This provides extra credits to support our model findings.

Discussion
Fitness has always been a fundamental concept both in empirical and theoretical evolutionary biology following the seminal idea of Sewall Wright about a hypothetical adaptive fitness landscape 49 . The metaphor of climbing uphill in an evolutionary landscape to reach a local peak is so popular in the literature because it is compelling and easily graspable. However, the initial idea has been criticized because the shape of the fitness landscape in real biological systems is often non-stationary and constantly varies in the course of evolution due to a permanent dynamical feedback between individuals and their environment 50 . Here we revise the concept of evolutionary fitness using the original idea of S. Wright but allowing for dynamical feedbacks. As a result, the theory still predicts that long-term selection for a best strategy or strategies should maximise some function or functional which we understand as the fitness. The revised definition of fitness takes into account dynamical feedbacks from individuals of the same population and the environment consisting of other species and abiotic factors, and is based on introducing a ranking order of strategies, defined as the long-term ratio of generalised densities (1) of some measure in an arbitrary function space. Such densities may not have the same meaning as population densities, because they can generally be functions of the 'true' population densities (cf. SM1(ii)). We argue that the mathematical formalism proposed here is generic and can be equally applied to modelling the processes of selection and evolution in chemistry, turbulence theory, chemical kinetics, cultural evolution or economics.
The population fitness J introduced in this paper allows us to formulate the variational principle of evolutionary modelling: the strategies remaining in the system after long-term selection should correspond to the maximum of the fitness functional across the strategies initially present. An example of derivation of the variational equations providing optimal parameters is shown for the DVM model in SM4. Note that the framework allows us to consider both scalar-and function-valued strategies v: in the latter case, the variational equations will be optimal control theory equations. Establishing variational principles in modelling biological evolution has a long history with various approaches proposed 6,8,51,52 . For example, optimisation principles in evolution have been earlier suggested in the adaptive dynamics framework, where it was found that evolution will optimise the invasion fitness (the per-capita rate of growth of a rare mutant introduced into a resident population) in the case where the environment affects fitness in an effectively monotone one-dimensional manner 52 . However, generic classes of models where this property holds are still poorly understood 53 . The concept of evolutionary optimization used here is different from that used in the adaptive dynamics framework; in particular, it does not consider an invasion-replacement paradigm 52,53 . Overall, unlike the situation in mechanics, thermodynamics or optics, a unified variational principle for modelling selection in biological systems is still in its mathematical infancy and we believe that this paper makes a further step in this direction.
The proposed theoretical approach allows us to more easily define evolutionary fitness for population models with structuring (including systems with delay) which can be a challenge in other modelling frameworks. As an important practical example, we derived the expression for evolutionary fitness J for the population model of von Foerster type with age-or size-structuring with an arbitrary mortality term. This type of equation is widely used in the modelling literature and possible applications of the theory can go well beyond the considered example of Scientific RepoRtS | (2019) 9:15020 | https://doi.org/10.1038/s41598-019-51310-w www.nature.com/scientificreports www.nature.com/scientificreports/ zooplankton DVM. The derived fitness may include as many developmental stages as possible (see SM1(ii)), or we can refine the behavior or life-history within a given developmental stage by introducing sub-stages. Using a similar approach one can obtain evolutionary fitness for a structured population model with n stages under some other modelling settings: for example, by considering a system of ODEs in which each stage is described via a separate differential equation 16 .
We apply our theoretical findings to modelling optimal diel vertical migration (DVM) of zooplankton. DVM is the largest synchronised animal movement on Earth and has tremendous consequences for marine dynamics, fisheries and the ocean carbon pump 19,20 . The need for creating new models of DVM is dictated by the fact that in previous models, patterns of optimal behaviour of zooplankton were obtained by a direct maximization of some criterion such as venturous revenue, reproductive value or predation pressure [7][8][9][10]12,14,15 , and in each case the choice of criterion was based on conventional wisdom or the personal preference of the researcher. Here are we not claiming, however, that all previous models of DVM are necessarily wrong, but trying to point out the main problem of setting some initially prescribed common sense criterion of optimality. Such a criterion often does not demonstrate a straightforward connection between the maximization or minimization of some factor and the long-term population growth, especially in the presence of various nonlinear feedbacks and trade-offs. On the contrary, according to the introduced definition of fitness, the resultant patterns of optimal DVM would be a product of long-term selection of the subpopulation employing the given strategy. As such the optimal DVM is only determined by the underlying population model and the parameterisation of the model coefficients.    Table 1. www.nature.com/scientificreports www.nature.com/scientificreports/ our model confirms that the intense variation of visual predation due to the periodic change in light intensity is the main trigger of DVM, a fact which was reported in earlier DVM models. Interestingly, however, the predation pressure should be supercritical: sufficiently low pressure results in cessation of DVM. Secondly, our model predicts that predation does not entirely determine the amplitude of DVM. Namely, we find that other components of the environment such as oxygen concentration and the temperature may play a crucial role, almost tripling the amplitude of DVM compared to the case with depth-independent metabolic costs (cf. Fig. 3B,C).
The main model-based explanation of unexpectedly large amplitudes DVM in the Black Sea in summer is that staying in deep waters with low metabolic costs is beneficial for grazers since this allow them to substantially reduce their energy losses (see Fig. 3A). As a result, the deepest migration depth throughout the day can be as high as 120-130 m whereas the light intensity is reduced up to 1-2% already at depths of 50-60 m 54 . The model therefore provides a theoretical basis to explain the large amplitude DVM observed in the Black Sea, which was previously a matter of discussion in the literature. Note that the key role of metabolic costs on the amplitude of DMV has previously been suggested by direct laboratory measurements of the metabolic rates of grazers under different oxygen conditions 31 . We argue, however, that the conclusion made based on the laboratory studies 31 that staying in deep waters with low basal metabolic costs during the daytime should always be beneficial for grazers is not quite as straightforward as it seems. The model-based computation shows that fitness is maximized by staying at the low oxygen boundary only if (i) the cost of migration is not high and (ii) the drop in metabolic costs at the edge of the anoxic zone should be supercritical, otherwise large amplitude DVM becomes counterproductive. On the other hand, predation rate, natural mortality and food availability are crucial factors affecting the amplitude of DVM (see Fig. 5). These conclusions could not be made solely through empirical observation and laboratory experiments.
Another counter-intuitive conclusion from the DVM model concerns the reason that early copepodite stages (CI-CIII) do not migrate. In some previous studies 7 , this absence of migration was postulated, and has been confirmed in a large number of study cases 23 . A widespread opinion in the literature is that the main reason for the absence of DVM in early stages is the lack of sufficiently large amounts of energy-per-biomass in stages CI-CIII to be able to perform DVM. However, in our model this behaviour is an emergent property (see also Fig. 1), and we found that it is the high mortality rates rather than a lack of energy that keeps the earlier stages in shallow waters. For example, we found that reducing the natural mortality in the model results in emergence of DVM for CI-CIII even if we increase the costs of migration (we do not show this result for brevity).
Finally, we would like to give a general warning about the interpretation of the influence of several factors on DVM predicted by the model, for example, the alteration and even disappearance of DVM when the food supply in shallow waters is high enough or when predation is sufficiently low (Fig. 5). We argue that we should interpret these theoretical predictions as the optimal responses of grazers to some long-term trends rather than to fast change in the environmental conditions. Our data show that the vertical profile of chlorophyll in the considered ecosystem is highly variable from year to year, even within a single season (the data are not shown), and the same concerns the predation pressure by visual predators. A number of studies of empirical observation of DVM in the given ecosystem, however, demonstrate a constant pattern of regular migration across both days and seasons. This supports the idea that the DVM pattern of grazers observed in the model is only optimal on average and can be suboptimal on a daily basis. On the other hand, one should also consider the presence of other species dynamically. Including dynamical predators, for example, with abundance depending on that of zooplankton herbivores, may result in DVM being beneficial in ecosystems with abundant food for zooplankton 16 . In the current model settings with a single population model, this can be modeled via a trade-off between predation rates and food density, for example.
Among future perspectives for modelling selection processes based on the proposed framework we can cite the following. It would be straightforward, wherever possible, to establish classes of models where the existence of evolutionary fitness can be demonstrated. We predict that some complications may occur in the case where fitness depends on initial conditions and where the ranking order can not be properly defined. Moreover, in the case where it is not possible to analytically derive the expression for evolutionary fitness, such as in multi-species models, proper computational methods should be developed. Finally, regarding the particular ecological study case of DVM, it would be interesting to extend the considered model of zooplankton migrations to some more realistic settings by explicitly describing each of the 6 developmental stages as well as including variation in physiological conditions in each stage such as lipid reserves, moulting and starvation, as was done in some previous studies 7 . We are planning to address the above issues in our future works.

Data availability
All data generated or analysed during this study are included in this published article (and its Supplementary Material).