Evolutionary branching patterns in predator-prey structured populations

Predator-prey ecosystems represent, among others, a natural context where evolutionary branching patterns may arise. Moving from this observation, the paper deals with a class of integro-differential equations modeling the dynamics of two populations structured by a continuous phenotypic trait and related by predation. Predators and preys proliferate through asexual reproduction, compete for resources and undergo phenotypic changes. A positive parameter $\varepsilon$ is introduced to model the average size of such changes. The asymptotic behavior of the solution of the mathematical problem linked to the model is studied in the limit $\varepsilon \rightarrow 0$ (i.e., in the limit of small phenotypic changes). Analytical results are illustrated and extended by means of numerical simulations with the aim of showing how the present class of equations can mimic the formation of evolutionary branching patterns. All simulations highlight a chase-escape dynamics, where the preys try to evade predation while predators mimic, with a certain delay, the phenotypic profile of the preys.


1.
Introduction. Selective pressures exerted by the surrounding environment can make a population initially composed of individuals expressing the same phenotype (i.e., a phenotypically monomorphic population) to split into two or more distinct populations, or phenotypic clusters. Such a dynamical process may occur several times along the history of a living specie leading to the formation of evolutionary branching patterns.
Predator-prey ecosystems represent, among others, a natural context where evolutionary branching patterns may arise. In fact, preys are parts of the predators' environment and, in turn, predators shape the environment around the preys, so that a change in the relative distribution over the phenotypic traits of the preys can induce a change in the predominant traits within the predators' population. For instance, one can dynamically observe the selection for those preys that are able to escape predation and those predators whose phenotypic traits allow them to catch a large number of preys.
The derivation of models able to reproduce branching patterns has fascinated mathematicians to such an extent that several different papers dealing with this topic, both from an analytical and a numerical point of view, have been published.
For instance, see [2,8,9,10,11,15,18,22] and references therein. In this work, we focus on a set of coupled integro-differential equations modeling the dynamics of two populations structured by a continuous phenotypic trait and related by facultative predation, i.e., the growth of predators may require the consumption of preys. Individuals proliferate through asexual reproduction, compete for resources and may give birth to offspring expressing their same trait (i.e., renewal) or a slightly different one (i.e., mutations) [6,7,13]. A positive parameter ε is introduced to model the average size of phenotypic changes.
Taking advantage of the approach used in [5,6], we study the asymptotic behavior of the solutions of the mathematical problem linked to the model in the limit ε → 0 (i.e., in the limit of small phenotypic changes). Numerical simulations are performed with the aim of illustrating and extending analytical results as well as highlighting the ability of the model to mimic the formation of evolutionary branching patterns.
It is worth noting that there is a substantial difference between the patterns here presented and the ones that we obtained in [5]. In fact, there the fittest traits (i.e., the ones that are selected over time) could be, in a certain sense, foreseen a-priori, while here they emerge from the interactions among the two populations and their evolution is, at least to a certain extent, not predictable. Therefore, the branching patterns reproduced by the present model are closer to the ones that can be found in Nature.
The contents of this work are organized into three further sections as follows: -Section 2 introduces the class of equations under consideration and the related underlying assumptions. -Section 3 presents the initial value problem of interest, recalls its well-posedness and studies the asymptotic behavior of its solution in the limit ε → 0. -Section 4 summarizes the results of some numerical simulations, which illustrate and extend the asymptotic results established by Section 3 and depict the formation of branching patterns.
Notations given hereafter are used throughout the paper: -Defining U as a compact subset of R + , on account of brevity, · 1 stands for the L 1 -norm over U , while · ∞ stands for the L ∞ -norm over U × U . -Given h 1 , h 2 : R + × U → R, the set X + is the positive cone of the Banach space 2. The model. This section aims at presenting the class of equations defining the model under consideration as well as at summarizing the related underlying assumptions.
The model here considered is defined by the following system of integro-differential equations, where function f (t, u) ≥ 0 stands for the density of preys expressing a given phenotypic trait, which is modeled by the independent variable u ∈ U := [0, 1], while function g(t, u) ≥ 0 is the density of predators whose phenotype allows them to mainly catch those preys that express the trait u: mutations and renewal proliferation and competition proliferation and competition Independent variable t is the time variable belonging to the interval (0, T ], with an arbitrary T > 0.
The following notations and assumptions hold: -Kernel K(u − u * ; ε) models the probability that an individual with trait u * gives birth to individuals with trait u and it is such that We make the assumption that mutations are rare, so that they preserve, net of renewal, the total number of individuals inside the system. Furthermore, since we allow only small variations in the phenotypic trait to occur from parent to offspring, we introduce a small parameter ε such that K is negligibly small for u outside an ε-neighborhood of u * . In particular, we define kernel K as follows: where δ is the Dirac's delta distribution and parameter α ∈ R + describes the average probability for mutations. -Parameters κ 1 , κ 2 ∈ R + describe, respectively, the average rate of preys' and predators' proliferation, net of death by natural causes, due to the consumption of resources offered by the surrounding environment. Parameter κ 2 embodies the facultative nature of the predation under consideration. -Parameters µ 1 , µ 2 ∈ R + model, respectively, the average rate of preys' and predators' death due to intra-population competition. -Function η(|u * − u|) (or η(|u * − u|)) stands for the rate of interactions leading a prey expressing the phenotypic trait modeled by u (or u * ) to be caught by a predator that is mainly able to catch those preys that express the trait identified by u * (or u), so that the prey dies, and thus the predator proliferates. As a result, we make the following assumption to hold: In order to capture the slow evolutionary process leading to substantial changes in the predominant traits, we can use the same time rescaling introduced in [7] and rewrite Eqs. (1) as follows: In the following section, we are going to analyze the behavior of functions f ε and g ε in the limit ε → 0. Such asymptotic analysis are meant to verify the weak-convergence to some limit measures f and g and, in particular, to identify their supports. In fact, as it has been highlighted in [6], the points belonging to the supports of such measures are equivalent to the Evolutionary Stable Strategies studied by the Evolutionary Game Theory [14] and thus, within this mathematical context, looking for branching patterns calls attention to the location of these points.
3. Analytical results. This section is meant to provide analytical properties for the Cauchy Problem derived by endowing Eqs. (1), or Eqs. (4), with suitable initial conditions. In more detail, we recall two theorems ensuring the well-posedness of such initial value problems and prove a third theorem characterizing the asymptotic behavior of the related solutions in the limit ε → 0.
We focus on the Cauchy Problem for functions f and g given hereafter, which can be derived by endowing Eqs. (1) with suitable initial conditions: where h(t, u) = (f (t, u), g(t, u)), h 0 (u) = (f 0 (u), g 0 (u)) and operator H[h] is componentwise defined by Eqs. (1).
Problem (5) is well-posed in the sense of Hadamard (i.e., the solution exists, it is unique and depends continuously on the initial data), as stated by the following theorems, which can be proved by standard fixed point arguments and so are left without proof: Theorem 3.1 (Well-posedness and non-negativity). For any h 0 ∈ X + , Problem (5) admits a unique local in time solution h ∈ C([0, T ], X + ) that satisfies where a 0 and T are positive constants depending on the initial data as well as on the model parameters. (5), which satisfies the following a priori estimate where C T is a given constant depending on T , on the initial data and on the parameters of the model.
In particular, the proof of Theorem 3.2 relies on the a priori estimates given hereafter, which will also be used in the asymptotic analysis developed in the following: Lemma 3.1. Let f (t, u) and g(t, u) be the components of the unique solution of Problem (5). Then, the following upper bounds hold true: for all t ≥ 0.
Proof. Integrating over U Eqs. (1) and making use of assumptions (2) and (3), as well as of the non-negativity of functions f and g, we achieve the following inequalities which imply: ∀t ≥ 0. This concludes the proof.
Let h(t, u) be the unique solution of Problem (5), whose existence is guaranteed by Theorem 3.1 and Theorem 3.
where H[h ε ] is componentwise defined by Eqs. (4). The results above stated for Problem (5) apply to Problem (6) as well, for any ε > 0. In particular, Problem (6) is well-posed, its unique solution is non-negative and satisfies componentwise the upper bounds established by Lemma 3.1. Now let us introduce the following definitions: Thus, the components of H ε (t, u) can be written as follows and the following preliminary result holds: Let f ε (t, u) and g ε (t, u) be the components of the unique solution of Problem (6). Then, for all smooth test function ϕ belonging to the completion of for any ε > 0 and t ≥ 0.
Proof. The definition of M 1ε (t, u) for u ∈ (ε, 1 − ε) implies that for any smooth test function ϕ ∈ C ∞ c ((ε, 1 − ε)). Hence, the following identity holds true: since t 0 f ε (s, ·)ds is bounded for any t ∈ [0, T ] and, in the limit ε → 0, the difference quotient approximates the second derivative of ϕ. By approximation, this also holds for the completion of C ∞ c ((ε, 1 − ε)) in W 2,∞ (U ). Since e The result related to g ε (t, u) can be proved in the same way.
Due to Lemma 3.2, the asymptotic behavior of f ε (t, u) and g ε (t, u) for ε → 0, i.e., in the limit of large times and small mutations, can be characterized by means of techniques similar to those that have been introduced in [6], as highlighted by the following Theorem 3.3 (Large times and small mutations asymptotics). Let the following additional assumptions hold: Then, there exist a subsequence of f ε (g ε ), denoted again as f ε (g ε ), and a subsequence of R 1ε (R 2ε ), denoted again as R 1ε (R 2ε ), such that: i) Establishing convergence.
where: f, g ∈ L ∞ ((0, T ), M 1 (U )) and ii) Characterizing the support of the limits f and g.
The same reasoning can be developed for R 2 (t, u) and g(t, u). In fact, assume that there existst s.t.
Definitions (10) imply where only the third terms of the right hand sides depend on u. Then, making also use of assumption (8), we can draw the conclusion given hereafter: Remark 3.1. From a biological perspective, the first one of assumptions (9) translate into mathematical terms the idea that proliferation of predators due to the intake of other resources occurs at a lower rate than proliferation of preys. On the other hand, the second one mimics a scenario where competition among predators occurs at a higher rate compared to predation. We focus on the initial value problem defined by endowing Eqs. (1) with initial conditions: where A 0 and B 0 are positive real constants such that assumptions (8) are satisfied. Since f (t, u) stands for the density of preys expressing a given phenotypic trait u, these initial conditions correspond, respectively, to a scenario where the preys population is initially composed of individuals that mostly express a single phenotype (i.e., a monomorphic population), or two different phenotypes (i.e., a dimorphic population). As predators and preys share the same initial distribution, predation can occur from the beginning of observations. Throughout all the simulations, parameters related to mutations (i.e., the average size of phenotypic changes ε and the average mutation rate α) are defined as ε := 0.001, α := 0.5.
In fact, we consider an evolutionary scenario where phenotypic changes are small. On the other hand, function η is defined as: so that assumptions (3) are satisfied. This definition embodies the selectivity of the predator-prey interactions under consideration. At first, we aim at illustrating the analytical results about the asymptotic behaviors of functions f (t, u) and g(t, u) in the limit ε → 0 established by Theorem 3.3. We set the κ-type parameters related to proliferation because of the consumption of resources and the µ-type parameters related to death caused by intra-population competition as follows: in agreement with assumptions (9). The numerical results summarized by Figures 1-4 illustrate the analytical ones presented in the previous section and highlight, through the selection of different phenotypic traits, the formation of evolutionary branching patterns. Let us note that, as expected on the basis of analytical results, at time t, function f is mainly concentrated in the minimum points of I 1 (t, u), while function g concentrates in the maximum points of I 2 (t, u). Functions I 1 and I 2 are defined by (11).
From a biological perspective, Figure 1-4 highlight how, starting from few phenotypic clusters, the occurrence of mutation, proliferation and competition phenomena under the effects of predation can cause the emergence of new different traits, so leading the phenotypic expression of preys and predators to become more heterogeneous over time.  (16), function η defined in (19) and parameter setting (18) and (20). Right panel: trend of I 1 (t, u), defined in (11), at t = 1000dt under the same parameter setting.
After that, in order to extend the results established by Theorem 3.3, we perform simulations setting the κ-type and µ-type parameters as follows or as so that the assumptions of Theorem 3.3 are not satisfied. Figure 5 and Figures 6-7 refer to parameter setting (21) and (22), respectively, and show that branching patterns arise even if the assumptions of Theorem 3.3 are not satisfied.
In more detail, making a comparison between the numerical results summarized by the previous figures and the ones here presented in Figure 5, we are led to conclude that the dynamics of g(t, u) can be slowed down if κ 2 is chosen sufficiently higher than κ 1 (i.e., if the proliferation of predators is sufficiently faster than the  (16), function η defined in (19) and parameter setting (18) and (20). Right panel: trend of I 2 (t, u), defined in (11), at t = 1000dt under the same parameter setting.  (17), function η defined in (19) and parameter setting (18) and (20). Right panel: trend of I 1 (t, u), defined in (11), at t = 1000dt under the same parameter setting.
one of the preys). In fact, this makes the environment to be harsher for predators and cause their phenotypic profiles to evolve more slowly. On the other hand, the right panel of Figure 7 shows that the solution of the initial value problem at hand is characterized by oscillations in the population sizes if the net proliferation rate of predators is negative and no death due to intra-population competition  (17), function η defined in (19) and parameter setting (18) and (20). Right panel: trend of I 1 (t, u), defined in (11), at t = 1000dt under the same parameter setting. occurs (i.e., if assumptions (22) are satisfied), while the assumptions of Theorem 3.3 prevent oscillatory trends, at least on the same time interval (see left panel). The arising of oscillations in the total densities of the two populations can be prima facie justified by noting that, under assumptions (22), definition (19) makes the evolution equations for 1 (t) and 2 (t) to resemble a kind of Lotka-Volterra system in the regime of small ε, save for a proper time line rescaling.
It is worth noting that all simulations highlight a dynamical adaptation of the phenotypic profile of the predators to the one of the preys. In particular, chaseescape dynamics arise, since the preys try to evade predation while predators mimic, with a certain delay, the phenotypic profile of the preys to make their attack more effective.
These chase-escape dynamics can be interpreted, at least to a certain extent, as collective learning processes. Collective learning is a self organization ability playing a crucial role in the evolutionary dynamics of several complex living systems, both in the biological and in the socio-economic context. For instance, making reference to cancer dynamics [20], cells belonging to the immune system evolve, namely through learning and adaptation, to become able to recognize and effectively attack cancer cells; on the other hand, the selective pressure exerted by the immune system against cancer cells leads, through a process called immunoediting [5], to the selection of those clones that are able to evade the predation exerted by immune cells [12]. In opinion formation [21] or in financial markets [16], the agents have to learn from information sources and to adapt to a changing environment where, in some cases, the option selected by the minority part of the agents can be the winning strategy (see minority games [3]). Finally, in crime dynamics [17], police seek for criminals who try resist the attacks and escape.