# The origin of species by means of mathematical modelling

Nikolai Bessonov, Natalia Reinberg, Malay Banerjee, Vitaly Volpert

► To cite this version:

Nikolai Bessonov, Natalia Reinberg, Malay Banerjee, Vitaly Volpert. The origin of species by means of mathematical modelling. Acta Biotheoretica, 2018, 66 (4), pp.333-344. 10.1007/s10441-018-9328-9. hal-01941452

HAL Id: hal-01941452

https://hal.science/hal-01941452

Submitted on 6 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The origin of species by means of mathematical modelling

Nikolai Bessonov · Natalia Reinberg ·
Malay Banerjee · Vitaly Volpert

**Abstract** Darwin described biological species as groups of morphologically similar individuals. These groups of individuals can split into several subgroups due to natural selection, resulting in the emergence of new species. Some species can stay stable without the appearance of a new species, some others can disappear or evolve. Some of these evolutionary patterns were described in our previous works independently of each other. In this work we have developed a single model which allows us to reproduce the principal patterns in Darwin's diagram. Some more complex evolutionary patterns are also observed. The relation between Darwin's definition of species, stated above, and Mayr's definition of species (group of individuals that can reproduce) is also discussed.

**Keywords** nonlocal reaction-diffusion equations · travelling waves · pulses · Darwin's diagram

N. Bessonov

Institute of Problems of Mechanical Engineering, Russian Academy of Sciences, 199178 Saint Petersburg, Russian Federation

N. Reinberg

Institute of Problems of Mechanical Engineering, Russian Academy of Sciences, 199178 Saint Petersburg, Russian Federation

M. Banerjee

Department of Mathematics and Statistics, IIT Kanpur, Kanpur, India

V. Volpert

Institut Camille Jordan, UMR 5208 CNRS, University Lyon 1, 69622 Villeurbanne, France
INRIA Team Dracula, INRIA Lyon La Doua, 69603 Villeurbanne, France

Peoples Friendship University of Russia (RUDN University, 6 Miklukho-Maklaya St, Moscow, 117198, Russian Federation)

Poncelet Center, UMI 2615 CNRS, 11 Bolshoy Vlasyevskiy, 119002 Moscow, Russian Federation

E-mail: volpert@math.univ-lyon1.fr
# 1 Evolution of species and Darwin's diagram

Emergence and evolution of biological species continue to attract attention from both the biological point of view and from the modelling point of view. Modern theory distinguishes three main mechanisms of speciation. Allopatric speciation implies the existence of geographic or genetic isolation where there is no gene exchange between different taxa; parapatric speciation admits partial exchange; and sympatric speciation occurs without geographic or genetic barriers [6], [11]. It is generally accepted that allopatric speciation is biologically realistic, though the mechanism which leads to the appearance of isolating barriers may sometimes be unclear [6] (Chapter 2). Sympatric speciation continues to instigate intensive discussions. In spite of the big body of experimental data, observations in nature, and theoretical models, this question remains controversial [6] (Chapter 4). Various models of sympatric speciation are developed in [8], [9], [10], [12], [13], [14], [15].

Competition for resources and selection of the most successful in this competition are among the fundamental properties of biological systems. Once the elements of such a system begin to compete, they try to get as many resources as possible. For this purpose they consume resources not only at the spatial location, where they are located, but also in some areas around it. This results in nonlocal consumption of resources and intraspecific competition. Nonlocal consumption of resources, reproduction, and variation of functional traits (diffusion) determine the process of speciation [13].

In this work we will address the question about the emergence and evolution of biological species from the point of view of Darwin's diagram. It is the only figure featured in his famous book "The origin of species" [7], and it is the quintessence of his theory, which corresponds, in the modern terminology, to sympatric speciation. We will attempt to give a more precise mathematical interpretation of his theory and to point out the assumptions which are necessary in order to reproduce Darwin's diagram. He described this figure as follows. "Let *A* to *L* represent the species of a genus large in its own country; these species are supposed to resemble each other in unequal degrees, as is so generally the case in nature, and is represented in the diagram by the letters standing at unequal distance ... The little fan of diverging dotted lines of unequal length proceeding from (*A*), may represent its varying offspring. The variations are supposed to be extremely slight, but of the most diversified nature, they are not supposed all to appear simultaneously, but often after long intervals of time; not are they all supposed to endure for equal periods. Only those variations which are in some way profitable will be preserved or naturally selected. And here the importance of principle of benefit being derived from divergence of character comes in; for this will generally lead to the most different or divergent variations (represented by the outer dotted lines) being preserved and accumulated by natural selection. ... After a thousand generations, species (*A*) is supposed to have produced two fairly well-marked varieties, namely *a*¹ and *m*¹. These two varieties will generally continue to be exposed to the same conditions which made their parents variable, and the tendency to variability
is in itself hereditary, consequently they will tend to vary, and generally to vary
in nearly the same manner as their parents varied. Moreover, these two vari-
eties, being only slightly modified forms, will tend to inherit those advantages
which made their common parent (A) more numerous than most of the other
inhabitants of the same country” ([7], p. 104).

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Left: Darwin's diagram explaining the emergence and evolution of biological species (adapted from [7]). Right: A typical population density distribution with respect to a morphological characteristics. The coordinates of the peaks of such distributions on the (*x*, *t*)-plane corresponds to the line in Darwin's diagram.</figcaption></figure>

