Global stability with selection in integro-differential Lotka-Volterra systems modelling trait-structured populations

ABSTRACT

We analyse the asymptotic behaviour of integro-differential equations modelling N populations in interaction, all structured by different traits. Interactions are modelled by non-local terms involving linear combinations of the total number of individuals in each population. These models have already been shown to be suitable for the modelling of drug resistance in cancer, and they generalize the usual Lotka-Volterra ordinary differential equations. Our aim is to give conditions under which there is persistence of all species. Through the analysis of a Lyapunov function, our first main result gives a simple and general condition on the matrix of interactions, together with a convergence rate. The second main result establishes another type of condition in the specific case of mutualistic interactions. When either of these conditions is met, we describe which traits are asymptotically selected.

KEYWORDS:

• Integro-differential equations
• population dynamics
• asymptotic analysis
• non-local model
• steady-states

1. Introduction

1.1. Biological motivations

We are interested in the evolution of N populations of individuals, each of which is structured by a continuous phenotypic trait, also called trait. In each species the phenotype models some continuous biological characteristics (such as the size of the individual, the concentration of a protein inside it, etc). We shall consider both interactions inside a given population and between the populations and we assume that mutations can be neglected. Mathematical modelling and analysis of such ecological scenarios is one purpose of the field of adaptive dynamics, a branch of mathematical biology which aims at describing evolution among a population of individuals, see [15,29,32] for an introduction to deterministic models.

The basis for our model stems from the logistic ODE $dN/dt=(r-dN)N$ where r is the net growth rate, dN the logistic death rate due to competition for nutrients and for space by direct or indirect inhibition of proliferation between individuals. Its natural extension to a density $n(t,x)$ of individuals of phenotype x (say in $[0,1]$) is a non-local logistic model (1) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}n(t,x)=\left(r\left(x\right)-d\left(x\right)\rho \left(t\right)\right)n(t,x),$(1) with $\rho \left(t\right):={\int }_0^{1}n(t,x)dx$ the total number of individuals.

Following [37], one might consider that these two terms combine both Darwinism selection through the intrinsic growth rate r, and Lamarckism induction through the logistic death rate since it depends on the environment. Note that these models can be derived from stochastic models at the individual level [7,16,21], and more generally measure-valued functions n can be considered [6,22]. The asymptotic behaviour of the previous model (1(1) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}n(t,x)=\left(r\left(x\right)-d\left(x\right)\rho \left(t\right)\right)n(t,x),$(1) ) and variants is analysed in [20,26,32], and one important property among others is that solution typically tend to concentrate on a few phenotypes, a convergence to Dirac masses in mathematical terms. These models are thus successful at representing the survival of only a few phenotypes, which we will refer to as selected.

To account for more complex interactions, one may want to consider a more general non-local logistic term than $d\left(x\right)\rho \left(t\right)={\int }_0^{1}d\left(x\right)n(t,y)dy$, in the form ${\int }_0^{1}K(x,y)n(t,y)dy$. The behaviour of such equations strongly depends on how localized the kernel is, and therefore so do the mathematical techniques to analyse them. Indeed, with an added diffusion term, a special case of this situation is the non-local Fisher-KPP equation. When the kernel is localized (small as soon as $|x-y|$ is large), then the solutions typically remain bounded independently of the mutation rate [23]: selection is no longer a feature. This property highlights how differently the solutions behave depending on the kernel, and that some choices are not appropriate for the ecological situation we are concerned with.

The non-locality $d\left(x\right)\rho \left(t\right)$ actually implies that interaction of an individual of phenotype x with other individuals of phenotype y has the same strength $d\left(x\right)$ regardless of y, because individuals do not only necessarily compete with those which have close phenotypes. As such, our choice can serve as a case study to understand the effect of a blind competition across individuals, essentially mediated by the total density.

In [20,26], the biological motivation to use this type of models comes from drug resistance in cancer: the phenotype represents the level of resistance to a given drug and the authors argue that this might be a better approach than a discrete description of the phenotype taking only a finite number of values. Indeed, it can be correlated to some continuous biological characteristics, such as the intracellular concentration of a detoxication molecule, the activity of detoxifying enzymes in metabolizing the administered drug, or drug efflux transporters eliminating the drug [9].

This model is further developed in [26] to incorporate the healthy cell population. Neglecting mutations, it becomes a system of two integro-differential equations of the form $\begin{array}{rl}\frac{\mathrm{\partial }{n}_1}{\mathrm{\partial }t}(t,x)& =\left(r_1\left(x\right)-d_1\left(x\right)\left(a_{11}{\rho }_1\left(t\right)+a_{22}{\rho }_2\left(t\right)\right)\right)n_1(t,x),\\ \frac{\mathrm{\partial }{n}_2}{\mathrm{\partial }t}(t,x)& =\left(r_2\left(x\right)-d_2\left(x\right)\left(a_{22}{\rho }_2\left(t\right)+a_{21}{\rho }_1\left(t\right)\right)\right)n_2(t,x),\end{array}$ with $a_{11}>0$, $a_{22}>0$, ${\rho }_1$, ${\rho }_2$ the total number of individuals in the cancer cell and healthy cell populations. The interspecific competition (between the two populations) is taken to be competitive, that is $a_{12}>0$, $a_{21}>0$, but below the intraspecific competition because each cell population is considered to belong to a different ecological niche: (2) $a_{12}<a_{11},a_{21}<a_{22}.$(2) In the context of a system arises the central question of persistence (whether asymptotically both species remain), complementing that of identifying which phenotypes are selected. With the additional difficulty of control terms to represent chemotherapeutic drugs, the asymptotic properties of this model are elucidated in [34], and assumption (2(2) $a_{12}<a_{11},a_{21}<a_{22}.$(2) ) happens to be crucial.

In the aforementioned work, carrying out the asymptotic analysis is instrumental for the application, which is the optimal control problem of minimizing the number of cancer cells through chemotherapy. Simulations suggest that the optimal strategy is to first let the cancer cell density concentrate on a sensitive phenotype, before using the maximum tolerated doses. The convergence to Dirac masses is at the cornerstone of the theoretical analysis of why this strategy happens to be optimal. It is also in sharp contradiction with the classical approach in the clinic, which relies on constant infusion of high doses.

These integro-differential models have therefore already proved their efficiency at helping understanding phenotypic heterogeneity in cancer. The mathematical results available for N=1 and N=2 for competitive interactions naturally call for generalizations on systems of interacting species with such non-local logical terms based on the total number of individuals. For instance, to study resistance in cancer, one may think also of different cancer subpopulations interacting with healthy cells and between them, each one of them being endowed with a specific drug resistance phenotype in a tumour 'bet hedging' strategy [4]. These generalizations, in turn, might help both unravel general principles about the underlying ecological processes, and develop new mathematical techniques to analyse them.

1.2. The model

We consider N populations structured by respective phenotypes $x\in X_i$, where $X_i$ is some compact subset of ${\mathbb{R}}^{p_i}$, with $p_i\in {\mathbb{N}}^{\ast }$, for $i=1,\dots ,N$. Although they model distinct quantities, we abusively denote all variables x to improve readability.

The model writes (3) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=\left(r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right)n_i(t,x),i=1,\dots ,N,$(3) where, for $i=1,\dots ,N$, $r_i$ and $d_i>0$ are functions in $L^{\mathrm{\infty }}\left(X_i\right)$, ${\rho }_i\left(t\right):={\int }_{X_i}{n}_i(t,x)dx$ is the total number of individuals in species i, and $a_{ij}\in \mathbb{R}$ are fixed (interaction) coefficients.

The initial conditions are (4) $n_i(0,\cdot )=n_i^{0}i=1,\dots ,N$(4) where each initial density $n_i^0$ is taken to be a non-negative function in $L^1\left(X_i\right)$. From now on, we will call these equations integro-differential Lotka-Volterra equations.

The matrix $A:=(a_{ij}{)}_{1\le i,j\le N}$, called matrix of interactions, describes the interactions between the populations: if $a_{ij}>0$, the species j acts positively on the species i, and negatively if $a_{ij}<0$. We will also say that the interaction between species i and j for $i\ne j$ is:

• mutualistic if $a_{ij}>0$ and $a_{ji}>0$,

• competitive if $a_{ij}<0$ and $a_{ji}<0$,

• predator-prey like if $a_{ij}{a}_{ji}<0$.

