# Food webs and the principle of evolutionary adaptation

Alexander S. Bratus1,2,111e-mail<EMAIL_ADDRESS>, Anastasiia V.
Korushkina${}^{3},$ Artem S. Novozhilov${}^{4},$222e-mail:
<EMAIL_ADDRESS>
${}^{\textrm{\emph{1}}}$Russian University of Transport, Moscow 127994, Russia
${}^{\textrm{\emph{2}}}$Moscow Center of Fundamental and Applied Mathematics,
Lomonosov Moscow State University, Moscow 119992, Russia
${}^{\textrm{\emph{3}}}$Faculty of Computational Mathematics and Cybernetics,
Lomonosov Moscow State University, Moscow 119992, Russia
${}^{\textrm{\emph{4}}}$Department of Mathematics, North Dakota State
University, Fargo, ND, 58108, USA

###### Abstract

A principle of evolutionary adaptation is applied to the Lotka–Volterra
models, in particular to the food webs. We present a relatively simple
computational algorithm of optimization with respect to a given criterion.
This algorithm boils down to a sequence of easy to solve linear programming
problems. As a criterion for the optimization we use the total weighted
population size of the given community and an ecological fitness, which is an
analogue of the potential energy in physics. We show by computational
experiments that it is almost always possible to substantially increase the
total weighed population size for an especially simple food web — food chain;
we also show that food chains are evolutionary unstable under the given
optimization criteria and, if allowed, evolve into more complicated structures
of food webs.

### Keywords:

Lotka–Volterra model, food chain, food web, fitness maximization

### AMS Subject Classification:

92D15, 92D25, 92D40

##  1 Controlling food chains

The replicator equations and Lotka–Volterra models of ecological communities
are among the most basic and classical mathematical models in mathematical
biology, e.g., [9]. The methods of physics and statistical physics in
particular have a long history being applied to these models, e.g. [7]. In the
vast majority of cases, however, it is assumed that the parameters in these
systems of ordinary differential equations are constant for each particular
system realization. That is, the evolution of distributions in the replicator
equation and the evolution of the population sizes in the Lotka–Volterra
equations is usually conditioned on the premise that the given set of
parameters is fixed for given circumstances. Clearly, the conditions may
change, and hence the evolution of the quantities we are interested in also
changes. However, such case is usually considered as a different possible
realization of the same system.

On the other hand one of the most basic and foundational principles in physics
is that the observed regularities are the consequence of some extremal
principle (e.g., least action principle). For many years there was a search
for a similar principle(s) to explain the observed regularities in ecological
communities (see, e.g., a historical discussion in [14]). Another possible
approach to this problem is first to formulate a mathematical model, postulate
an extremal principle (define a quantity to be maximized or minimized) and
then study the regularities which will form under this principle. In the
present text we follow this rout.

In [2] it was suggested to look at the replicator equation from an
evolutionary prospective, allowing for the changes of the parameters within
the same system realization; the basic assumption was the separation of the
time scales whereas during the fast time dynamics the system settles at an
asymptotic regime, and during the slow time dynamics the system is allowed to
change its parameter values according to some prescribed evolutionary
principle (usually the fitness maximization). This approach was further
extended and illustrated by other examples of replicator equations in [6] and
[10]. The Lotka–Volterra equations, being a close relative to the replicator
equation, are equally well suited for this approach of evolutionary
adaptation, and in the present text we show how to implement it using a
specific example of the Lotka–Volterra equations — the so-called food chains
and food webs (e.g., [5]).

Our starting point is the dynamical system that models food chain, i.e., we
assume that there is a basal species at the bottom of the food chain, there is
a predator for this basal species, which is the only prey for another
predator, and so on. Taking the intra-specific competition into account we end
up with the system of ordinary differential equations

$\dot{u}_{i}=u_{i}\Bigl{(}\rho_{i}-\bigl{(}{\boldsymbol{Au}}\bigr{)}_{i}\Bigr{)},\quad
i=1,\ldots,n,$ (1.1)

where $\rho_{1}=r_{1}>0,\rho_{j}=-r_{j},r_{j}>0,\,j=2,\ldots,n,$ are the
Malthusian parameters, in particular $\rho_{2},\ldots,\rho_{n}$ are the death
rates of the predators in the chain; matrix ${\boldsymbol{A}}$ has the form

${\boldsymbol{A}}=\begin{bmatrix}a_{11}&a_{12}&0&\ldots&0&0\\\
-a_{21}&a_{22}&a_{23}&\ldots&0&0\\\ \vdots&&&&&\vdots\\\
0&0&\ldots&0&-a_{n,n-1}&a_{nn}\\\ \end{bmatrix},$ (1.2)

where all the parameters are positive; finally, $({\boldsymbol{u}})_{i}$ is
the $i$-th element of the vector ${\boldsymbol{u}}=(u_{1},\ldots,u_{n})$.