Thus, each line corresponds to one species or subspecies. The horizontal
axis measures their resemblance; that is, some morphological characteristics
or phenotype. The vertical axis is time. Different species are represented as
lines (trajectories) in the (*x*, *t*)-space.

The appearance of new species is shown as several lines starting from the
same point. In the beginning, these are small variations inside the same species.
After some time they can become different species. Among many possible
variations, only favorable ones are chosen and reinforced by natural selection.
The others will disappear. One more important feature of the diagram is that
the lines (species) can disappear when two of them approach each other. We
can suppose that this happens because the corresponding species begin to
compete for the same resources.

We will model the main elements of Darwin’s diagram: existence of species
as localized groups of individuals, their evolution (motion in the morphological
space), emergence of new species, preservation of the most fitted species and
disappearance of others due to natural selection, and competition of species
for resources. We will introduce the model in the next section and will analyze
it in Section 3. We will then discuss the results in Section 4.
## 2 Model

We consider the diagram in Figure 1 (left) from the bottom upwards (in the direction of growing time). It begins with 11 lines noted by the letters A, B, ..., L. We will mainly consider the left part of the diagram (A, B, ..., F), though some other parts will also be mentioned when it is appropriate.

Following Darwin's definition, we understand biological species as groups of morphologically similar individuals. A typical population density distribution for one species is shown in Figure 1 (right). The variable $x$ here is a morphological characteristic or a functional trait of the population (height, weight, beak size, and so on) and the function $u(x, t)$ describes the density distribution with respect to this characteristic at some given moment of time $t$. The species are characterized by narrow distributions, such that they can be described by the coordinate $x_m$ of the maximum of the function $u(x, t)$. This coordinate $x_m(t)$, considered as a function of time, represents a line on the $(x, t)$ plane. Darwin's diagram corresponds to the ensemble of such lines for different species.

Thus, our aim is to develop a model such that a) it has solutions in the form of pulses (or peaks), which describe population density distributions, and b) the coordinates of the maxima of these pulses on the $(x, t)$-plane form the patterns qualitatively similar to those in the diagram.

In [5] we reproduced some parts of the diagram independently of each other. In this work we will develop a new model which will allow us to describe these different parts together and to reproduce the left half of the diagram. We consider 6 species in accordance with the letters A, B, ..., F and introduce the corresponding population densities $u_i(x, t)$, $i = 1, 2, ..., 6$. In order to describe the desired behavior of solutions, we consider the integro-differential system of equations:

$$ \frac{\partial u_i}{\partial t} = d_i \frac{\partial^2 u_i}{\partial x^2} + a_i u_i^2 W_i - \sigma u_i, \quad i = 1, 2, ..., 6. \quad (1) $$

Here, the space variable $x$ corresponds to some morphological characteristics of the populations, $t$ is time, and $d_i$, $i = 1, ..., 6$ are the diffusion coefficients, $a_i$ and $\sigma$ are some positive constants. The first terms in the right-hand side of these equations describe small variations of the morphology due to random mutations, the second terms describe the reproduction of the population and the last terms describe species' mortality. Let us note that the reproduction terms are proportional to the square of the population density (Supplementary materials) and to available resources $W_i$, considered in the following general form:

$$ W_i = (1 - J_i(u))(1 - I_i(u)), \quad i = 1, ..., 6, $$

where

$$ J_i(u) = \sum_{j=1}^{6} b_{ij} J_j^0, \quad I_i(u) = \sum_{j=1}^{6} s_{ij} I_j^0, $$
$$J_j^0 = \int_{-\infty}^{\infty} \phi(x-y)u_j(y,t)dy, \quad I_i^0 = \int_{-\infty}^{\infty} u_i(y,t)dy.$$

Here, $b_{ij}$ are non-negative numbers. The coefficients $s_{ij} \ge 0$ depend, in general, on the space variable $x$, though we will begin below with the simpler case of constant coefficients. The kernel $\phi(x)$ of the integral determines the efficacy of this consumption as a function of the distance $x-y$. We take it as a piece-wise constant function,

$$\phi(x) = \frac{1}{2N} \begin{cases} 1, & |x| < N \\ 0, & |x| \ge N \end{cases}.$$