Finally, we will say that the equations are competitive (resp., mutualistic) if $a_{ij}<0$ (resp., $a_{ij}>0$) for all $i\ne j$.

Another interpretation of the equations is to see them as coupled logistic equations of the form (5) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=(r_i\left(x\right)-d_i\left(x\right)I_i\left(t\right))n_i(t,x),i=1,\dots ,N.$(5) In other words, the species i reacts to its environment through the non-local variable $I_i$ defined for $i=1,\dots ,N$ by (6) $I_i:=-\sum _{j=1}^N{a}_{ij}{\rho }_j.$(6) The terms $r_i\left(x\right)$ and $d_i\left(x\right)I_i$ respectively stand for the net proliferation rate and logistic death rate of individuals in species i, of phenotype x.

We will also use the notation $R_i(x,{\rho }_1,\dots ,{\rho }_N):=r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j$, with which the equations rewrite: (7) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=R_i\left(x,{\rho }_1\left(t\right),\dots ,{\rho }_N\left(t\right)\right)n_i(t,x),i=1,\dots ,N.$(7) These models generalize Lotka-Volterra ordinary differential equation (ODE) models [1]: if the functions $r_i$, $d_i$ are all constant (say equal to some $r_i$, and $d_i=1$), then after integration with respect to $x\in X_i$, the equations boil down to (8) $\frac{d}{dt}{\rho }_i\left(t\right)=\left(r_i+\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right){\rho }_i\left(t\right),i=1,\dots ,N,$(8) which we will from now on refer to as classical N-dimensional Lotka-Volterra equations. Thus, another advantage of a logistic term directly defined by ρ is that it makes our model more tractable with respect to the corresponding already well understood ODE models. Conversely, the integro-differential model can be seen as a perturbation of the ODE one where the individuals among a given population are allowed to have different proliferation and death rates depending on their phenotype x.

Our goal is to understand the asymptotic behaviour of the solutions to these equations, both in terms of convergence at the level of the total number of individuals ${\rho }_i$, and in terms of concentration at the level of the densities $n_i$. The first problem is usual in population dynamics while the second is specific to adaptive dynamics and consists of determining which traits asymptotically survive, taking over the whole population. These are then called Evolutionary Attractors, and the fact that it is the generic situation has been coined exclusion principle. Mathematically, this corresponds to a given density $n_i$ converging to a sum of Dirac masses. For one Dirac mass only, concentration writes, for some $x_0\in X_i$: (9) $n_i(t,\cdot )-{\rho }_i\left(t\right){\delta }_{x_0}⟶0$(9) as $t\to +\mathrm{\infty }$, in the weak sense of measures.

The more precise aim of this paper is to study the global asymptotic stability (GAS) of what we will call coexistence steady-states, namely of possible ${\rho }^{\mathrm{\infty }}$ with positive components (all species asymptotically survive) such that ρ converges to ${\rho }^{\mathrm{\infty }}$, because we will see how it determines on which phenotypes the densities concentrate. When it is possible, we will investigate the speed at which convergence and concentration occur. An interesting question within the scope of this paper is also to see if a result of that type is sharp, i.e. to compare the assumptions needed to obtain global asymptotic stability in our generalized setting to those known for classical Lotka-Volterra equations.

At this stage, we did not make any restrictive assumptions on the matrix A. However, it will be clear from the results recalled below in the ODE case and the ones presented in Section 2, that answers to the previous questions are available when interspecific interactions are low compared to the intraspecific ones (reminiscent of (2(2) $a_{12}<a_{11},a_{21}<a_{22}.$(2) )). Thus, we are covering the ecological scenario of each species i having its own niche, but inside which competition (if $a_{ii}<0$) is blind because of the term $a_{ii}{\rho }_i$.

Notations. In what follows, ${\mathbb{R}}_{>0}^N$ will stand for the positive orthant in ${\mathbb{R}}^N$, the set of vectors whose components are all positive, and we will write x>y when $x-y\in {\mathbb{R}}_{>0}^N$. We will also use the usual ordering on the set of symmetric matrices: for A a real symmetric matrices, we denote $A\ge 0$ (resp., A>0) when A is positive semidefinite (resp., positive definite). Finally, ${\mathcal{M}}^1\left(X\right)$ will denote the set of Radon measures supported in X.

1.3. State of the art

Classical Lotka-Volterra equations. The ODE system (8(8) $\frac{d}{dt}{\rho }_i\left(t\right)=\left(r_i+\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right){\rho }_i\left(t\right),i=1,\dots ,N,$(8) ) has been extensively studied, dating back to the pioneering works of Lotka and Volterra for two populations of preys and predators [28,36]. Since then, many contributions to the analysis of steady-states and their stability have been made, and we refer to [31] for an introduction and to [1] for a review.

Regarding the global asymptotic stability of coexistence steady-states, a very well-known result due to Goh [19] states a simple and very general condition on the matrix $A=(a_{ij}{)}_{1\le i,j\le N}$ which ensures convergence to the (unique) coexistence steady-state:

Theorem 1.1

[19]

Assume that the equation $A\rho +r=0$ $($where $r\in {\mathbb{R}}^N$ and $\rho \in {\mathbb{R}}^N$ are the vectors $(r_i{)}_{1\le i\le N}$ and $\left({\rho }_i{)}_{1\le i\le N}\right)$ has a solution ${\rho }^{\mathrm{\infty }}$ in ${\mathbb{R}}_{>0}^N$. If there exists a diagonal matrix D>0 such that $A^{T}D+DA<0,$ then ${\rho }^{\mathrm{\infty }}$ is GAS in ${\mathbb{R}}_{>0}^N$ $($and hence is the unique coexistence steady-state$)$ for system (8(8) $\frac{d}{dt}{\rho }_i\left(t\right)=\left(r_i+\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right){\rho }_i\left(t\right),i=1,\dots ,N,$(8) ).

A result also worth stating is that the mere existence of a unique coexistence steady-state is not enough for it to be GAS. Other steady-states on the boundary of ${\mathbb{R}}_{>0}^N$ can attract trajectories even in dimension N=2. Another possibility is the occurrence of chaotic behaviour even in low dimension as evidenced in [35] for N=3. Finally, we mention the more recent work [13], where the authors tackle the problem of GAS for some type of Lotka-Volterra ODEs with mutations. In particular, they obtain GAS of the coexistence steady-state in the case where the logistic variables $I_i$, $i=1,\dots ,N$ all coincide, that is, when they are equal to some variable $I:=\sum _{j=1}^N{a}_j{\rho }_j\left(t\right)$. In such a case, it is proved that the convergence to the equilibrium is exponential. The result of GAS is also extended to perturbations of this specific case.

Integro-differential Lotka-Volterra equations. The first question for such equations is the existence of a solution for all positive times. This obviously does not hold true in full generality since the ODE $y^{\prime }=y^2$ is a particular case. Let us first state an existence and uniqueness theorem.

Theorem 1.2

Assume that for a given $n^0\in \prod _{i=1}^N{L}^1\left(X_i\right),$ $n^0\ge 0,$ there exists $0<{\rho }^{sup}$ such that we have an a priori upper bound $\rho \left(t\right)\le {\rho }^{sup}$ for the functions ${\rho }_i$ whenever they are defined. Then the Cauchy problem (3(3) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=\left(r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right)n_i(t,x),i=1,\dots ,N,$(3) ) and (4(4) $n_i(0,\cdot )=n_i^{0}i=1,\dots ,N$(4) ) has a unique solution $n=(n_i{)}_{1\le i\le N},$ $n\ge 0,$ in $C\left(\right[0,+\mathrm{\infty }),\prod _{i=1}^N{L}^1(X_i\left)\right)$.

The proof is a straightforward generalization of that given in [32, Theorem 2.4] for N=1, relying on the Banach Fixed Point Theorem.

In the case of a single equation, the asymptotic behaviour is well understood. For N=1, assuming $a_{11}<0$ to avoid blow-up, the equation is simply $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_1(t,x)=\left(r_1\left(x\right)-d_1\left(x\right){\rho }_1\left(t\right)\right)n_1(t,x)$where, without loss of generality, we have set $a_{11}=-1$. The first result is that ${\rho }_1$ converges.

Theorem 1.3