It is well known [8, 11, 9] that if there exists an equilibrium
${\boldsymbol{\hat{u}}}$ of (1.1) in the interior of $\mathbf{R}^{n}_{+}$ then
this equilibrium is globally stable in the sense that any orbit starting in
the interior of $\mathbf{R}^{n}_{+}$ converges to ${\boldsymbol{\hat{u}}}$.

First we consider the problem to control a given food chain through the choice
of the death rates $r_{2},\ldots,r_{n}$. From the practical point of view this
is a natural thing to look into, since in this case we can imagine that we
assign quotes to harvest or hunt on particular species, which directly leads
to changes in the death rates. If we, say, increase $r_{j}$, it means that we
increase the quotes to harvest the $j$-th species; similarly, if we decrease
$r_{j}$, it implies that we employ some strategies to boost the $j$-th species
reproduction.

Within given above interpretation it is a reasonable goal to maximize the
(weighted) total population size at a fixed time moment $t=T$:

$S_{p}({\boldsymbol{u}}(t))=\sum_{i=1}^{n}p_{i}u_{i}(t)\longrightarrow\max_{t=T},$
(1.3)

where $p_{i}\geq 0$ are given weights.

To the best of our knowledge, this problem was analyzed mathematically for the
first time in [13]. In this work the functional $S_{p}({\boldsymbol{u}})$ was
maximized for the predator–prey model for the given constant $T$ and given
initial conditions, which poses some computational difficulties. In the
present text we solve the same problem using an additional hypothesis of time
separation in our model, similar to what was done for certain classes of
replicator equations in [2, 6, 10]. This hypothesis allows a significant
simplification in solving problem (1.3).

To wit, we assume that there are two time scales in our problem: fast system
dynamics when all the parameters are assumed constant, and slow dynamics, in
which we allow changes in the death rates (more general case will be
considered below). Mathematically, we have

$\frac{\,\mathrm{d}u_{i}}{\,\mathrm{d}t}(t;\epsilon t)=u_{i}(t;\epsilon
t)\Bigl{(}\rho_{i}(\epsilon t)-\bigl{(}{\boldsymbol{Au}}(t;\epsilon
t)\bigr{)}_{i}\Bigr{)},\quad i=1,\ldots,n,$ (1.4)

where $\epsilon>0$ is a small parameter. Introducing $\tau=\epsilon t$ yields

$\frac{\,\mathrm{d}u_{i}}{\,\mathrm{d}\tau}(\tau/\epsilon;\tau)=u_{i}(\tau/\epsilon;\tau)\Bigl{(}\rho_{i}(\tau)-\bigl{(}{\boldsymbol{Au}}(\tau/\epsilon;\tau)\bigr{)}_{i}\Bigr{)},\quad
i=1,\ldots,n.$ (1.5)

For the finite values $0\leq t\leq T$, due to the smallness of $\epsilon$, the
dynamics of (1.1) is close to the dynamics of (1.4). Taking the limit
$\epsilon\to 0$ implies

$\hat{u}_{i}(\tau)\Bigl{(}\rho_{i}(\tau)-\bigl{(}{\boldsymbol{A\hat{u}}}(\tau)_{i}\bigl{)}\Bigr{)}=0,\quad
i=1,\ldots,n,$ (1.6)

where $\hat{u}_{i}(\tau)=\lim_{\epsilon\to 0}u_{i}(\tau/\epsilon;\tau)$. If we
impose the restriction that $\hat{u}_{i}\neq 0$ for any $i$, the process of
evolutionary changes of the system boils down to consideration of the system
of linear equations

${\boldsymbol{A\hat{u}}}(\tau)={\boldsymbol{\rho}}(\tau),\quad{\boldsymbol{\rho}}(\tau)=\bigl{(}r_{1}(\tau),-r_{2}(\tau),\ldots,-r_{n}(\tau)\bigr{)},$
(1.7)

which depends on the slow time $\tau$. As a result, we obtain a related but
different maximization problem

$S_{p}({\boldsymbol{\hat{u}}}(\tau))=\sum_{i=1}^{n}p_{i}\hat{u}_{i}(\tau)\longrightarrow\max_{\rho(\tau)},$
(1.8)

on all possible solutions to (1.7), where vector ${\boldsymbol{\rho}}$ must
satisfy certain constraints. Specifically, we assume that parameter $r_{1}$ is
fixed, i.e., we do not allow harvesting the basal species. The rest of the
parameters must satisfy

$m_{i}\leq r_{i}(\tau)\leq M_{i},\quad i=2,\ldots,n,$ (1.9)

where $m_{i},M_{i}$ are given positive constants.

To further represent problem (1.7)–(1.9) in the form of evolutionary
adaptation, we divide this problem into a sequence of steps, at each of which
a linear programming problem should be solved. Strictly speaking this step is
not necessary for the problem at hands, but will be essential for the more
complicated problems, considered in the following.

