## Abstract

Phenotypically structured equations arise in population biology to describe the interaction of species with their environment that brings the nutrients. This interaction usually leads to the selection of the fittest individuals. Models used in this area are highly nonlinear, and the question of long-term behaviour is usually not solved. However, there is a particular class of models for which convergence to an evolutionary stable distribution is proved, namely when the quasi-static assumption is made. This means that the environment, and thus the nutrient supply, reacts immediately to the population dynamics. One possible proof is based on a total variation bound for the appropriate quantity. We extend this proof to cases where the nutrient is regenerated gradually. A simple example is the chemostat with a rendering factor, then our result does not use any smallness assumption. For a more general setting, we can treat the case with a fast equilibration of the nutrient concentration.

## 1. Introduction

In population biology, long-term behaviour for phenotypically structured models is a difficult question related to interaction with environmental conditions, selection of fittest trait and lack of dissipation principles. The *competitive exclusion principle* is a famous general result, and states that, with a single type of ‘niche’ or substrate, a single trait is selected.

A typical example, where this can be proved rigorously, is the chemostat model
*n*(*t*,*x*) of individuals which at time *t* have the trait *x*. The substrate, whose concentration is denoted by *S*, is used with a trait-dependent uptake coefficient *η*(*x*,*S*) and a rendering factor *a*(*x*). The renewal of the reactor, with fresh nutrient *S*_{0}, occurs with the rate *R*_{0}.

The simplest situation is when there is a unique evolutionary stable distribution (ESD; a term coined in [1]) which concentrates in a single Dirac mass. That means there is a unique trait *η* is increasing with *S*. And, the second equality gives

Then, it is known when *a*≡1, see [2], that the competitive exclusion principle can be expressed as

However, we do not know general assumptions on *η*^{n}, *η*^{S} which would lead to a similar result for the more general chemostat model
*η*^{n}(*x*,*S*(*t*)) and *η*^{S}(*x*,*S*(*t*)) depends not only on the trait *x*, but also on *S*. A general method is to use a Lyapunov functional (entropy), but this requires a particular structure on the system [1,3,4].

The laws for nutrient delivery and consumption may differ for other models [5,6], but similar questions still arise. A mathematical model, which contains (1.1)–(1.2) with *η*(*x*,*S*)=*η*(*S*) as a particular case, can be written as
*η*(*x*,*S*)=*η*_{1}(*x*)*η*_{2}(*S*), can be handled. Here, *R*(*x*,*S*) denotes a trait-dependent birth–death rate, *S* is still the nutrient concentration, and *ρ*(*t*) is a measure of the pressure exerted by the total population for nutrient consumption with rate *Q*. The parameter *β*, which obviously could be included in *Q* is used here for a simple mathematical purpose. It gives a time scale which, in the limit *β*=0, just gives 0=*Q*(*S*(*t*),*ρ*(*t*)). Under suitable assumptions, this equation can be inverted in *S*=*q*(*ρ*). In this case, the long-term selection of the ESD, (1.4), is known to hold [7–9].

Our aim is to prove the same convergence result to an ESD, (1.4), when *β* is small. Section 3 is devoted to prove the result and a precise statement is given in theorem 3.1. In order to make the proof more intuitive, we begin with the simpler case of the chemostat system (1.1)–(1.2); this is developed in §2. As notation, we define for a function *g* depending on the variable *y* the partial derivative *g*_{y}:=(∂/∂*y*)*g*.

## 2. The chemostat with rendering factor

The model of the chemostat with a rendering factor is defined by the system (1.1)–(1.2). We complete it with initial data *n*^{0}(*x*), *S*^{0} that satisfy

In order to analyse the long-term behaviour, we need assumptions on the problem parameters and coefficients. Namely we need to ensure first non-extinction which follows from the assumptions
*a*(*x*)≤1, but here we use only that for some constants *a*_{m}, *a*_{M}>0

Then, we have the following generalization of the case *a*≡1 which is treated in [2].

### Theorem 2.1

*With assumptions (2.1)–(2.4), there are constants ρ*_{m}*, ρ*_{M}*, such that
**Assuming also (1.3), as* *we obtain
*

### Proof.