In the limit of small $N$, we get the case of local consumption. In the case of global consumption, the consumption of resources is proportional to the total population density (integrals $I_i^0$) [5], [16].

In conventional models of population dynamics, available resources are considered in the form $(1 - u_i(x, t)/K)$, where $K$ is carrying capacity. Consumption of resources is proportional to the population density $u_i(x, t)$ taken at the space point $x$. In the case of nonlocal consumption of resources, it is replaced by the integral $J_i(u)$ [13]. This implies that the individuals of the population consume resources in some areas around their average location. If there are several species, then all of them (or only some of them) can contribute to the consumption of some particular resources.

The model considered in this work contains two resources: one with nonlocal consumption and another with global consumption. We will see below that both are necessary in order to describe the patterns shown in the diagram. It is a generic model of the interaction of 6 species. It generalizes the previously considered models with local and nonlocal consumption of resources (see Supplementary materials).

### 3 Modelling of the diagram

Darwin's diagram illustrates different evolution patterns, such as speciation, survival of the better adapted, competition of species, among others. In the model presented above, we set $b_{11} = b_1, b_{1j} = 0, j = 2, ..., 6, b_{ij} = b_2, i = 2, ..., 6, j = 1, ..., 6, s_{11} = s_1, s_{ii} = s_2, i = 2, ..., 6, s_{ij} = 0$ for $i \neq j$. These assumptions on the parameters signify that the first equation for the variable $u_1(x, t)$ is independent of the other equations but it does influence the other variables. Other values of parameters will be specified below.

*Speciation.* The main issue of Darwin's theory is the emergence of new species. This is shown in the diagram on many occasions: points $A, a^1, a^2, ..., m^1, m^2, ..., I, z^1, z^2, ...$. For each point, there are several lines starting from a single point. These lines do not correspond to fully formed species, but only to small variations in the existing species, which can lead to new species in the future.
Among the species A, B, ..., F, only the first one manifests speciation and leads to the appearance of new species. The corresponding function $u_1(x, t)$ satisfies the equation

$$ \frac{\partial u_1}{\partial t} = d_1 \frac{\partial^2 u_1}{\partial x^2} + a_1 u_1^2 (1 - b_1 J_1^0(u_1))(1 - s_1 I_1^0(u_1)) - \sigma u_1. \quad (1) $$

In order to simplify the model and the analysis of the results, we chose this equation to be independent of other variables. Clearly, similar results can be obtained if the dependence on other variables is weak enough.

As discussed above, we characterized solutions of this equation by the positions of their maxima at the $(x, t)$-plane. They correspond to the lines (species) in Darwin's diagram. A typical branching pattern obtained in numerical simulations of equation (1) is shown in Figure 2 (the left part of the graph). The initial condition for the first species, $u_1(x, 0)$, is a piece-wise constant function different from zero in a narrow interval at the center of the computation domain. This initial population splits into two sub-populations (Figure 2, right), then again and again, expanding in space as a periodic travelling wave and leaving stationary localized peaks behind it. These peaks correspond to the new species.

Speciation can be described by the model with nonlocal consumption of resources (see Supplementary materials) and by the model with two resources: nonlocal and global. The question about the possibility of sympatric speciation was intensively debated in the literature, and various mathematical models were developed [1], [8], [9], [10], [13], [14].

<figure><img src="image_3.png" /><figcaption>**Fig. 2 Left:** Numerical simulations of system (1) with $d_i = 0.01, a_i = 1, b_1 = 1, b_2 = 20, s_1 = 0.01, s_2 = 0.1, \sigma = 0.01, N = 5$. The lines show positions of maxima of each component of the solution on the $(x, t)$-plane. The first species (variable $u_1(x, t)$) shows a typical branching pattern, while the other species do not change until they meet the descendants of the first species and then become extinct. Different colors correspond to different components $u_i(x, t)$ of the solution. **Right:** Initially the population density had unimodal distribution. After some time, it split into two subpopulations (reprinted with permission from [13]).</figcaption></figure>

*Survival of the better adapted.* We observe in the diagram that among several new sub-species only one or two survive, while the others disappear after a
relatively short time. The model of speciation with nonlocal consumption does not describe this effect: once a species appears, it persists. This contradicts Darwin’s description that only better adapted species survive. Usually, the outer lines (species) persist, while the species located inside become extinct (cf. Section 1).

If we accept the assumption that some species (or sub-species) are better adapted than the others, then their natality and mortality rates depend on the phenotype, that is on the space variable $x$. Hence the coefficients of the equations depend on the space variable. This approach was suggested in [5]. We will use it here in a more complete model.