Consider an allowable perturbation of vector ${\boldsymbol{\rho}}$ in the form

${\boldsymbol{\rho}}+\Delta{\boldsymbol{\rho}}=\bigl{(}r_{1},-(r_{2}+\Delta
r_{2}),\ldots,-(r_{n}+\Delta r_{n})\bigr{)},$

assuming

$\sum_{i=2}^{n}\Delta r_{i}\leq\varepsilon,\quad\varepsilon>0,$

where $\varepsilon$ is a sufficiently small given constant. Then the
equilibrium vector ${\boldsymbol{\hat{u}}}$ will get a perturbation
$\Delta{\boldsymbol{\hat{u}}}$. Clearly, the equation connecting
$\Delta{\boldsymbol{\rho}}$ and $\Delta{\boldsymbol{\hat{u}}}$ is the linear
system

${\boldsymbol{A}}\Delta{\boldsymbol{\hat{u}}}=\Delta{\boldsymbol{\rho}}.$
(1.10)

Consider an auxiliary system

${\boldsymbol{A}}^{\top}{\boldsymbol{v}}={\boldsymbol{p}}.$ (1.11)

After taking the scalar product of both sides of (1.10) with
${\boldsymbol{v}}$ we obtain

$\left\langle\Delta{\boldsymbol{\rho}}\,,{\boldsymbol{v}}\right\rangle=\left\langle{\boldsymbol{A}}\Delta{\boldsymbol{\hat{u}}}\,,{\boldsymbol{v}}\right\rangle=\left\langle\Delta{\boldsymbol{\hat{u}}}\,,{\boldsymbol{A}}^{\top}{\boldsymbol{v}}\right\rangle=\left\langle\Delta{\boldsymbol{\hat{u}}}\,,{\boldsymbol{p}}\right\rangle,$

from where

$\Delta
S_{p}=\left\langle\Delta{\boldsymbol{\hat{u}}}\,,{\boldsymbol{p}}\right\rangle=-\sum_{i=2}^{n}v_{i}\Delta
r_{i}.$

Let $K^{+}$ ($K^{-}$) be the subsets of the set of indexes $\\{2,\ldots,n\\}$,
for which the components of ${\boldsymbol{v}}$ that solves (1.11) are positive
(negative) and, assuming we have no zero components, $K^{+}\bigcup
K^{-}=\\{2,\ldots,n\\}$.

The form of the expression $\Delta S_{p}$ immediately implies that if both of
the sets $K^{+}$ and $K^{-}$ are nonempty then there is a potential way to
increase the value of $S_{p}$ by either choosing $\Delta r_{i}<0$ for $i\in
K^{+}$ or $\Delta r_{i}>0$ for $i\in K^{-}$. Moreover, if $K^{+}=\emptyset$
($K^{-}=\emptyset$) then increase of the total population size at the expense
of increasing harvest of some specific species (decreasing it) is impossible.
We say “a potential way to increase the value of $S_{p}$” because the analysis
of perturbations $\Delta r_{i}$ and $\Delta{\hat{u}}_{i}$ should also include
the analysis of given constraints. In particular, recall that $M_{i}\geq
r_{i}\geq m_{i}$ for all $i$. So, if we find that, for instance, $r_{i}\leq
m_{i}$ then one must add that $\Delta r_{i}\geq 0$ with certain changes for
the expression for $\Delta S_{p}$. Similarly, since our hypothesis of time
separation is based on the fact that the internal equilibrium of the food
change is globally stable, if it exists, we also require that all
$\hat{u}_{i}\geq\delta>0$, and if for some index $i$ we already approached the
coordinate plane, i.e., $\hat{u}_{i}=\delta$ then one must additionally
require that

$\Delta\hat{u}_{i}=\Bigl{(}{\boldsymbol{A}}^{-1}\Delta{\boldsymbol{\rho}}\Bigr{)}_{i}\geq
0.$

Taking everything together, we obtain an efficient numerical procedure to
solve problem (1.8), and hence, with high degree of accuracy, approximate
significantly more computationally involved problem (1.3).

Here is the basic algorithm. We start with fixing the values of
$m_{i},M_{i},i=2,\ldots,n,\varepsilon>0,\delta>0$. Let
${\boldsymbol{\rho}}^{0}$ satisfies the given constraints. We first solve the
auxiliary problem (1.11), find ${\boldsymbol{v}}$, and consider the linear
programming problem

$\Delta S_{p}=-\sum_{i=2}^{n}v_{i}\Delta r_{i}^{0},\quad\sum_{i=2}^{n}\Delta
r_{i}^{0}\leq\varepsilon,\quad i=2,\ldots,n,$ (1.12)