Assume some regularity on $X_1,$ $r_1,$ $d_1,$ and $r_1>0,$ Then, for any positive continuous initial condition $n_1^0,$ ${\rho }_1$ the function $t\mapsto {\rho }_1\left(t\right)$ is well defined on $[0,+\mathrm{\infty })$ and converges to ${\rho }_1^M:=\underset{x\in X_1}{max}\left(r_1\left(x\right)/d_1\left(x\right)\right)$ as $t\to +\mathrm{\infty }$.

This, in turn, completely determines where $n_1$ concentrates.

Corollary 1.1

Under the previous hypotheses, $n_1\left(t\right),$ viewed as a Radon measure on $X_1,$ concentrates on the set $\left\{x\in X_1,r_1\left(x\right)-d_1\left(x\right){\rho }_1^M=0\right\}$ as $t\to +\mathrm{\infty }$. If this set is reduced to some $x_1^{\mathrm{\infty }}$, we obtain in particular $n_1(t,\cdot )⟶{\rho }_1^M{\delta }_{x_1^{\mathrm{\infty }}}$ as $t\to +\mathrm{\infty }$ in ${\mathcal{M}}^1\left(X_1\right)$ equipped with its usual weak star topology.

A proof of this result can be found in [34] in the special case $X_1=[0,1]$, and it relies on proving that ${\rho }_1$ is a bounded variation (BV) function on $[0,+\mathrm{\infty })$. Let us stress that when the set on which $n_1$ concentrates is not reduced to a singleton, the steady-state (at the level of $n_1$) is not unique. For example, if the set is made of two points, the repartition of the limiting density on these two points depends on the initial condition, see for example [12]. This is why for this equation and the general equations considered here, there is no hope in proving general GAS results directly at the level of the densities $n_i$.

For a general logistic term $\left({\int }_{X}K\right(x,y\left)n\right(t,y\left)dy\right)n(t,x)$ and a single equation, the asymptotic behaviour is also analysed in detail in both [14] and [24]. In the latter, under some suitable assumptions on the kernel K, a Lyapunov functional is used to prove that some measure is GAS, in a very specific sense depending on K. Similar results can be found in [8], where their counterpart for competitive classical Lotka-Volterra equations are also discussed.

In the case of integro-differential systems, however, much less is known about the asymptotic behaviour. A Lyapunov functional inspired by [24] has been used successfully in [34] to prove GAS for a competitive system of two populations which writes exactly as (3(3) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=\left(r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right)n_i(t,x),i=1,\dots ,N,$(3) ). We also mention [5] where an integro-differential system of two populations is analysed, and whose form does not fit in our framework. A particular triangular coupling structure allows the authors to perform an asymptotic analysis.

The paper is organized as follows. In Section 2, we explain how coexistence steady-states can be identified, allowing us to state rigorously what we mean by GAS for system (3(3) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=\left(r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right)n_i(t,x),i=1,\dots ,N,$(3) ). We explain why, under the hypothesis of GAS, only some phenotypes are generically selected, and how to compute them. Then, we present the two main results about GAS for such equations. Section 3 is devoted to the proof of the first result, which applies for any type of interactions and relies on analysing a suitably designed Lyapunov functional. In the specific case of mutualistic interactions, our second main result gives alternative conditions sufficient for GAS. It is presented in Section 4. In Section 5, we conclude with several comments and open questions.

2. Possible coexistence steady-states and main results

For the rest of the article, we will work with the following assumptions: (10) $r_i,d_i,n_i^0\in C\left(X_i\right),n_i^0>0\mathrm{for}i=1,\dots ,N.$(10) This will simplify statements, but we will be more specific below as to which data our results generalize.

2.1. Analysis of coexistence steady-states

In the context of this system of integro-differential equations, the expression ‘GAS in ${\mathbb{R}}_{>0}^N$’ must be defined. By that, we mean that there exists ${\rho }^{\mathrm{\infty }}>0$ such that, whatever the positive continuous initial conditions $n_i^0$ are, ${\rho }_i$ converges to ${\rho }_i^{\mathrm{\infty }}$ for all i.

First, let us explain how to compute the possible steady-states at the level of ρ, i.e. possible limits ${\rho }^{\mathrm{\infty }}>0$ for positive continuous initial conditions. We will work with the following topological assumption on the sets $X_i$: (11) $\mathrm{\forall }x\in \mathrm{\partial }{X}_i,\mathrm{\forall }\eta >0,{\lambda }_{p_i}\left(B(x,\eta )\cap X_i\right)>0,$(11) where ${\lambda }_{p_i}$ stands for the Lebesgue measure on ${\mathbb{R}}^{p_i}$ and $B(x,\eta )$ for the open ball of centre x and radius η.

Assume that each ${\rho }_i$ converges to some ${\rho }_i^{\mathrm{\infty }}>0$, in which case the exponential behaviour of $n_i(t,x)$ is asymptotically governed by $r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j^{\mathrm{\infty }}$, the sign of which we can analyse as follows.

• If this quantity is positive for some $x_0$, let us prove that $n_i(t,x)$ blows up in its neighbourhood, leading to the explosion of ${\rho }_i$. If $r_i\left(x_0\right)+d_i\left(x_0\right)\sum _{j=1}^N{a}_{ij}{\rho }_j^{\mathrm{\infty }}>0$, there exists $\eta >0$ such that by continuity $r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j^{\mathrm{\infty }}>0$ on $\left(B\right(x_0,\eta )\cap X_i)$, and then ${\lambda }_{p_i}\left(B\right(x_0,\eta )\cap X_i)>0$ whether $x_0\in int\left(X_i\right)$ or also if $x_0\in \mathrm{\partial }{X}_i$ thanks to (11(11) $\mathrm{\forall }x\in \mathrm{\partial }{X}_i,\mathrm{\forall }\eta >0,{\lambda }_{p_i}\left(B(x,\eta )\cap X_i\right)>0,$(11) ). For $\epsilon >0$ small enough and t large enough (say $t\ge t_0$) such that $r_i\left(x_0\right)+d_i\left(x_0\right)\sum _{j=1}^N{a}_{ij}{\rho }_j>\epsilon$, we can write: $\begin{array}{rl}{\rho }_i\left(t\right)& \ge {\int }_{B(x_0,\eta )\cap X_i}{n}_i(t,x)dx\\ & \ge {\int }_{B(x_0,\eta )\cap X_i}{n}_i(t_0,x)e^{{\int }_{t_0}^t{R}_i\left(x,{\rho }_1\left(s\right),\dots ,{\rho }_N\left(s\right)\right)ds}dx\\ & \ge {\lambda }_{p_i}\left(B(x_0,\eta )\cap X_i\right)\left(\underset{B(x_0,\eta )\cap X_i}{inf}{n}_i(t_0,x)\right)e^{\epsilon (t-t_0)},\end{array}$ with the right-hand side blowing up as $t\to +\mathrm{\infty }$, which cannot be compatible with the convergence of ${\rho }_i$.

• If $r_i+d_i\sum _{j=1}^N{a}_{ij}{\rho }_j^{\mathrm{\infty }}$ is negative globally on $X_i$, this clearly implies the extinction of species i, which is also incompatible with the convergence of ${\rho }_i$ to a positive limit.

Remark 2.1

It is possible to relax the regularity on both the sets $X_i$ and the data $r_i$ and $d_i$ by working only with points which are both Lebesgue points of $r_i/d_i$ and of Lebesgue density 1 for $X_i$, see [17]. If the functions $r_i/d_i$ are in $L^1\left(X_i\right)$, one can indeed check that $r_i+d_i\sum _{j=1}^N{a}_{ij}{\rho }_j^{\mathrm{\infty }}\le 0$ $a.e.$ on $X_i$.

The previous results motivate the following definition: $I_i^{\mathrm{\infty }}:=\underset{x\in X_i}{max}\frac{r_i\left(x\right)}{d_i\left(x\right)},i=1,\dots ,N.$ With this definition, a steady-state ${\rho }^{\mathrm{\infty }}>0$ exists if and only if the following assumption holds: (12) $\mathrm{the}\mathrm{equation}A\rho +I^{\mathrm{\infty }}=0\mathrm{has}\mathrm{a}\mathrm{solution}{\rho }^{\mathrm{\infty }}\mathrm{in}{\mathbb{R}}_{>0}^N,$(12) which we assume from now on.

The previous computations also show that $n_i$ vanishes where $r_i\left(x\right)-d_i\left(x\right)I_i^{\mathrm{\infty }}<0$, which implies the following result:

Proposition 2.1