We illustrate modelling with space-dependent coefficients for another part of the diagram starting with $z^6$ [5]. Figure 3 shows an example of numerical simulations of equation (1) with a space dependent function $s_1(x)$. The initial condition is a narrow pulse located at the center of the interval. Due to specification, it gives six pulses. The outer two of these pulses reach the points where the function $s_1(x)$ decreases (step-wise constant function in the same figure). As a consequence, these two outer pulses begin to grow while the others disappear. A similar behavior repeats several times at each step of the function $s_1(x)$. We obtain a specific V-pattern, which is one of the main repeating patterns in the diagram. Similar patterns start at the points $A, a^3, a^5, m^4, z^4, z^6$.

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Modelling of the part of the diagram which begins at the point $z^6$ with equation (1). The green lines show the trajectories of the maxima of the function $u(x, t)$, the red curve is a snapshot of solution, and the violet line shows the function $s_1(x)$. The values of parameters are as follows: $d_1 = 0.1, a_1 = 20, b_1 = 0.5, s_1 = 0.018$ (lowest value of the function $s_1(x)$), $\sigma = 0.1, N = 3$.</figcaption></figure>

Let us emphasize the total size of the population is limited due to global consumption. Therefore, if one of the species grows in size, then the others will decline. This is the mechanism which ensures development of the most adapted to the prejudice of the others. If we consider for example only non-
local consumption, where the size of the population is not limited, then this
mechanism does not work.

*The absence of speciation.* Species *B*, *E*, ..., *F* persist without speciation. Such behavior is observed in the case of global consumption (see Supplementary materials). Hence, the behavior of the variables *u*₂, ..., *u*₆ is primarily determined by global consumption while nonlocal consumption takes into account competition with other species (the next paragraph).

*Competition of species.* Another repeating pattern in the diagram is extinction of one of the species when two of them approach each other. Among them, species *B*, *C* and *D* disappear when they meet the line *m* descendant from *A*. It also happens for *G* and *u* and in several other cases. Thus, we need to describe the situation where two species coexist if they consume different resources and one of them disappears when they begin to compete for the same resources. Since we associate species to the coordinates of the maxima of the corresponding density distributions, two species do not compete if these coordinates are different, and they begin to compete when they are close to each other.

Thus we need to describe a) the displacement of species in the morpholog-
ical space, b) their competition when two species (coordinates of the maxima)
are close and the absence of competition when they are far from each other
in this space, c) the disappearance of some of the species due to this competi-
tion. Clearly, these properties are not described by the conventional model of
competition of species which does not take into account their location in the
space of phenotypes.

Let us begin with the case where the coefficients $s_i$ are constant. The first
component $u_1$ of the solution expands in space and begins to influence other
components of the solution (Figure 2). First, it meets the second component
$u_2$, which is bounded from the other side by the third species and is rapidly
going to extinction. Similarly, the third, fourth and fifth species disappear one
after another when pushed from one side by the expanding first species and
from the other side by the remaining species. The last, sixth species persists
during the time interval shown in the diagram.

The introduction of space dependent coefficients allows us to take into
account the survival of the better adapted discussed above. Figure 4 (upper
left) shows the initial conditions and the function $s_1(x)$. The results of the
numerical simulation are shown in Figure 4 (upper right).

*Time-dependent coefficients.* In order to finish modelling Darwin's diagram, we need to introduce time dependent coefficients. Let us explain why it is necessary. The pulse moves in such a way that it tends to improve living conditions of the population: decrease the mortality rate, increase natality and the total size of the population. These conditions determine the direction of its motion, and it cannot move in the opposite direction because it would be unfavorable for the population. These observations are in agreement with
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Numerical simulations of system (1) with $d_1 = 0.2, d_2 = ... = d_6 = 0.01, a_1 = 20, a_2 = ... = a_6 = 1, s_2 = ... = s_6 = 0.1, b_1 = 0.5, b_2 = 20, \sigma = 0.1, N = 3$. The function $s_1(x)$ (magenta) changes in time and is shown for each time interval. The upper left image shows the initial conditions for all components of a solution and the initial form of the function $s_1(x)$. The result of this simulation at $t = 167$ is shown in the upper right image. The snapshots of solutions for $t = 265$ and $t = 375$ are shown in the lower figures. Each of the figures shows the form of the function $s_1(x)$ used at the next stage of the simulation.</figcaption></figure>

the biological meaning of the considered models, and they are confirmed by all obtained results.

Consider the line from the point A to $a^5$ of the diagram which moves to the left. However, the line *f*, which starts at the point $a^5$, turns to the opposite direction. If the parameters of the model remain the same, this should not happen because either the line *a* or the line *f* would deteriorate the life conditions of the corresponding species. Similar behavior is shown in other parts of the diagram ($m^4, z^4, z^6$). Hence we need to conclude that at the moment when the curve *a* reaches the point $a^5$, the coefficients of the equation change. Biologically, this means that the new species possess new properties which are reflected through the corresponding changes in the coefficients of the model.