and such that the vectors
${\boldsymbol{\rho}}^{0}+\Delta{\boldsymbol{\rho}}^{0}$ and
${\boldsymbol{\hat{u}}}$ are allowable. It means that
${\boldsymbol{\rho}}^{0}+\Delta{\boldsymbol{\rho}}^{0}$ within the given
constants $m_{i},M_{i}$, and
${\boldsymbol{\hat{u}}}^{0}+\Delta{\boldsymbol{\hat{u}}}^{0}$ is away from the
boundary of $\mathbf{R}^{n}_{+}$ at least $\delta$ units. If we are capable to
find such $\Delta{\boldsymbol{\rho}}^{0}$ we compute

${\boldsymbol{\rho}}^{1}={\boldsymbol{\rho}}^{0}+\Delta{\boldsymbol{\rho}}^{0}$

and repeat the previous step, until we can make an allowable change. This
problem is convex because of the convexity of the set of imposed constraints
and linearity of problem (1.7), and therefore its solution always exists,
although it could be non unique.

Here is an example. We note that we use this example only for illustrative
purposes to show how our algorithm of evolutionary adaptation works, without
pretending that a real food chain is considered. We consider the case of
$n=10$ species in the food chain.

Let

${\boldsymbol{A}}=\left[\begin{array}[]{cccccccccc}1&4&0&0&0&0&0&0&0&0\\\
-4&1&6&0&0&0&0&0&0&0\\\ 0&-7&1&5&0&0&0&0&0&0\\\ 0&0&-6&1&3&0&0&0&0&0\\\
0&0&0&-8&1&6&0&0&0&0\\\ 0&0&0&0&-7&1&9&0&0&0\\\ 0&0&0&0&0&-2&1&3&0&0\\\
0&0&0&0&0&0&-4&1&5&0\\\ 0&0&0&0&0&0&0&-2&1&9\\\ 0&0&0&0&0&0&0&0&-3&1\\\
\end{array}\right]$ (1.13)

and

${\boldsymbol{\rho}}=(16,-1,-1,-1,-1,-1,-1,-1,-1,-2)$ (1.14)

in the problem $\eqref{eq1:1}$. We choose the following weights

${\boldsymbol{p}}=(10,20,30,40,50,60,10,10,10,10).$ (1.15)

For the given parameter values we find, keeping one digit after the decimal
point, that

${\boldsymbol{\hat{u}}}=(4.2,2.9,2.2,3.5,2.8,4.0,1.6,1.8,0.7,0.2)\in\mathbf{R}^{n}_{+}$

solves (1.10), and, solving (1.11),

${\boldsymbol{v}}=(48.1,9.5,26.0,8.9,12.4,-1.6,6.3,-4.5,2.2,-10.0).$

Therefore it looks like there is a potential to improve the value of $S_{p}$,
which is initially $S_{p}(0)=731$.

The constraints for the problem are

$0.1\leq r_{i}\leq 3,\quad i=2,\ldots,n,$ $\hat{u}_{i}>\delta=0.001,\quad
i=1,\ldots,n,$

and, finally, $\varepsilon=0.01.$

Now we can run a sequence of steps, at each of which a linear programming
problem is solved, until our solution moves outside of the allowable set. For
the given set of parameters the algorithm runs for approximately 550 steps, at
which the value of $S_{p}$ is 805, see Fig. 1.1.

Figure 1.1: The change in the value of functional $S_{p}$ given by (1.8) at
each step of evolutionary algorithm in Section 1.

As we see in the coordinates of vector ${\boldsymbol{v}}$ for the given
example, the algorithm tries to increase the sixth, eighth, and tenth death
rates and decrease all others, see Fig. 1.2.

Figure 1.2: The change in the death rates ${\boldsymbol{r}}$ at each step of
evolutionary algorithm in Section 1.

In general what this simple numerical experiment shows is that if one is
capable to modulate the death rates of a given food chain within a reasonable
interval of allowable values, it is almost always possible (and was always
possible in our numerical experiments) to increase the weighted total
population size at the equilibrium. Moreover, an extremely simple numerical
procedure of finding the auxiliary vector ${\boldsymbol{v}}$ in (1.11)
indicated which coefficients of the vector ${\boldsymbol{\rho}}$ in (1.1)
should be increased and which decreased during the evolutionary steps.

As a final remark here, we note that from a practical point of view it it the
harvesting of the top predator in a given food chain that of the most
interest. This specific situation can be easily modeled with the described
approach by setting vector ${\boldsymbol{p}}$ to $(0,\ldots,0,1)$.

##  2 From food chains to food webs

In the previous section we analyzed only food chains, i.e., communities, for
which the interaction graph of matrix ${\boldsymbol{A}}$ has the simplest
linear structure, where any node, except for the first and the last ones, has
only two neighbors, one is a prey and another one is a predator. In reality
these graphs are certainly more complicated and represent what is generally
called wood web. Let us ask therefore the following question: Is it
advantageous for a given food chain to modify itself (i.e., to add new trophic
interactions) under given circumstances (evolutionary constraints)?