*First step. A conserved quantity.* For future use, we define
*a*, integrating and adding equation (1.2), we obtain

The *a priori* bounds follow easily. Because *n*>0 (from (2.1)), we find *S*≤*S*_{0} and because *η*(*x*,0)=0 from assumption (2.2), we find *S*>0. For an upper bound on *ρ*, we use that *u*(*t*) is bounded and assumption (2.4). We find
*ρ*_{m} can be derived in the same way, using *a*_{m}.

*Second step. Bounded variation (BV) estimates of* *J* in (2.5), we have, using (2.6),

We define the negative part of *J* by *ρ*, it follows that *J* has bounded variation and *u*(*t*) converges to 0, we conclude that *S*(*t*) has a limit
*Third step. The limits.* At this stage, we can identify *ρ*, we immediately conclude that the growth rate should vanish on the long term, that is written *S* of *η*, this tells us that *n*(*t*,*x*) concentrates as a Dirac mass at the point *S*(*t*) and *u*(*t*), we know that

Theorem 2.1 is proved. ▪

## 3. The general setting

In the general setting of the system (1.5)–(1.7), the same proof does not apply *per se*. This is because we do not dispose of a quantity, as *u*(*t*) in the previous proof, which is easy to control and brings us back to the quasi-steady state where *S* is a function of *ρ*. For this reason, we need a smallness condition which is well expressed in terms of *β*. With this condition, we can build a quantity which belongs to *J*(*t*) in the previous proof.

### (a) Assumptions and main result

We complete the system (1.5)–(1.7) with initial conditions *S*^{0}, *n*^{0}, which are compatible with some invariant region of interest
*S*_{m} and *S*_{0} in the assumptions below, this assumption for *S*^{0} is made to simplify the statements and can be seen as a generalization of those for the chemostat in §2).

Next, for the Lipschitz continuous functions *R* and *Q*, we assume that there are constants *S*_{0}>0, *K*_{Q}>0… such that

Note that from assumption (3.2), we directly obtain the bounds

With these assumptions, the smallness condition on *β* can be written as
*ρ*_{M}, *S*_{m} below, which depends only on the assumptions above).

### Theorem 3.1

*With assumptions (3.1)–(3.4), there is a constant ρ*_{M} *such that
**where the value S*_{m}*<S*_{0} *is defined by Q*(*S*_{m}*,ρ*_{M})=0 *and this exists, thanks to assumption* (*3.2*).

*Assuming also (3.6), ρ*^{2}*(t) has bounded total variation. Consequently, there are limits* *and
*

As for the chemostat, the solution can go extinct, that means *n*(*t*) concentrates on the maximum points of

The end of this section is devoted to prove theorem 3.1. This requires to adapt the method introduced in [2,7,15] which is to prove that *ρ*(*t*) has a bounded total variation. This method works well in the quasi-static case, that is *β*=0. The adaptation is not as direct as one could think in view of §2.

### (b) An upper bound for *ρ*

This step is not as simple as usual. Integrating equation (1.5) with respect to *x* and using assumption (3.2), yields that
*Q*, we have *Q*(*S*,*ρ*)≤−*K*_{Q}*ρ*+*Q*(0,0), adding equation (1.6), we obtain the inequality
*C*_{2} the root in *ρ*_{M} for *ρ*(*t*).

From this upper bound, we obtain the lower bound on *S*(*t*) because

### (c) Bounded variation estimates

Our next goal is to find a quantity which converges as

*First step. Equations on * With these definitions, from equations (1.5) and (1.6), we can write

*t*and

*x*sometimes. With the definitions

*n*and

*S*, we obtain

*Second step. Bound on a linear combination of P and J.* Now, we consider a linear combination of

*P*and

*J*, where

*μ*(

*t*) is a function to be determined later. We write

*μ*(

*t*) such that the expression between the second pair of parentheses in equation (3.12) is zero. In other words,

*μ*(

*t*) solves the differential equation

*μ*(

*t*)>0 of (3.13) for all times. To do so, we note that the zeroes of −

*β*|

*Q*

_{ρ}|

*μ*

^{2}+

*μ*|

*Q*

_{S}|−

*α*are