Figure 4 shows the results of numerical simulations of system (1). The initial functions $u_i(x, 0)$ and the initial function $s_1(x)$ are shown in the upper left image. The result of this first stage of the simulation is shown in the upper right image together with the function $s_1(x)$, which is used in the next stage. Further results and the function $s_1(x)$ are shown in the other images. The final result is shown in Figure 5 (left). It qualitatively repeats the main features of this part of the diagram (Figure 5, right).
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Darwin's diagram (right) and the numerical simulations (left) show its main patterns: coexistence of several species, V-pattern, extinction of species due to competition, time dependent coefficients allowing the species *f* to return back in the space of phenotypes.</figcaption></figure>

## 4 Discussion

Darwin's diagram represents the essence of his theory, and it contains several different evolutionary patterns. First of all, the emergence of new species, the question largely discussed in the modelling literature, and also several other models including survival of the most adapted, existence of species without speciation, and competition of species. Even the question of existence of species, which might seem obvious, appears to be quite involved as far as mathematical models are concerned. All these questions stimulated the development of the new models that were brought together in this work. We discuss some of these below as well as some further possible developments of these models.

*Existence of species.* Biological species correspond to localized population density distributions in the form of pulses. Reaction-diffusion equation (1) (Supplementary materials) with a bistable nonlinearity has such stationary solutions. However, they are unstable. They become stable in the case of global consumption of resources (3). The bistable case corresponds to sexual reproduction where the reproduction rate is proportional to the square of the population density $u^2$. The existence of such solutions is also proved for the system of two equations with different densities of males and females [17]. In the monostable case, where the reproduction rate is proportional to the first power of *u*, single stationary pulses do not exist. However, there are periodic travelling waves which can leave multiple pulses after propagation. Thus, the existence of biological species, considered as particular types of solutions, imposes certain conditions on the model and on the parameters.

*Emergence of species and survival of most adapted.* One of the main issues of Darwin's theory concerns sympatric speciation. In our previous work [13] we have developed a minimal reaction-diffusion model that describes this process.
It includes nonlocal consumption of resources, which leads to intraspecific competition, reproduction, and small random mutations. In the case of sexual reproduction, the phenotypes, or the functional traits, related to consumption of resources of males and females are supposed to be the same. Their difference can influence the emergence of species (see the next paragraph).

The minimal model of the emergence of species is not sufficient to describe some properties of this process stressed by Darwin; namely, survival of the most adapted subspecies and disappearance of the others. In order to describe this effect, we introduce a model with two resources: nonlocal and global. It is interesting that, similar to Darwin’s description, we observed that “the most different or divergent variations” survive while the others go to extinction.

*Darwin versus Mayr.* Darwin considered biological species as groups of morphologically similar individuals. According to Mayr, these are groups of individuals that can reproduce among themselves and have fertile offspring. We will briefly discuss here how these two definitions are related to each other from the point of view of the modelling presented in this work.

The reproduction terms in equations (1) are proportional to the second power of the population density. This term appears in the case of sexual reproduction where the reproduction rate depends on male and female densities which are supposed to be equal to each other (see Supplementary materials). Therefore, it is implicitly assumed that males and females having offspring possess the same functional trait that is responsible for consumption of resources.

Suppose now that males and females having common offspring can possess different phenotypes. Then, instead of $u^2$, we have the nonlocal term $u(x,t) \int_{-\infty}^{\infty} \psi(x-y)u(y,t)dy$, where the kernel $\psi(x-y)$ shows the probability of phenotypes $x$ and $y$ to have common offspring. This model was studied in [4], and it was shown that the new species will appear if the support of $\psi$ is small enough, such that the phenotypes $x$ and $y$ with common offspring belong to the same localized solution (species). If we allow inter-breeding, that is, the support of $\psi$ is larger than the distance between the species in the space of phenotypes, then new species do not appear.

Hence, Mayr’s condition is necessary for the emergence of species that are understood in the sense of groups of morphologically similar individuals. Indeed, if individuals from different species can have common offspring, then the groups of morphologically similar individuals will not emerge.

*Multi-level speciation.* Darwin’s diagram contains many important elements of his theory. However, it is only an illustration, and many other cases and generalizations are possible. We will present here one of them in order to show an example of multi-level speciation.

Biological species are organized in genera, families, and so on. In the modelling presented above, we considered only species and not higher levels of this hierarchy. In order to introduce two levels of speciation, we need to consider
double nonlocal consumption of resources with two different kernels. In this case, the modified model is as follows:

$$ \frac{\partial u}{\partial t} = d \frac{\partial^2 u}{\partial x^2} + au^2(b_1 - c_1J_1(u))(b_2 - c_2J_2(u)) - \sigma u, \quad (1) $$

where

$$ J_i(u) = \int_{-\infty}^{\infty} \phi_i(x - y)u(y, t)dy, \quad \phi_i(x) = \begin{cases} 1/(2N_i), & |x| \le N_i \\ 0, & |x| > N_i \end{cases}, \quad i = 1, 2. $$

Figure 6 shows an example of numerical simulations of equation (1). The maxima of solutions, which correspond to different species, are organized in groups that correspond to genera. Genera have a larger range of consumption of resources than species because of a greater variation of the total phenotype.

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Numerical simulations of equation (1) with double nonlocal consumption. The values of parameters are as follows: $d = 0.2$, $a = 1$, $b_1 = 2$, $b_2 = 10$, $c_1 = 1$, $c_2 = 15$, $N_1 = 3$, $N_2 = 20$, $\sigma = 0.1$. The red curve shows a snapshot of the solution and the green curve shows the positions of maxima at the $(x, t)$-plane.</figcaption></figure>

*Nonlocal prey-predator models.* Competition of species leads to a quite simple dynamic where one of the species disappears or the two species coexist. It is interesting to note that the coexistence of species is not shown in Darwin's diagram. The situation is quite different if the interaction of species is determined by prey-predator dynamics with nonlocal consumption of prey by predator [2] or with nonlocal consumption of resources by prey [3]. In these cases, numerous stationary, time periodic, or aperiodic regimes are observed.

## Acknowledgements

The authors are grateful to R. Penner for his help in the preparation of the manuscript. N. Bessonov was supported by Russian Foundation of Basic Re-
search grant 16-01-00068, 2016-2018. V. Volpert was supported by the "RUDN University Program 5-100".

## References

1. S. Atamas. Self-organization in computer simulated selective systems. Biosystems, 39 (1996), 143-151.

2. M. Banerjee, V. Volpert. Prey-predator model with a nonlocal consumption of prey. Chaos, 26, 083120 (2016).

3. M. Banerjee, V. Volpert. Spatio-temporal pattern formation in Rosenzweig-Macarthur model: effect of nonlocal interactions. Ecol. Comp., 30 (2017) 210.

4. M. Banerjee, V. Vougalter, V. Volpert. Doubly nonlocal reaction-diffusion equations and the emergence of species. Appl. Math. Model., 42 (2017), 511-599.

5. N. Bessonov, N. Reinberg, V. Volpert. Mathematics of Darwins Diagram. Math. Model. Nat. Phenom., Vol. 9, No. 3, 2014, 5-25.

6. J.A. Coyne, H.A. Orr. Speciation. Sinauer Associates, Sunderland, 2004.

7. C. Darwin. The origin of species by means of natural selection. Barnes & Noble Books, New York, 2004. Publication prepared on the basis of the first edition appeared in 1859.

8. L. Desvillettes, P.E. Jabin, S. Mischler, G. Raoul. On selection dynamics for continuous structured populations. Commun. Math. Sci., 6 (2008), no. 3, 729-747.

9. U. Dieckmann, M. Doebeli. On the origin of species by sympatric speciation. Nature, 400 (1999), 354-357.

10. M. Doebeli, U. Dieckmann. Evolutionary branching and sympatric speciation caused by different types of ecological interactions. The American Naturalist, 156 (2000), S77-S101.

11. S. Gavrilets. Fitness Landscape and the Origin of Species. Princeton University Press, Princeton, 2004.

12. S. Genieys, N. Bessonov, V. Volpert. Mathematical model of evolutionary branching. Math. Comp. Model., 2008, doi: 10/1016/j.mcm.2008.07.023

13. S. Genieys, V. Volpert, P. Auger. Pattern and waves for a model in population dynamics with nonlocal consumption of resources. Math. Model. Nat. Phenom., 1 (2006), no. 1, 65-82.

14. S. Genieys, V. Volpert, P. Auger. Adaptive dynamics: modelling Darwin's divergence principle. Comptes Rendus Biologies, 329 (2006) no. 11, 876-879.

15. S. Méléard. Random modeling of adaptive dynamics and evolutionary branching. The mathematics of Darwin's legacy, 175192, Math. Biosci. Interact., Birkhäuser/Springer Basel AG, Basel, 2011.

16. V. Volpert. Elliptic partial differential equations. Volume 2. Reaction-diffusion equations. Birkhäuser, 2014.

17. V. Volpert, N. Reinberg, M. Benmir, S. Boujena. On pulse solutions of a reaction-diffusion system in population dynamics Nonlinear Analysis, 120 (2015), 76-85.
Supplementary materials

5 Reaction-diffusion equations with local, nonlocal and global consumption of resources