Mathematically it means that now we allow evolutionary changes in the elements
of matrix ${\boldsymbol{A}}$. For simplicity, we first assume that the
predator with number $k$ has the ability to feed on any species with numbers
$k-j,j=2,\ldots,k-1$, where the basal species has the number $1$. For this
computation we fix vector ${\boldsymbol{\rho}}$, and the elements of matrix
${\boldsymbol{A}}$ are now chosen from the allowable set
$\sum_{i,j=1}^{n}a_{ij}^{2}\leq K$ for some prescribed constant $K>0$. In this
case the elements of the perturbation
${\boldsymbol{A}}+\Delta{\boldsymbol{A}}$ of the original matrix $\Delta A$
must be negative, if they are under the main diagonal, and positive if they
are above the main diagonal. As before, to use the principle of the time scale
separation, we assume that each perturbation $\Delta a_{ij}$ satisfies

$|\Delta a_{ij}|\leq\varepsilon$ (2.1)

for some given constant $\varepsilon$.

We also call any perturbation $\Delta{\boldsymbol{A}}$ from the set (2.1)
allowable, if the spherical norm of matrix ${\boldsymbol{A}}$ is non-
increasing, which, up to the second order, means that

$\sum_{i,j=1}^{n}a_{ij}\Delta a_{ij}\leq 0.$ (2.2)

Assuming that all the perturbations are proportional to $\varepsilon>0$ we
have, up to the order $O(\epsilon^{2})$,

${\boldsymbol{A}}\Delta{\boldsymbol{\hat{u}}}+\Delta{\boldsymbol{A}}{\boldsymbol{\hat{u}}}=0.$
(2.3)

Multiplying the equality (2.3) by ${\boldsymbol{v}}$ that solves (1.11), we
find that

$\Delta
S_{p}=\left\langle\Delta{\boldsymbol{\hat{u}}}\,,{\boldsymbol{p}}\right\rangle=-\left\langle\Delta{\boldsymbol{A}}{\boldsymbol{\hat{u}}}\,,{\boldsymbol{v}}\right\rangle,$
(2.4)

which implies that the total weighted population size can be potentially
increased if it is possible to find the allowable changes of
$\Delta{\boldsymbol{A}}$, which make the expression (2.4) positive. Hence, at
each evolutionary step our goal is to maximize (2.4). We note that in this
case we are not guaranteed to reach the global maximum.

Before we illustrate the process of evolutionary adaptation in this case, we
should pause and note two complications that arise here. These complications
are also relevant to the application of the described methods to general
Lotka–Volterra models.

In the previous section we were able to use the process of evolutionary
adaptation because the food chain is globally stable if there exists an
internal unique equilibrium. If we allow, as in this section, for the matrix
${\boldsymbol{A}}$ stops being of the form (1.2), and assume that $|a_{ij}|>0$
if $|i-j|>2$, then the presence of the internal equilibrium does not guarantee
that the system during the fast time scale settles at it; it is actually
possible to have much more complicated scenarios in this case, see [9]. To
make sure that the process of evolutionary adaptation works as desired we must
require one more condition at each step: that the system that we obtain after
our perturbation is permanent.

Recall that general Lotka–Volterra system (1.1) is permanent, if its orbits,
after a sufficiently large initial time, end up in a compact set
$L\subseteq\operatorname{int}\mathbf{R}^{n}_{+}$, where
$\operatorname{int}\mathbf{R}^{n}_{+}$ is the interior of
$\mathbf{R}^{n}_{+}$. A necessary condition for (1.1) to be permanent is the
existence of unique internal equilibrium ${\boldsymbol{\hat{u}}}\in L$. From
the discussion above it follows that for the food chains this condition is
also sufficient, but it is no longer true for general ${\boldsymbol{A}}$. What
is most important for us is the fact that if system (1.1) is permanent then

$\lim_{T\to\infty}{\boldsymbol{\bar{u}}}=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}{\boldsymbol{u}}(t)\,\mathrm{d}t={\boldsymbol{\hat{u}}},$
(2.5)

i.e., the averages ${\boldsymbol{\bar{u}}}$ along the orbits are exactly the
coordinates of the internal equilibrium in the limit $T\to\infty$. This means
that for the permanent systems the coordinates of the unique equilibrium still
contain significant information about system behavior even if the system
itself is not attracted to ${\boldsymbol{\hat{u}}}$.

Therefore, for our evolutionary algorithm to work we also at each step should
check whether the obtained allowable perturbation implies a permanent food
web. There are efficient, but computationally involved, algorithms to perform
this task (e.g., [12]), but for simplicity we shall resort to only checking
the condition (2.5), as our numerical experiments show this works well for
most of the time. Specifically, we run our system from a number of random
initial conditions and calculate time averages along the orbits; if they are
sufficiently close to the coordinates of the equilibrium
${\boldsymbol{\hat{u}}}$ we assume that the system is permanent and we can
continue out adaptation.