Assume that assumption (12(12) $\mathrm{the}\mathrm{equation}A\rho +I^{\mathrm{\infty }}=0\mathrm{has}\mathrm{a}\mathrm{solution}{\rho }^{\mathrm{\infty }}\mathrm{in}{\mathbb{R}}_{>0}^N,$(12) ) holds, and that ρ converges to the coexistence steady-state ${\rho }^{\mathrm{\infty }}$. Then, $n_i\left(t\right),$ viewed as a Radon measure, concentrates on the set $K_i:=\{x\in X_i,r_i\left(x\right)-d_i\left(x\right)I_i^{\mathrm{\infty }}=0\}$ as $t\to +\mathrm{\infty },$ for all $i=1,\dots ,N$.

If, for some i, $K_i$ is reduced to some $x^{\mathrm{\infty }},$ we obtain in particular $n_i(t,\cdot )⟶{\rho }_i^{\mathrm{\infty }}{\delta }_{x_i^{\mathrm{\infty }}}$ as $t\to +\mathrm{\infty }$ in ${\mathcal{M}}^1\left(X_i\right)$.

Densities $n_i$ of total mass ${\rho }_i^{\mathrm{\infty }}$ and concentrated on $K_i$ are called Evolutionary Stable Distributions (ESD) in [24].

Remark 2.2

In the hypothesis of global existence and convergence of ρ towards ${\rho }^{\mathrm{\infty }}$, the previous reasoning actually shows that the concentration is ensured as soon as $n_i^0\in L^1\left(X_i\right)$ is bounded by below by a positive constant on a neighbourhood of one of the points of $K_i$. For more general hypotheses ensuring concentration, we refer to [24].

Remark 2.3

If all the sets $K_i$ are reduces to some singletons $x_i^{\mathrm{\infty }}$, then the dynamics of ρ are asymptotically governed by classical Lotka-Volterra equations concentrated in $(x_1^{\mathrm{\infty }},\dots ,x_N^{\mathrm{\infty }})$, namely $\frac{d}{dt}{\rho }_i\left(t\right)\simeq \left(r_i(x_i^{\mathrm{\infty }})+d_i(x_i^{\mathrm{\infty }})\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right){\rho }_i\left(t\right),i=1,\dots ,N,$ as t goes to $+\mathrm{\infty }$. For a precise statement, see [34].

2.2. Analysis of coexistence steady-states

Our first approach to prove GAS is to further generalize the approach of [34] in dimension N. The main idea is to mix a Lyapunov functional which is inspired by the one designed in [24] and the Lyapunov functional used for classical Lotka-Volterra equations [19], which is the key tool to obtain Theorem 1.1. With some mild regularity assumptions on the data, we obtain the following result:

Theorem 2.1

Assume (12(12) $\mathrm{the}\mathrm{equation}A\rho +I^{\mathrm{\infty }}=0\mathrm{has}\mathrm{a}\mathrm{solution}{\rho }^{\mathrm{\infty }}\mathrm{in}{\mathbb{R}}_{>0}^N,$(12) ) and that there exists a diagonal matrix D>0 such that DA is symmetric and DA<0. Then the solution to the Cauchy problem (3(3) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=\left(r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right)n_i(t,x),i=1,\dots ,N,$(3) ) and (4(4) $n_i(0,\cdot )=n_i^{0}i=1,\dots ,N$(4) ) is globally defined. Furthermore, the solution ${\rho }^{\mathrm{\infty }}$ to $A\rho +I^{\mathrm{\infty }}=0$ is GAS $($and hence, unique$)$.

We emphasize that there is no assumption on the type of interactions, i.e. on the sign of the coefficients of A. However, a consequence of this result is that A must be such that $a_{ii}{a}_{jj}>a_{ij}{a}_{ji}$ for all i, j. For this result to apply, interactions must therefore be stronger inside species than between them.

We also remark that our hypothesis is exactly the one exhibited in [8] for competitive classical Lotka-Volterra equations. The analysis of the Lyapunov functional allows to determine a speed at which convergence to ${\rho }^{\mathrm{\infty }}$ and concentration on a given set of phenotypes occur. In dimension 2, we also analyse more deeply the link between this condition and the one for classical Lotka-Volterra equations, which in most interesting cases happen to be equivalent.

Our second main result focuses on the special case of mutualistic interactions, and an informal statement of the theorem is the following.

Theorem 2.2

Assume (12(12) $\mathrm{the}\mathrm{equation}A\rho +I^{\mathrm{\infty }}=0\mathrm{has}\mathrm{a}\mathrm{solution}{\rho }^{\mathrm{\infty }}\mathrm{in}{\mathbb{R}}_{>0}^N,$(12) ), that for $i=1,\dots ,N,$ $r_i>0$ and that for some explicit $0<c_i<C_i,$ the matrix $\hat{A}$ defined by ${\hat{a}}_{ii}:=c_i{a}_{ii}$ and ${\hat{a}}_{ij}=C_i{a}_{ij}$ for $i\ne j$ is Hurwitz. Then the solution to the Cauchy problem (3(3) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=\left(r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right)n_i(t,x),i=1,\dots ,N,$(3) ) and (4(4) $n_i(0,\cdot )=n_i^{0}i=1,\dots ,N$(4) ) is globally defined. Furthermore, the solution ${\rho }^{\mathrm{\infty }}$ to $A\rho +I^{\mathrm{\infty }}=0$ is GAS.

Again, this applies to the case of interspecific interactions being higher that intraspecific ones, because a Hurwitz matrix is a matrix such that all its eigenvalues have negative real part and it has to do with diagonally dominant matrices (see Section 4).

Because of the hypothesis of mutualism, the system is cooperative, and sub and supersolution techniques can succeed. More precisely, it is possible to prove that all functions ${\rho }_i$ are BV on $[0,+\mathrm{\infty })$ and this implies their convergence.

3. General interactions

3.1. Proof of the main theorem

In this section, we need slightly more regularity for the data, namely that the functions are Lipschitz continuous: (13) $\mathrm{for}i=1,\dots ,N,r_i,d_i\in C^{0,1}\left(X_i\right).$(13) We can now restate the first theorem:

Theorem 3.1

Assume (12(12) $\mathrm{the}\mathrm{equation}A\rho +I^{\mathrm{\infty }}=0\mathrm{has}\mathrm{a}\mathrm{solution}{\rho }^{\mathrm{\infty }}\mathrm{in}{\mathbb{R}}_{>0}^N,$(12) ) and (13(13) $\mathrm{for}i=1,\dots ,N,r_i,d_i\in C^{0,1}\left(X_i\right).$(13) ). Assume that there exists a diagonal matrix D>0 such that DA is symmetric and DA<0. Then the solution to the Cauchy problem (3(3) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=\left(r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right)n_i(t,x),i=1,\dots ,N,$(3) ) and (4(4) $n_i(0,\cdot )=n_i^{0}i=1,\dots ,N$(4) ) is globally defined.

Furthermore, the solution ${\rho }^{\mathrm{\infty }}$ to $A\rho +I^{\mathrm{\infty }}=0$ is GAS with (14) $\rho \left(t\right)-{\rho }^{\mathrm{\infty }}=\mathrm{O}\left({\left(\frac{\ln \left(t\right)}{t}\right)}^{1/2}\right).$(14) Concentration of a given $n_i$ occurs at speed $\mathrm{O}\left(\ln \left(t\right)/t\right),$ in the following sense: (15) ${\int }_{X_i}{m}_i\left(x\right)R_i\left(x,{\rho }_1^{\mathrm{\infty }},\dots ,{\rho }_N^{\mathrm{\infty }}\right)n_i(t,x)dx=\mathrm{O}\left(\frac{\ln \left(t\right)}{t}\right).$(15) In particular, if $K_i$ is reduced to a singleton $x_i^{\mathrm{\infty }},$ then (16) $\mathrm{\forall }\epsilon >0,{\int }_{X_i\setminus B\left(x_i^{\mathrm{\infty }},\epsilon \right)}{n}_i(t,x)dx=\mathrm{O}\left(\frac{\ln \left(t\right)}{t}\right).$(16)

Proof.

First step: definition of the Lyapunov functional.