We are going to find two constants 0<*μ*_{m}<*μ*_{M} such that, choosing initially *μ*_{m}<*μ*(0)<*μ*_{M}, then we have for all times
*μ*_{m} defined by the condition

We first show how to enforce the inequality (3.15). We use that, for 0≤*x*≤1, the concavity inequality holds, *β* (3.6) is enough to obtain the inequality (3.15).

The lower bound in (3.14) is because the condition (3.15) imposes *μ*_{m}∈(*μ*_{−}(*t*),*μ*_{+}(*t*)) and thus *t*≥0.

For the upper bound in (3.14), we obtain
*μ*(*t*) and coming back to equation (3.12), we arrive at

*Third step. L*^{1}-*bound on P.* From the above inequality, we wish to prove that

*P*(

*t*) is integrable on the half-line. Adding

*α*(

*P*/

*βμ*) to (3.10), we find the ODE

*P*is bounded combined with the estimate (3.16), we have for some constant

*C*

With the lower bound on *α* in (3.9), we conclude that
*ρ*(*t*) is bounded, we finally find that (*d*/*dt*)*ρ*^{2} is bounded in *ρ*^{2} has a limit for *ρ* has a limit

*Fourth step. Conclusion.* Because *ρ*(*t*) has a long-term limit *Q*, more precisely *Q*_{S}<0 in (3.2), shows *S*(*t*) also has a limit

As usual, [7,10,11], we can conclude that *x*. Otherwise, *n*(*t*,*x*) would have an exponential growth for *x* in an open set, which would imply exponential growth for large times, a contradiction with the upper bound on *ρ*(*t*).

This gives the statements for theorem 3.1.

### (d) Numerical considerations

For *β* large, we could expect that the system could become unstable and that solutions can be periodic. This is the case for inhibitory integrate-and-fire models; these are PDEs describing neural networks, with strong relaxation properties to a steady state. It is well-known, see [16] for instance, that delays can generate a spontaneous activity, i.e. periodic solutions.

However, we did not observe such a behaviour in numerical simulations we conducted. This is confirmed by the stability analysis of a simplified equation.

The numerics have been performed in Matlab with parameters as follows. As initial data, we use *S*(*t*=0)=5 and *n*(*t*=0)=*C*_{mass} *e*^{−200(x−0.5)2}, where *C*_{mass} is chosen such that the initial mass in the computational domain is equal to 5. We set *R*:=20(−0.6+0.2*S*−(*x*−.5)^{2}) and *Q*(*ρ*,*S*):=8.5−(0.5+*ρ*)*S*. The equation is solved by an implicit–explicit finite-difference method on a grid consisting of 1000 points (time step *dt*=4⋅10^{−4}). Figure 1 shows oscillations of *ρ*(*t*) and *I*(*t*). Moreover, numerically, it seems that

### Remark 3.2

We can rewrite (3.10) and (3.8) as
*β*, small *A* has real eigenvalues, whereas for *β* large, it has complex eigenvalues. Therefore, our method cannot work for *β* large. One direction to extend the result would be to work directly on the system (3.10)–(3.8).

## 4. Perspectives and open questions

We have proved long-term convergence to an ESD for a general model of a chemostat where the nutrient concentration does not immediately equilibrate to the population dynamics. Our proof extends the proof based on total variation (TV) bounds developed in [2,7,15] and uses a fast (but not infinite) nutrient production measured by the small parameter *β*.

Surprisingly, the proof does not seem to give directly uniform TV bounds for *β*≈0. It does not seem to be possible with this approach to prove uniform bounds for the full range *β*∈[0,*β*_{0}] for some small *β*_{0}, which could be a first step to prove uniform convergence of *S*(*t*) for *β*→0.

There are several related problems which, usually, can be approached with the same method. One of them is the rare mutations/long-term behaviour described by the following extension of (1.5)–(1.7)
*ε*.

From the modelling side, the TV bounds giving long-term behaviour is not known in several examples of chemostat systems. An example is the quasi-stationary case with general uptake rate and rendering factor,

## Funding statement

This research was supported by the French ‘ANR blanche’ project Kibord: ANR-13-BS01-0004.

- Received February 2, 2014.
- Accepted April 4, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.