The second complication is that we consider such changes in the matrix
${\boldsymbol{A}}$ that non necessarily lead to the true predator–prey
relation. Namely, we allow, for instance, $a_{25}$ to be positive and $a_{52}$
be zero, or $a_{52}$ be negative and $a_{25}$ be zero, hence including more
general trophic interactions.

Now we are ready to summarize the algorithm for the evolutionary adaptation if
matrix ${\boldsymbol{A}}$ is allowed to be changed. As usual, we start with an
allowable matrix ${\boldsymbol{A}}$, and maximize $\Delta S_{p}$ in (2.4)
assuming that $\Delta a_{ij}$ satisfy (2.1) and (2.2). Moreover, we also
assume additionally that elements of ${\boldsymbol{A}}+\Delta{\boldsymbol{A}}$
are such that all the entries above the main diagonal are non-positive, and
all the entries below the main diagonal are non-negative. Finally, we check
that the new coordinates ${\boldsymbol{\hat{u}}}+\Delta{\boldsymbol{\hat{u}}}$
are at distance at least $\delta$ units from the boundary of
$\mathbf{R}^{n}_{+}$ and that the resulting system is permanent. If all the
conditions are met, we can make the next evolutionary step in maximizing the
total weighted population size.

To see what may happen in this scenario, we take ${\boldsymbol{A}}$,
${\boldsymbol{\rho}}$, ${\boldsymbol{p}}$ as in (1.13), (1.14), and (1.15). As
before we take $\delta=0.001$ and $\varepsilon=0.01$. Following the steps
outlined above, we find that our evolutionary algorithm is capable to increase
$S_{p}$ approximately 1.8 times, the final value is $S_{p}(end)=1300$ (recall
that $S_{p}(0)=731$), Fig. 2.1.

Figure 2.1: The change in the value of functional $S_{p}$ given by (1.8) at
each step of evolutionary algorithm in Section 2.

This time however, the algorithm stops because one of the equilibrium
coordinates approaches zero, see Fig. 2.2.

Figure 2.2: The changes in the equilibrium coordinates in the example of
evolutionary algorithm in Section 2.

For comparison we also show the averages along the orbits in this particular
case, Fig. 2.3.

Figure 2.3: The averages along the orbits at each evolutionary step of the
evolutionary algorithm in Section 2. Compare with Fig. 2.2.

Finally, the resulting graph is shown in Fig. 2.4, where the numbering of the
vertexes corresponds to the linear order of the species in the original food
chain, number 12 being the basal species, and number 10 being the top
predator. We have both arrows connecting $i$ and $j$ is they are in the
relation “predator–prey” and only one arrow if the relation of amensalism or
commensalism.

Figure 2.4: The trophic graph of the resulting food web in the numerical
example in Section 2.

##  3 Maximizing the mean population fitness in the Lotka–Volterra systems

In the previous two sections we considered applications of the principle of
evolutionary adaptation to the food webs using as an optimization criterion
the total (weighted) population size. Such an approach is very natural from
the point of view of artificial management of the resources because it is the
total population sizes of some of the species in a given food web that are
frequently the most important quantities from the point of view of harvesting
or survival of some species. From the point of view of natural evolution this
criterion is doubtful at best. Hence we would like to consider a different
extreme principle for the possible evolution of food webs.

As in many other sciences the search for the extreme principle that shapes the
population communities has a long history, see some discussion in [14]; for
our purposes we define the population fitness in the Lotka–Volterra models
(1.1) as

$F({\boldsymbol{u}})=\sum_{i=1}^{n}\rho_{i}u_{i}-\frac{1}{2}\left\langle{\boldsymbol{Au}}\,,{\boldsymbol{u}}\right\rangle,$
(3.1)

for the given population vector ${\boldsymbol{u}}$. To the best of our
knowledge this function was used to find the correlations of system’s dynamics
with some extremal principle (namely, maximizing (3.1)) for the first time in
[14]. As some justification we note that if ${\boldsymbol{A}}$ is symmetric,
then $F$ is increasing along the orbits of (1.1). Indeed,

$\displaystyle\dot{F}({\boldsymbol{u}})$
$\displaystyle=\sum_{i=1}^{n}\rho_{i}^{2}u_{i}-2\sum_{i=1}^{n}({\boldsymbol{Au}})_{i}\rho_{i}u_{i}+\sum_{i=1}^{n}({\boldsymbol{Au}})_{i}^{2}=$
$\displaystyle=\sum_{i=1}^{n}(\rho_{i}-({\boldsymbol{Au}})_{i})^{2}u_{i}\geq
0.$