From (12(12) $\mathrm{the}\mathrm{equation}A\rho +I^{\mathrm{\infty }}=0\mathrm{has}\mathrm{a}\mathrm{solution}{\rho }^{\mathrm{\infty }}\mathrm{in}{\mathbb{R}}_{>0}^N,$(12) ), Evolutionary Stable Densities exist and we can pick one of them: for $i=1,\dots ,N$, we choose any measure $n_i^{\mathrm{\infty }}$ in ${\mathcal{M}}^1\left(X_i\right)$ satisfying $n_i^{\mathrm{\infty }}\left(X_i\right)={\rho }_i^{\mathrm{\infty }}$, which is furthermore concentrated on $K_i$, i.e. (17) $\mathrm{supp}\left(n_i^{\mathrm{\infty }}\right)\subset K_i.$(17) We abusively write integration of functions g against measures μ as ${\int }_{X}g\left(x\right)\mu \left(x\right)dx$. We also set $m_i:=1/d_i$ and define N functions $V_i$ by $V_i\left(t\right):={\int }_{X_i}{m}_i\left(x\right)\left[n_i^{\mathrm{\infty }}\left(x\right)\ln \left(\frac{1}{n_i(t,x)}\right)+\left(n_i(t,x)-n_i^{\mathrm{\infty }}\left(x\right)\right)\right]dx.$ In what follows, we consider the following Lyapunov functional: $V\left(t\right):=\sum _{i=1}^N{\lambda }_i{V}_i\left(t\right)$ where the positive constants ${\lambda }_i$ are to be chosen later. The diagonal matrix of diagonal entries ${\lambda }_1,\dots ,{\lambda }_N$ is denoted by D.

Second step: computation and sign of the derivative.

For any i, we compute $\begin{array}{rl}\frac{d{V}_i}{dt}& ={\int }_{X_i}{m}_i\left(x\right)\left[-n_i^{\mathrm{\infty }}\left(x\right)\frac{{\mathrm{\partial }}_t{n}_i(t,x)}{n_i(t,x)}+{\mathrm{\partial }}_t{n}_i(t,x)\right]dx\\ & ={\int }_{X_i}{m}_i\left(x\right)R_i\left(x,{\rho }_1,\dots ,{\rho }_N\right)\left[n_i(t,x)-n_i^{\mathrm{\infty }}\left(x\right)\right]dx\\ & ={\int }_{X_i}{m}_i\left(x\right)\left(R_i\left(x,{\rho }_1,\dots ,{\rho }_N\right)-R_i\left(x,{\rho }_1^{\mathrm{\infty }},\dots ,{\rho }_N^{\mathrm{\infty }}\right)\right)\left[n_i(t,x)-n_i^{\mathrm{\infty }}\left(x\right)\right]dx\\ & +{\int }_{X_i}{m}_i\left(x\right)R_i\left(x,{\rho }_1^{\mathrm{\infty }},\dots ,{\rho }_N^{\mathrm{\infty }}\right)\left[n_i(t,x)-n_i^{\mathrm{\infty }}\left(x\right)\right]dx.\end{array}$ The definition of $m_i$ implies that the first term simplifies as follows ${\int }_{X_i}{m}_i\left(x\right)d_i\left(x\right){\left[A\left(\rho -{\rho }^{\mathrm{\infty }}\right)\right]}_i\left[n_i(t,x)-n_i^{\mathrm{\infty }}\left(x\right)\right]dx={\left[A\left(\rho -{\rho }^{\mathrm{\infty }}\right)\right]}_i\left({\rho }_i-{\rho }_i^{\mathrm{\infty }}\right).$ For the second term, the choice (17(17) $\mathrm{supp}\left(n_i^{\mathrm{\infty }}\right)\subset K_i.$(17) ) leads to $B_i\left(t\right):={\int }_{X_i}{m}_i\left(x\right)R_i\left(x,{\rho }_1^{\mathrm{\infty }},\dots ,{\rho }_N^{\mathrm{\infty }}\right)n_i(t,x)dx.$ The functions $B_i$ are all non-positive by definition of ${\rho }^{\mathrm{\infty }}$.

Defining the vector $u:=\rho -{\rho }^{\mathrm{\infty }}$, we arrive at: $\begin{array}{rl}\frac{dV}{dt}& =\sum _{i=1}^N{\lambda }_i{\left[A\left(\rho -{\rho }^{\mathrm{\infty }}\right)\right]}_i\left({\rho }_i-{\rho }_i^{\mathrm{\infty }}\right)+\sum _{i=1}^N{\lambda }_i{B}_i\\ & =\sum _{i=1}^N{\lambda }_i(Au{)}_i{u}_i+\sum _{i=1}^N{\lambda }_i{B}_i.\end{array}$ Thus, we end up with the expression (18) $\frac{dV}{dt}=u^T\left(DA\right)u+\sum _{i=1}^N{\lambda }_i{B}_i.$(18) Since the antisymmetric part of DA does not play any role, this can also be expressed $\frac{dV}{dt}=\frac{1}{2}{u}^T(A^{T}D+DA)u+\sum _{i=1}^N{\lambda }_i{B}_i.$ By assumption, $M:=A^{T}D+DA<0$, from which we deduce $dV/dt\le 0$. Furthermore, the convergence of the term $u^{T}Mu$ to 0 is equivalent to that of ρ to ${\rho }^{\mathrm{\infty }}$. However, we do not have the usual property $V\ge 0$ for Lyapunov functions, so that we cannot yet conclude.

Third step: estimates on $dV/dt$.

Let $G:=\frac{1}{2}{u}^{T}Mu+2\sum _{i=1}^N{\lambda }_i{B}_i.$ We are going to show that G is non-decreasing.

We denote by $⟨u,v⟩$ the canonical scalar product of two vectors u and v in ${\mathbb{R}}^N$. $\begin{array}{rl}\frac{d}{dt}\left(u^T\left(DA\right)u\right)& =\frac{d}{dt}⟨\left(DA\right)u,u⟩\\ & =⟨\left(DA\right)\frac{du}{dt},u⟩+⟨\left(DA\right)u,\frac{du}{dt}⟩.\end{array}$ For $i=1,\dots ,N$, the derivative of $B_i$ is given by $\begin{array}{rl}\frac{d{B}_i}{dt}& ={\int }_{X_i}{m}_i\left(x\right)R_i\left(x,{\rho }_1^{\mathrm{\infty }},\dots ,{\rho }_N^{\mathrm{\infty }}\right)R_i\left(x,{\rho }_1,\dots ,{\rho }_N\right)n_i(t,x)dx\\ & ={\int }_{X_i}{m}_i\left(x\right)R_i^2\left(x,{\rho }_1,\dots ,{\rho }_N\right)n_i(t,x)dx\\ & +{\int }_{X_i}{m}_i\left(x\right)\left[R_i\left(x,{\rho }_1^{\mathrm{\infty }},\dots ,{\rho }_N^{\mathrm{\infty }}\right)-R_i\left(x,{\rho }_1,\dots ,{\rho }_N\right)\right]R_i\left(x,{\rho }_1,\dots ,{\rho }_N\right)n_i(t,x)dx\\ & \ge [A\left({\rho }^{\mathrm{\infty }}-\rho \right)]{\int }_{X_i}{R}_i\left(x,{\rho }_1,\dots ,{\rho }_N\right)n_i(t,x)dx\\ & =-(Au{)}_i{\left(\frac{du}{dt}\right)}_i\end{array}$ leading to the bound $\begin{array}{rl}\frac{d}{dt}\left(\sum _{i=1}^N{\lambda }_i{B}_i\right)& \ge -\sum _{i=1}^N{\lambda }_i(Au{)}_i{\left(\frac{du}{dt}\right)}_i\\ & =-⟨\left(DA\right)u,\frac{du}{dt}⟩.\end{array}$ Put together, these estimates yield: $\begin{array}{rl}\frac{dG}{dt}& \ge ⟨\left(DA\right)\frac{du}{dt},u⟩+⟨\left(DA\right)u,\frac{du}{dt}⟩-2⟨\left(DA\right)u,\frac{du}{dt}⟩\\ & =⟨\left(DA\right)\frac{du}{dt},u⟩-⟨\left(DA\right)u,\frac{du}{dt}⟩.\end{array}$ The last expression is equal to 0 if DA is symmetric, in which case G is non-decreasing as claimed.

The assumptions that DA is symmetric and that $A^{T}D+DA<0$ are equivalent to the assumption made for the theorem: DA is supposed to be a symmetric negative definite matrix.