We present here some properties of local and nonlocal reaction-diffusion equations and their solutions used in this work.

*Scalar reaction-diffusion equation.* Let us recall that the usual (local) scalar reaction-diffusion equation

$$ \frac{\partial u}{\partial t} = d \frac{\partial^2 u}{\partial x^2} + F(u) \qquad (1) $$

can have travelling wave solutions, that is, solutions of the form $y(x, t) = w(x - st)$, where $s$ is the wave speed. Their existence and properties depend on the function $F(u)$. Some typical examples are given by the functions:

$$ F_1(u) = au(b - cu) \quad \text{and} \quad F_2(u) = au^2(b - cu) - \sigma u. $$

The coefficients $a, b, c, d$, and $\sigma$ are positive numbers. The number of parameters in this equation can be reduced. However, it is convenient to keep them in order to analyze their influence on the behavior of solutions. For the first function $F_1(u)$ (monostable case), monotone travelling waves exist for all $s \ge s_0 = 2\sqrt{abd}$. These waves are stable in the corresponding weighted spaces. If $c\sigma/a < b^2/4$, then the second function $F_2(u)$ has three zeros: $u_+ = 0$, $u_0$ and $u_-$ are solutions of the equation $au(b - cu) = \sigma$, $u_0 < u_-$. The wave with the limits $w(\pm\infty) = u_\pm$ (bistable case) exists for a unique value of $s$. It is globally asymptotically stable. In population dynamics, equation (1) describes the distribution of the population density. The term $au^k(b - cu)$ is the reproduction rate proportional to the density $u$ of the population in the first or second power and to the available resources $(b - cu)$. The term $-\sigma u$ describes mortality of the population. Let us note that $k = 1$ corresponds to asexual reproduction or to sexual reproduction, where the reproduction rate depends only on the density of females. The value $k = 2$ corresponds to sexual reproduction where the reproduction rate depends on both densities, males and females. Assuming that these densities are equal to each other, we obtain the second power of the population density. A more general case without this assumption is studied in [22].

*Nonlocal consumption of resources.* In the case of nonlocal consumption of resources, equation (1) with the function $F_1(u)$ is replaced by the equation

$$ \frac{\partial u}{\partial t} = d \frac{\partial^2 u}{\partial x^2} + au(b - cJ(u)), \qquad (2) $$

where
$$J(u) = \int_{-\infty}^{\infty} \phi(x - y)u(y, t)dy.$$

In population dynamics, nonlocal consumption of resources means that the individuals located at the space point y consume resources not only at this point, but also in some areas around it. The consumption rate depends on the distance to resources. Therefore, the kernel of the integral depends on $(x - y)$. If we replace it by the $\delta$-function, then we will return to the local equation.

From the modelling point of view, the main interest of this equation is that the homogeneous in space stationary solution $u = b/c$ can be unstable. This is not the case for the local reaction-diffusion equation. As a result of this instability, spatial structures and periodic travelling waves can appear. In the bistable case, the corresponding equation becomes

$$\frac{\partial u}{\partial t} = d \frac{\partial^2 u}{\partial x^2} + au^2(b - cJ(u)) - \sigma u. \quad (3)$$

As in the monostable case, spatial structures and periodic travelling waves can be observed for this equation. Propagation of simple and periodic waves for the nonlocal reaction-diffusion equation was studied in [1], [4], [12], [13], [14], [9], [10], [12], [14] for $k = 1$ (monostable case) and in [1], [2], [3] for $k = 2$ (bistable case).

*Global consumption of resources.* Along with the local and nonlocal consumption of resources, we will also consider the global consumption [16] where the integral $J(u)$ in equation (3) is replaced by the integral $I(u)$:

$$\frac{\partial u}{\partial t} = d \frac{\partial^2 u}{\partial x^2} + au^2(b - cI(u)) - \sigma u, \quad I(u) = \int_{-\infty}^{\infty} u(y, t)dy. \quad (4)$$

In this case, individuals have equal access to resources independent of their location. This equation has pulse solutions. These are positive stationary solutions on the whole axis decaying at infinity.

Such solutions exist for equation (1) with the function $F_2(u)$ if $\int_0^{u_+} F_2(u)du > 0$. They can easily be found analytically. However, these solutions are unstable because the principle eigenvalue of the corresponding linearized operator is positive [17]. In the case of equation (4), existence of pulse solutions can also be easily verified. However, these solutions become stable. Their stability is observed in numerical simulations [16], but it is not yet proved analytically. If the support of the function $\phi(x)$ in equation (3) is sufficiently large, then it can also describe stable pulses or even multiple pulses where several pulses exist and slowly move from each other [21]. Their number depends on the initial condition.