This property is no longer true if ${\boldsymbol{A}}$ is not symmetric, but we
still consider (3.1) as the optimality criterion that we expect to be
maximized during the long time evolution of ecological community. Certainly,
it is well known that for many mathematical models of evolution the mean
population fitness is actually not increasing and not monotone (see, e.g., [1,
3, 4] for a detailed discussion), but as a first approximation of the action
of natural selection on the ecological communities it is arguably the most
logical extremal principle.

To summarize, in what follows we apply the principle of evolutionary
adaptation to the permanent Lotka–Volterra model (1.1) under the goal to
maximize the functional (3.1) at the internal equilibrium
${\boldsymbol{\hat{u}}}$, allowing now for both the vector
${\boldsymbol{\rho}}$ and matrix ${\boldsymbol{A}}$ to be changed during the
slow time evolution. The only requirement, as before, is the constraint to
keep them in the form to describe food webs.

As before we consider the linear system (1.7), where now we assume that both
the elements of ${\boldsymbol{A}}$ and ${\boldsymbol{\rho}}$ depend smoothly
on the slow time $\tau$, and also satisfy the constraints

$-m_{i}\geq\rho_{i}\geq-M_{i},\quad
i=2,\ldots,n,\quad\sum_{i,j=1}^{n}a_{ij}^{2}\leq K,$ (3.2)

for given positive constants $m_{i},M_{i},K$.

Denote
$\Delta{\boldsymbol{A}},\Delta{\boldsymbol{\hat{u}}},\Delta{\boldsymbol{\rho}}$
the main parts of perturbations at each evolutionary step of
${\boldsymbol{A}},{\boldsymbol{\hat{u},{\boldsymbol{\rho}}}}$ respectively (we
omit the dependence on $\tau$). Our first goal is to find the perturbation the
fitness functional (3.1) gets in this case.

From (1.7) it follows that

$\Delta{\boldsymbol{A}}{\boldsymbol{\hat{u}}}+{\boldsymbol{A}}\Delta{\boldsymbol{\hat{u}}}=\Delta{\boldsymbol{\rho}}.$
(3.3)

Now,

$\Delta
F=\left\langle{\boldsymbol{\rho}}\,,\Delta{\boldsymbol{\hat{u}}}\right\rangle+\left\langle\Delta{\boldsymbol{\rho}}\,,{\boldsymbol{\hat{u}}}\right\rangle-\frac{1}{2}\left\langle({\boldsymbol{A}}+{\boldsymbol{A}}^{\top}){\boldsymbol{\hat{u}}}\,,\Delta{\boldsymbol{\hat{u}}}\right\rangle-\frac{1}{2}\left\langle\Delta{\boldsymbol{A}}{\boldsymbol{\hat{u}}}\,,{\boldsymbol{\hat{u}}}\right\rangle.$

From (3.3) we have

$\Delta{\boldsymbol{\hat{u}}}={\boldsymbol{A}}^{-1}(\Delta{\boldsymbol{\rho}}-\Delta{\boldsymbol{A}}{\boldsymbol{\hat{u}}}),$

therefore,

$\Delta
F=\frac{1}{2}\left\langle{\boldsymbol{\rho}}\,,\Delta{\boldsymbol{\hat{u}}}\right\rangle+\frac{1}{2}\left\langle\Delta{\boldsymbol{\rho}}\,,{\boldsymbol{\hat{u}}}\right\rangle=\frac{1}{2}\left(\left\langle({\boldsymbol{A}}^{-1})^{\top}{\boldsymbol{\rho}}\,,\Delta{\boldsymbol{\rho}}-\Delta{\boldsymbol{A}}{\boldsymbol{\hat{u}}}\right\rangle+\left\langle{\boldsymbol{\hat{u}}}\,,\Delta{\boldsymbol{\rho}}\right\rangle\right).$

If we introduce an auxiliary linear system

${\boldsymbol{A}}^{\top}{\boldsymbol{v}}={\boldsymbol{\rho}}\implies{\boldsymbol{v}}=({\boldsymbol{A}}^{\top})^{-1}{\boldsymbol{\rho}},$
(3.4)

we finally find

$\Delta
F=\frac{1}{2}\left(\left\langle{\boldsymbol{\hat{u}}}+{\boldsymbol{v}}\,,\Delta{\boldsymbol{\rho}}\right\rangle-\left\langle{\boldsymbol{v}}\,,\Delta{\boldsymbol{A}}{\boldsymbol{\hat{u}}}\right\rangle\right),$
(3.5)

where the perturbations $\Delta{\boldsymbol{\rho}}$ and
$\Delta{\boldsymbol{A}}$ must be allowable, i.e., satisfy the conditions
(3.2), keep the structure of the system as for a food web, maintain the
uniqueness of the internal equilibrium
${\boldsymbol{\hat{u}}}+\Delta{\boldsymbol{\hat{u}}}$, and guarantee that the
resulting system is permanent, similar to what we discussed in more details in
the previous section.