As a consequence of the monotonicity of G, we get $u^T(-DA)u\le -G\left(0\right)$ for all t. The left-hand side is the square of some norm on ${\mathbb{R}}^N$, which means that ρ has to be bounded: these a priori bounds ensure the global definition of the solution to (3(3) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=\left(r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right)n_i(t,x),i=1,\dots ,N,$(3) ) and (4(4) $n_i(0,\cdot )=n_i^{0}i=1,\dots ,N$(4) ).

Fourth step: a lower estimate for V .

To estimate V from below, we need a uniform (in x) upper bound on the densities $n_i$. Because of the regularity assumption (13(13) $\mathrm{for}i=1,\dots ,N,r_i,d_i\in C^{0,1}\left(X_i\right).$(13) ), there exists C>0 such that: $\mathrm{\forall }i=1,\dots ,N,\mathrm{\forall }(x,y)\in X_i^2,R_i\left(y,{\rho }_1,\dots ,{\rho }_N\right)\ge R_i\left(x,{\rho }_1,\dots ,{\rho }_N\right)-C|x-y|.$ The constant C can be chosen to be independent of t since the functions ${\rho }_i$ are bounded.

This implies for a given i $\begin{array}{rl}{\int }_{X_i}{n}_i(t,y)dy& ={\int }_{X_i}{n}_i^0\left(y\right)\exp \left({\int }_0^t{R}_i\left(y,{\rho }_1,\dots ,{\rho }_N\right)ds\right)dy\\ & \ge {\int }_{X_i}\frac{n_i^0\left(y\right)}{n_i^0\left(x\right)}\left(n_i^0\left(x\right)\exp \left({\int }_0^t{R}_i\left(x,{\rho }_1,\dots ,{\rho }_N\right)\right)\right)\exp \left(-Ct|x-y|\right)dy\\ & \ge \frac{n_i(t,x)}{n_i^0\left(x\right)}{\int }_{X_i}\exp \left(-Ct|x-y|\right)dy.\end{array}$ Computing the integral, we write, thanks to the boundedness of ${\rho }_i$ $n_i^0$ and $n_i^0$ and (C has changed and is independent of t and x): for t large enough, $n_i(t,x)\le Ct$. The bound on V follows immediately: (19) $V\left(t\right)\ge -C\left(\ln \left(t\right)+1\right).$(19)

Fifth step: convergence.

G bounds $dV/dt$ from above: $dV/dt\le \frac{1}{2}G$. Thus $V\left(t\right)-V\left(0\right)\le \frac{1}{2}{\int }_0^{t}G\left(s\right)ds\le \frac{1}{2}tG\left(t\right)$ thanks to the third step.

We can now write $G\left(t\right)\ge -C\left(\right(\ln \left(t\right)+1\left)/t\right)$: $G=\mathrm{O}\left(\ln \left(t\right)/t\right)$, consequently, each non-positive term it is composed of also vanishes like $\mathrm{O}\left(\ln \left(t\right)/t\right)$ as $t\to +\mathrm{\infty }$.

In other words, $\frac{1}{2}{u}^{T}Mu$ and each $B_i$ onverge to 0 as as well $\mathrm{O}\left(\ln \left(t\right)/t\right)$. This is nothing but the two first statements (14(14) $\rho \left(t\right)-{\rho }^{\mathrm{\infty }}=\mathrm{O}\left({\left(\frac{\ln \left(t\right)}{t}\right)}^{1/2}\right).$(14) ) and (15(15) ${\int }_{X_i}{m}_i\left(x\right)R_i\left(x,{\rho }_1^{\mathrm{\infty }},\dots ,{\rho }_N^{\mathrm{\infty }}\right)n_i(t,x)dx=\mathrm{O}\left(\frac{\ln \left(t\right)}{t}\right).$(15) ).

For the last statement, we fix i and $\epsilon >0$. We denote $h_i:=-m_i{R}_i(\cdot ,{\rho }_1^{\mathrm{\infty }},\dots ,{\rho }_N^{\mathrm{\infty }})$, which is non-negative on $X_i$, and by assumption vanishes at $x_i^{\mathrm{\infty }}$ only. We choose a>0 small enough such that $a{1}_{X_i\setminus B(x_i^{\mathrm{\infty }},\epsilon )}\le h$ on $X_i$. This enables us to write ${\int }_{X_i\setminus B\left(x_i^{\mathrm{\infty }},\epsilon \right)}{n}_i(t,x)dx\le \frac{1}{a}{\int }_{X_i}{m}_i\left(x\right)R_i\left(x,{\rho }_1^{\mathrm{\infty }},\dots ,{\rho }_N^{\mathrm{\infty }}\right)n_i(t,x)dx=\mathrm{O}\left(\frac{\ln \left(t\right)}{t}\right).$

Remark 3.1

The first interesting fact is that the convergence rate of G to 0, in $\mathrm{O}\left(\ln \left(t\right)/t\right)$, is almost optimal in many cases. Indeed, if the sets $K_i$ are reduced to singletons, there cannot exist any $\alpha >1$ such that this sum vanishes like $\mathrm{O}\left(1/t^{\alpha }\right)$. This comes from the fact that if it were true, $dV/dt$ would be integrable on the half-line, which would imply the convergence of V . This is not possible since each $V_i$ has to go to $-\mathrm{\infty }$ as t goes to $+\mathrm{\infty }$.

This might seem contradictory with the exponential convergence rates obtained in [13] for some classical Lotka-Volterra equations, but the Lyapunov functional gives us information on the speed of both phenomena in the sense defined above (through the function G) and it does not say whether one of the two is faster.

3.2. Sharpness in dimension 2

It is clear that if we can find D>0 diagonal such that DA is symmetric and DA<0, then $A^{T}D+DA<0$. The condition that DA should be symmetric is constraining, especially if $N\ge 3$ in which case it imposes some polynomial equalities on the coefficients of the matrix A. In dimension 2, however, the result is sharp in various contexts, as stated in the following proposition.

Proposition 3.1

Assume N=2, $a_{11}<0,$ $a_{22}<0$ and $a_{12}{a}_{21}>0$. Then the following conditions are equivalent.

1. there exists D>0 diagonal such that DA is symmetric and $DA<0;$

2. there exists D>0 diagonal such that $A^{T}D+DA<0;$

3. $det\left(A\right)>0$.

Proof.

(i) implies (ii) as noticed before.

Now, let us assume (ii) and compute $M:=A^{T}D+DA=\left(\begin{array}{cc}2{\lambda }_1{a}_{11}& {\lambda }_1{a}_{12}+{\lambda }_2{a}_{21}\\ {\lambda }_1{a}_{12}+{\lambda }_2{a}_{21}& 2{\lambda }_2{a}_{22}\end{array}\right)$, which has positive determinant, i.e. $det\left(M\right)=4{\lambda }_1{\lambda }_2{a}_{11}{a}_{22}-({\lambda }_1{a}_{12}+{\lambda }_2{a}_{21}{)}^2>0$. If $det\left(A\right)\le 0$, $det\left(M\right)\le 4{\lambda }_1{\lambda }_2{a}_{12}{a}_{21}-({\lambda }_1{a}_{12}+{\lambda }_2{a}_{21}{)}^2=-({\lambda }_1{a}_{12}-{\lambda }_2{a}_{21}{)}^2\le 0$, a contradiction.

Now, if (iii) holds, we take ${\lambda }_1:=1/|a_{12}|$ and ${\lambda }_1:=1/|a_{21}|$ for which $DA=\left(\begin{array}{cc}\frac{a_{11}}{|a_{12}|}& \mathrm{sgn}\left(a_{12}\right)\\ \mathrm{sgn}\left(a_{21}\right)& \frac{a_{22}}{|a_{21}|}\end{array}\right)$ is clearly symmetric negative definite, whence (i).

4. Cooperative case

4.1. Some facts about Hurwitz matrices

We now focus on the cooperative case, i.e. on the case where all off-diagonal elements of A are non-negative. We will also assume that the diagonal elements are negative, since otherwise blow-up clearly occurs: there is intra-spectific competition inside any given species. We will say that such a matrix is cooperative.

In this case, we can hope for stronger results at the level of the integro-differential system because we can use sub and super-solution techniques. For our purpose, the following result on ODEs is sufficient.

Lemma 4.1