Periodic travelling waves, and simple and multiple pulses represent the main types of solutions used to describe the emergence and evolution of species. Periodic travelling waves describe the process of speciation [13], [14] while pulses describe stable species without speciation [5]. Their combination and
the introduction of space dependence in the coefficients of the equations can
describe the variety of patterns in Darwin’s diagram.

References

1. N. Apreutesei, N. Bessonov, V. Volpert, V. Vougalter. Spatial structures and generalized travelling waves for an integro-differential equation. DCDS B, 13 (2010), no. 3, 537-557.

2. N. Apreutesei, A. Ducrot, V. Volpert. Travelling waves for integro-differential equations in population dynamics. Discrete Contin. Dyn. Syst., Ser. B, 11 (2009), no. 3, 541–561.

3. N. Apreutesei, A. Ducrot, V. Volpert. Competition of species with intra-specific com-
petition. Math. Model. Nat. Phenom., 3 (2008), 1–27.

4. H. Berestycki, G. Nadin, B. Perthame, L. Ryzhik. The non-local Fisher-KPP equation: travelling waves and steady states. Nonlinearity, 22 (2009), no. 12, 2813–2844.

5. N.F. Britton. Spatial structures and periodic travelling waves in an integro-differential reaction-diffusion population model. SIAM J. Appl. Math., 6 (1990), 1663–1688.

6. A. Ducrot, M. Marion, V. Volpert. Spectrum of some integro-differential operators and stability of travelling waves. Nonlinear Analysis Series A: Theory, Methods and Applications, 74 (2011), no. 13, 4455-4473.

7. G.M. Edelman. Bright air, brilliant fire: on the matter of the mind. Basic Books, 1992.

8. N. Eymard, V. Volpert, V. Vougalter. Existence of pulses for local and nonlocal reaction-diffusion equations. Submitted.

9. S.A. Gourley. Travelling front solutions of a nonlocal Fisher equation, J. Math. Biol., 41 (2000), 272–284.

10. S.A. Gourley, M.A.J. Chaplain, F.A. Davidson. Spatio-temporal pattern formation in a nonlocal reaction-diffusion equation. Dynamical systems, 16 (2001), no. 2, 173–192.

11. D. Iron, M.J. Ward. A metastable spike solution for a nonlocal reaction-diffusion model. SIAM J. Appl. Math., 60 (2000), no. 3, 778–802.

12. G. Nadin, L. Rossi, L. Ryzhik, B. Perthame. Wave-like solutions for nonlocal reaction-diffusion equations: a toy model. Math. Model. Nat.Phenom., 8 (2013), no. 3, 33–41.

13. Y. Nec, M. J. Ward. The stability and slow dynamics of two-spike patterns for a class of reaction-diffusion system. Math. Model. Nat. Phenom., 8 (2013), no.5.

14. B. Perthame, S. Genieys. Concentration in the nonlocal Fisher equation: the Hamilton-Jacobi limit. Math. Model. Nat. Phenom., 4 (2007), 135–151.

15. B.L. Segal, V.A. Volpert, A. Bayliss. Pattern formation in a model of competing popu-
lations with nonlocal interactions. Physica D Nonlinear Phenomena, 253 (2013), 12-22.

16. F. Thomas et al. Applying ecological and evolutionary theory to cancer: a long and winding road. Evol. Appl., 6 (2013), 1-10.

17. A.I. Volpert, V. Volpert, V.A. Volpert. Traveling wave solutions of parabolic systems. Translation of Mathematical Monographs, Vol. 140, AMS, Providence, 1994.

18. V. Volpert. Elliptic partial differential equations. Volume 1. Fredholm theory of elliptic problems in unbounded domains. Birkhäuser, 2011.

19. V. Volpert, S. Petrovskii. Reaction-diffusion waves in biology. Physics of Life Reviews, 6 (2009), 267–310.

20. V. Volpert, V. Vougalter. Emergence and propagation of patterns in nonlocal reaction-diffusion equations arising in the theory of speciation. In: Dispersal, individual movement and spatial ecology. M. Lewis, Ph. Maini, S. Petrovskii. Editors. Springer Applied Interdisciplinary Mathematics Series. Lecture Notes in Mathematics, Volume 2071, 2013, 331-353.

21. V. Volpert. Pulses and waves for a bistable nonlocal reaction-diffusion equation. Applied Mathematics Letters, 44 (2015), 21-25.

22. V. Volpert, N. Reinberg, M. Benmir, S. Boujena. On pulse solutions of a reactiondiffu-
sion system in population dynamics Nonlinear Analysis, 120 (2015), 76-85.

23. J. Wei, M. Winter. Existence, classification and stability analysis of multiple-peaked solutions for the Gierer-Meinhardt system in $R^1$. Methods Appl. Anal., 14 (2007), no. 2, 119-163.