To illustrate the algorithm consider the same initial matrix
${\boldsymbol{A}}$ and initial vector ${\boldsymbol{\rho}}$ as in the example
of Section 1, i.e., (1.13) and (1.14). Solving the linear programming problem
on each step for maximizing $\Delta F$, we find that it is possible to
increase the mean population fitness almost 4.4 times, see Fig. 3.1.

Figure 3.1: Increase of the mean population fitness in the numerical example
of evolutionary adaptation to maximize (3.1).

For this specific calculation it turns out that at some step the system stops
being permanent, which we check in the same way as was described in Section 2.
The changes in the coordinates of the equilibrium ${\boldsymbol{\hat{u}}}$ are
shown in Fig. 3.2.

Figure 3.2: The changes in the equilibrium coordinates in the numerical
example of evolutionary adaptation to maximize (3.1).

As before, we find that the major changes come as a result of turning the
initial food chain into a much more complicated structure of a food web, see
Fig. 3.3. It is an interesting open question to find some regularities in the
final trophic web, obtained as a result of evolution under the principle of
maximizing the mean population fitness (3.1).

Figure 3.3: The trophic graph of the resulting food web in the example of
Section 3. See details in the text.

##  4 Concluding words

Here we presented a relatively simple computational algorithm to incorporate
an explicit time dependent community evolution in the classical Lotka–Volterra
models. This is an approximate algorithm, however it consists in a sequence of
steps at each of which a linear programming problem is solved, which can be
done very fast and accurately. We illustrate this algorithm by using two
different optimization criteria, namely, optimizing the total weighted
population size and what we call ecological fitness, which can be considered
as an analogue of the potential energy in physics. We show that our algorithm
performs well being able to increase the required functional in a significant
way.

First we apply our algorithm to the food chains and show that it is always
possible, by managing the death rates of the system, to increase the total
population size of the community. Moreover, a simple algebraic test indicated
which death rates should be increased, and which — decreased.

We also consider more complicated cases allowing for the appearance and
disappearance of new trophic interactions. An interesting observation here is
that the food chain, which is our usual initial step, is not evolutionary
stable under given extremal principals, and it always evolves into a more
complicated structure of food web.

It is an interesting open problem to identify the regularities in the obtained
under evolutionary adaptations systems, which is part of the ongoing work.

### Acknowledgements:

ASB is supported by the Russian Science Foundation Grant 19-11-00008 and by
the Ministry of Science and Higher Education Grant 075-15-2019-1621.

## References

  * [1] J. Birch.  Natural selection and the maximization of fitness.  Biological Reviews, 91(3):712–727, 2016.
  * [2] A. S. Bratus, S. Drozhzhin, and T. Yakushkina.  On the evolution of hypercycles.  Mathematical biosciences, 306:119–125, 2018.
  * [3] A. S. Bratus, A. S. Novozhilov, and Y. S. Semenov.  Adaptive fitness landscape for replicator systems: To maximize or not to maximize.  Mathematical Modelling of Natural Phenomena, 18(3):25, 2018.
  * [4] A. S. Bratus, A. S. Novozhilov, and T. Yakushkina.  Geometry of Fitness Surfaces and Dynamics of Replicator Systems, pages 69–77.  Springer International Publishing, Cham, 2021.
  * [5] J. E. Cohen, F. Briand, and Ch. M. Newman.  Community food webs: data and theory.  Springer-Verlag, 2012.
  * [6] S. Drozhzhin, T. Yakushkina, and A. S. Bratus.  Fitness optimization and evolution of permanent replicator systems.  Journal of Mathematical Biology, 82(3):1–26, 2021.
  * [7] N. S. Goel, S. C. Maitra, and E. W. Montroll.  On the volterra and other nonlinear models of interacting populations.  Reviews of modern physics, 43(2):231, 1971.
  * [8] G. W. Harrison.  Global stability of food chains.  The American Naturalist, 114(3):455–457, 1979.
  * [9] J. Hofbauer and K. Sigmund.  Evolutionary Games and Population Dynamics.  Cambridge University Press, 1998.
  * [10] I. Samokhin, T. Yakushkina, and A. S. Bratus.  Open quasispecies systems: New approach to evolutionary adaptation.  arXiv preprint arXiv:2011.11742, 2020.
  * [11] J. W. H. So.  A note on the global stability and bifurcation phenomenon of a lotka-volterra food chain.  Journal of theoretical biology, 80(2):185–187, 1979.
  * [12] P. F. Stadler and R. Happel.  The probability of permanence.  Mathematical biosciences, 113(1):25–50, 1993.
  * [13] Yu. M. Svirezhev and E. Ya. Yelizarov.  Mathematical modelling of biological systems.  Nauka, Moscow (in Russian), 1972.
  * [14] Yu.M. Svirezhev and D. O. Logofet.  The stability of biological communities.  Nauka, 1978.