For T>0 $($possibly $T=+\mathrm{\infty }),$ let $f:[0,T)\times {\mathbb{R}}^N⟶{\mathbb{R}}^N$ be a continous function on $[0,T)\times {\mathbb{R}}^N,$ locally Lipschitz in $x\in {\mathbb{R}}^N$ uniformly in $t\in [0,T)$. Denoting $f(t,x):=\left(f_i\right(t,x_1,\dots ,x_N){)}_{1\le i\le N},$ further assume that for all $i=1,\dots ,N,$ $f_i$ is non-decreasing with $x_j$ for all $j\ne i$.

Assume that we have a solution z on $[0,T)$ of the following Cauchy problem: (20) $\begin{array}{r}\begin{array}{rl}\frac{dz}{dt}& =f(t,z)\\ z\left(0\right)& =z_0,\end{array}\end{array}$(20) where $z_0\in {\mathbb{R}}^N,$ and a function y subsolution to the previous Cauchy problem, i.e. (21) $\begin{array}{r}\begin{array}{rl}\frac{dy}{dt}& \le f(t,y)\\ y\left(0\right)& \le z_0.\end{array}\end{array}$(21) Then $y\left(t\right)\le z\left(t\right)$ on $[0,T)$.

When the matrix A is cooperative, it is possible to give an equivalent condition to the one required in Theorem 1.1 for GAS in classical Lotka-Volterra equations. Let us explain how, starting with the three following lemmas, the two first of which can be found in [1].

Lemma 4.2

Let A be a cooperative matrix. Then it is Hurwitz if and only if it is negatively diagonally dominant, i.e. if and only if there exists a vector v>0 such that $a_{ii}{v}_i+\sum _{j\ne i}{a}_{ij}{v}_j<0$ for $i=1,\dots ,N$.

This first lemma will be useful in its own right in this section. A consequence is that

Lemma 4.3

If A is cooperative and $r>0,$ $A\rho +r=0$ has a unique solution in ${\mathbb{R}}_{>0}^N$ if and only if A is Hurwitz.

Finally, it comes from the theory of M-matrices (see [33] for a review) that

Lemma 4.4

Let A be cooperative. Then A is Hurwitz if and only if there exists $D>0$ diagonal such that $A^{T}D+DA<0$.

An important consequence of these three lemmas is the following rephrasing of Theorem 1.1 for classical Lotka-Volterra equations.

Proposition 4.1

Assume that A is cooperative, r>0 and that the equations (8(8) $\frac{d}{dt}{\rho }_i\left(t\right)=\left(r_i+\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right){\rho }_i\left(t\right),i=1,\dots ,N,$(8) ) have a unique steady-state in ${\mathbb{R}}_{>0}^N$. Then the equations are globally defined and this steady-state is GAS.

Thus, in the context of cooperation between the species, the requirement that A is Hurwitz is somehow optimal to obtain a GAS coexistence steady-state, since it is already required to have its mere existence, a fact mentioned in [18]. We will mainly work with this characterization (rather than the equivalent one given by Lemma 4.4 which we used for a general interaction matrix A) because the next results will lead us to modify the matrix A: analysing whether it is still Hurwitz or not is easier than checking this equivalent condition.

4.2. A priori bounds

For the remaining part of this section, we make the assumption that $r_i$ is positive on $X_i$ for $i=1,\dots ,N$, and we define the lower and upper bounds $0<d_i^m\le d_i\left(x\right)\le d_i^M$, $0<r_i^m\le r_i\left(x\right)\le r_i^M$.

Theorem 4.1

Assume that the matrix $\tilde{A}$ defined by ${\tilde{a}}_{ii}:=d_i^m{a}_{ii}$ and ${\tilde{a}}_{ij}:=d_i^M{a}_{ij}$ is Hurwitz. Then the solutions to (3(3) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{n}_i(t,x)=\left(r_i\left(x\right)+d_i\left(x\right)\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right)\right)n_i(t,x),i=1,\dots ,N,$(3) ) are globally defined and bounded.

Proof.

First remark that since $\tilde{A}$ is Hurwitz, then so is A from Lemma 4.2.

We integrate the equations with respect to x and bound them (through $r_i\left(x\right)\le r_i^M$) $\frac{d}{dt}{\rho }_i\left(t\right)\le \left(r_i^M+\sum _{j=1}^N{a}_{ij}{\rho }_j\left(t\right){\int }_{X_i}{d}_i\left(x\right)n_i(t,x)dx\right)i=1,\dots ,N.$ Since the diagonal elements are negative, the off-diagonal non-negative, we obtain $\frac{d}{dt}{\rho }_i\left(t\right)\le \left(r_i^M+a_{ii}{d}_i^m{\rho }_i+\sum _{j\ne i}^N{a}_{ij}{d}_i^M{\rho }_j\left(t\right)\right){\rho }_i\left(t\right),i=1,\dots ,N.$ Thus, the vector $({\rho }_1,\dots ,{\rho }_N)$ is a subsolution of the previous system which is nothing but classical Lotka-Volterra equations with interaction matrix $\tilde{A}$. Thanks to (4.1), the solutions to this system are bounded. Thus, so are those of the integro-differential one.

Remark 4.1

Note that the assumption that $\tilde{A}$ is Hurwitz reduces to A being Hurwitz in the case of constant coefficients. Thus, this result for boundedness is sharp, since it is exactly the one required for convergence to the coexistence steady-state when the equations at hand are classical Lotka-Volterra equations.

Using Theorem 4.1, we can thus define ${\rho }^M>0$ as the GAS steady-state for the system obtained in the previous proof, that is to the equations $\frac{d}{dt}{u}_i=\left(r_i^M+a_{ii}{d}_i^m{u}_i+\sum _{j\ne i}^N{a}_{ij}{d}_i^M{u}_j\left(t\right)\right)u_i\left(t\right),i=1,\dots ,N.$ In other words, ${\rho }^M:=-(\tilde{A}{)}^{-1}{r}^M$ where $r^M$ is the vector $(r_i^M{)}_{1\le i\le N}$. This means that we can write (22) $\underset{t\to +\mathrm{\infty }}{lim sup}{\rho }_i\le {\rho }_i^{M}i=1,\dots ,N.$(22) In a similar fashion to the previous proposition, bounding ${\rho }_i$ away from 0: $\frac{d}{dt}{\rho }_i\ge \left(r_i^m+a_{ii}{d}_i^M{\rho }_i+\sum _{j\ne i}^N{a}_{ij}{d}_i^m{\rho }_j\left(t\right)\right){\rho }_i\left(t\right),i=1,\dots ,N,$ leading to (23) $\underset{t\to +\mathrm{\infty }}{lim inf}{\rho }_i\ge {\rho }_i^{m}i=1,\dots ,N$(23) where ${\rho }^m:=-(B{)}^{-1}{r}^m>0$ with B, a Hurwitz matrix defined by $b_{ii}:=d_i^M{a}_{ii}$, $b_{ij}:=d_i^m{a}_{ij}$ for $i\ne j$ and $r^m:=(r_i^m{)}_{1\le i\le N}$.

4.3. GAS in the mutualistic case

We can now state the main result:

Theorem 4.2

Assume $r_i>0$ for all $i=1,\dots ,N,$ and that the matrix $\hat{A}$ defined by ${\hat{a}}_{ii}:=d_i^m{\rho }_i^m{a}_{ii}$ and ${\hat{a}}_{ij}:=d_i^M{\rho }_i^M{a}_{ij}$ for $i\ne j,$ is Hurwitz.

Then $\tilde{A},$ A and B are also Hurwitz. Furthermore, ${\rho }^{\mathrm{\infty }}:=-A^{-1}{I}^{\mathrm{\infty }}$ lies in ${\mathbb{R}}_{>0}^N$ and it is GAS.

Proof.

The fact that $\tilde{A}$, A and B are also Hurwitz is a direct consequence of Lemma 4.2.

$\tilde{A}$ being Hurwitz, the solutions are globally defined with ρ bounded thanks to Theorem 4.1.

A being Hurwitz, it is invertible and ${\rho }^{\mathrm{\infty }}:=-A^{-1}{I}^{\mathrm{\infty }}$ is in ${\mathbb{R}}_{>0}^N$ thanks to Lemma 4.3.

Let us now prove that it is GAS. The idea is to prove that each ${\rho }_i$ is BV on $[0,+\mathrm{\infty })$. Identifying the limit is straightforward, thanks to the reasoning made in Section 2.

For any i, we define $q_i:=d{\rho }_i/dt$ and write $R_i=R_i(x,{\rho }_1,\dots ,{\rho }_N)$ for readability. Since $q_i=d{\rho }_i/dt={\int }_{X_i}{R}_i{n}_i$, we obtain $\begin{array}{rl}\frac{d{q}_i}{dt}& ={\int }_{X_i}{R}_i^2{n}_i+{\int }_{X_i}\left(\sum _{j=1}^N\frac{\mathrm{\partial }{R}_i}{\mathrm{\partial }{\rho }_j}{q}_j\right)n_i\\ & \ge \sum _{j=1}^N{a}_{ij}\left({\int }_{X_i}{d}_i\left(x\right)n_i(t,x)dx\right)q_j.\end{array}$ Let $b_i\left(t\right):={\int }_{X_i}{d}_i\left(x\right)n_i(t,x)dx$. The idea is that ${\rho }_i$ is ‘mostly’ increasing, so we are interested in the negative part of $q_i$, denoted by $(q_i{)}_{-}$. For this quantity we have the (a.e.) bound $\frac{d(q_i{)}_{-}}{dt}\le b_i\sum _{j=1}^N{a}_{ij}{q}_j\left(-1_{q_i<0}\right)$ On the one hand, $b_i{a}_{ii}{q}_i\left(-1_{q_i<0}\right)=b_i{a}_{ii}(q_i{)}_{-}.$ On the other hand, for $i\ne j$, $b_i{a}_{ij}{q}_j\left(-1_{q_i<0}\right)\le b_i{a}_{ij}(q_j{)}_{-}.$ Combining these two, we get $\frac{d(q_i{)}_{-}}{dt}\le b_i{\left(A(q{)}_{-}\right)}_i.$ We fix $\epsilon >0$ small enough and t large enough (say $t\ge t_0$) such that $\hat{A}+\epsilon J$ is Hurwitz (where J is the matrix composed of ones only) and such that, for each $(i,j)$, $b_i\left(t\right)a_{ij}\le {\hat{a}}_{ij}+\epsilon$. The first requirement is easily derived from Lemma 4.2 since $\hat{A}+\epsilon J$ is clearly cooperative and negatively diagonally dominant for $\epsilon >0$ small enough. The second requirement comes from the lower and upper bounds for the functions ${\rho }_i$ as stated in (22(22) $\underset{t\to +\mathrm{\infty }}{lim sup}{\rho }_i\le {\rho }_i^{M}i=1,\dots ,N.$(22) ) and (23(23) $\underset{t\to +\mathrm{\infty }}{lim inf}{\rho }_i\ge {\rho }_i^{m}i=1,\dots ,N$(23) ).

The resulting inequality is $\frac{d(q_i{)}_{-}}{dt}\le {\left(\left(\hat{A}+\epsilon J\right)(q{)}_{-}\right)}_i,$ so that $\left(\right(q_1{)}_{-},\dots ,\left(q_N{)}_{-}\right)$ is a subsolution of the system with same initial conditions at $t_0$, given by $\frac{dy}{dt}=\left(\hat{A}+\epsilon J\right)y.$ The solutions to this system go exponentially to 0 since $\hat{A}+\epsilon J$ is Hurwitz.

For any i, we have thus proved that $(q_i{)}_{-}$ goes to 0 exponentially. Together with the fact that ${\rho }_i$ is bounded from above, we conclude that it is BV on $[0,+\mathrm{\infty })$. Indeed it holds true that a function $u$ which is both bounded from above and such that $u_{-}\in L^1\left(\right[0,+\mathrm{\infty }\left)\right)$ is BV on $[0,+\mathrm{\infty })$, see [32, Lemma 6.7]. Therefore, ρ converges (to ${\rho }^{\mathrm{\infty }}$).

5. Conclusion

We have analysed the global asymptotic stability properties for integro-differential systems of N species structured by traits x belonging to different trait spaces $X_i$. The coupling comes from a non-local logistic term, which is a linear combination of the total number of individuals ${\rho }_i$ in each species. Theses systems generalize the usual Lotka-Volterra ODEs for which many stability analyses are available in the literature. Our main focus has been on the asymptotic behaviour of the functions ${\rho }_i\left(t\right)$ as $t\to +\mathrm{\infty }$, especially towards equilibrium ${\rho }^{\mathrm{\infty }}$ with positive components, i.e. of persistence of all species. In Section 2, we explained how identifying it essentially determines the asymptotic behaviour of the underlying density $n_i$, namely the phenotypes on which the measures $n_i(t,\cdot )$ concentrate in large time.

In Section 3, an adequate Lyapunov functional allowed us to state a general result relying solely on an assumption on the matrix A, regardless of the type of interactions. For N=2, this is essentially a sharp result, but becomes more restrictive for $N\ge 3$. This tool also provided us with convergence rates to equilibrium. In Section 4, we presented another strategy based on a BV bound which yielded a second result of global asymptotic stability, this time for mutualistic equations.

The result of Theorem 4.2 is partly less general than the one given in Theorem 3.1 because it requires a sign on the coefficients of A. However, the set of matrices which satisfy the hypothesis given in the last theorem is an open subset of the set of real matrices ${\mathbb{R}}^{N\times N}$ in any dimension. This in sharp contrast with the hypothesis of Theorem 3.1, which, as already mentioned, imposes some polynomial equalities on the coefficients of A as soon as $N\ge 3$. In other words, for a small perturbation of a cooperative matrix for which GAS holds, GAS still holds. In particular, if one has weakly (but mutualistically) coupled equations, GAS holds, whereas Theorem 3.1 does not cover the case of any weakly coupled equations for general interactions, except for N=2.

In both cases, the assumptions fall within the class of matrices which cannot have off-diagonal coefficients which are too high compared to the diagonal ones. The present results thus apply to cases where interactions among individuals of a same species are not only blind because of the term $a_{11}{\rho }_1$, but also stronger than the interactions between species. In other words, each one of them has its own ecological niche inside which interactions are independent of how given phenotypes x and y are away from another.

Let us remark that the BV method would apply to more general functions $R_i(x,{\rho }_1,\dots ,{\rho }_N)$, as long as they are increasing in the variables ${\rho }_j$, $j\ne i$. However, the Lyapunov functional used in Theorem 3.1 seems to be dependent on the linear coupling chosen here and it is an open problem to generalize our results for other settings. Another open question is about finding whether there are matrices A for which the underlying classical Lotka-Volterra equations converge to the coexistence steady state (for example such that there exists D>0 with $A^{T}D+DA<0$), but for which there is no GAS for the integro-differential system. Numerically at least, we could not build any such case.

We finally mention a natural extension of these integro-differential systems, which has drawn much attention in adaptive dynamics: that is when one adds a mutation term $\mathcal{M}\left[n\right]$, and a typical single equation then writes (24) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}n(t,x)=\left(r\left(x\right)-{\int }_{X}K(x,y)n(t,y)dy\right)n(t,x)+\mathcal{M}\left[n\right](t,x).$(24) This is usually done either through a second order elliptic operator like $\mathcal{M}\left[n\right](t,x)=\beta \mathrm{\Delta }n(t,x),$ with Neumann boundary conditions, or using an integral operator allowing for long-range mutations like $\mathcal{M}\left[n\right](t,x)={\int }_{X}m(x,y)\left(n(t,y)-n(t,x)\right)dy.$ Finally, let us mention that in some more recent works an advection term is also considered [10,11]. The idea is to model stress-induced adaptation: individuals actively adapt to their environment and this can be thought of as an appropriate modelling of Lamarckism induction.

Most of the studies on model (24(24) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}n(t,x)=\left(r\left(x\right)-{\int }_{X}K(x,y)n(t,y)dy\right)n(t,x)+\mathcal{M}\left[n\right](t,x).$(24) ) have aimed at understanding how small mutations affect the dynamics of the surviving phenotype (when one expects a single Dirac mass) with a vanishing mutation term ${\mathcal{M}}_{\epsilon }$, after a proper rescaling of time [2,27,30]. Without smallness assumptions on the mutations, non-trivial steady-states have been investigated in detail in [3,12,25], and results of GAS have essentially been obtained for $K(x,y)=k\left(y\right)$. It would be interesting to investigate the asymptotic stability properties of steady states for more general kernels. To the best of our knowledge, general asymptotic results are indeed still elusive, even in the case explored in our paper, namely when $K(x,y)=d\left(x\right)$ for d non constant. It is not clear whether the techniques developed in the present paper can be extended to that setting nor to systems of this form.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Camille Pouchol  http://orcid.org/0000-0002-1152-9896