# Continuous limits of large plant-pollinator random networks and some
applications††thanks: The authors thank Mathilde Dufaÿ and Pauline Lafitte-
Godillon who participated to the early discussions on this project. They also
thank Emmanuel Grenier for his help in the proof of Proposition 4.5. S.B.,
H.L. and V.C.T. acknowledge support of the “Chaire Modélisation Mathématique
et Biodiversité (MMB)” of Veolia Environnement-Ecole Polytechnique-Museum
National d’Histoire Naturelle-Fondation X. T.R. is partially funded by Labex
CEMPI (ANR-11-LABX-0007-01) and ANR Project MoHyCon (ANR-17-CE40-0027-01).
V.C.T is supported by ANR EcoNet (ANR-18-CE02-0010), ANR CADENCE
(ANR-16-CE32-0007) and Labex Bézout (ANR-10-LABX-0058).

Sylvain Billiard111Univ. Lille, CNRS, UMR 8198 - Evo-Eco-Paleo, F-59000 Lille,
France<EMAIL_ADDRESS>Hélène Leman222Inria, ENSL, UMPA, CNRS
UMR 5669, 69364 Lyon, France<EMAIL_ADDRESS>Thomas Rey333Univ. Lille,
CNRS, Inria, UMR 8524 - Laboratoire P. Painlevé, F-59000 Lille, France.
<EMAIL_ADDRESS>and Viet Chi Tran444LAMA, Univ Gustave Eiffel, Univ
Paris Est Creteil, CNRS, F-77454 Marne-la-Vallée, France. chi.tran@univ-
eiffel.fr

###### Abstract

We study a stochastic individual-based model of interacting plant and
pollinator species through a bipartite graph: each species is a node of the
graph, an edge representing interactions between a pair of species. The
dynamics of the system depends on the between- and within-species
interactions: pollination by insects increases plant reproduction rate but has
a cost which can increase plant death rate, depending on pollinators density.
Pollinators reproduction is increased by the resources harvested on plants.
Each species is characterized by a trait corresponding to its degree of
generalism. This trait determines the structure of the interactions graph and
the quantity of resources exchanged between species. Our model includes in
particular nested or modular networks. Deterministic approximations of the
stochastic measure-valued process by systems of ordinary differential
equations or integro-differential equations are established and studied, when
the population is large or when the graph is dense and can be replaced with a
graphon. The long-time behaviors of these limits are studied and central limit
theorems are established to quantify the difference between the discrete
stochastic individual-based model and the deterministic approximations.
Finally, studying the continuous limits of the interaction network and the
resulting PDEs, we show that nested plant-pollinator communities are expected
to collapse towards a coexistence between a single pair of species of plants
and pollinators.

This work is dedicated to S. Méléard and has been presented at her 60th
birthday’s conference (2018).

Keywords: ecological mutualistic community, birth and death process,
interacting particles, limit theorem, kinetic limit, graphon, integro-
differential equation, stationary solution.
AMS Codes: 92D40, 92D25, 05C90, 60J80, 60F17, 47G20.

Communities of interacting species, such as plants and pollinators, are
systems characterized by a large diversity of both their components (dozens of
species are typically involved) and the topology of their interaction networks
[18]. How does the structure of the interaction networks affect the dynamics
of species and the stability of communities is a common and long-standing
question in the ecological literature, e.g. [40, 46, 38]. Addressing this
question is challenging because the number of species in a community is large,
and so is the number of possible edges. Describing ecological networks from
observations thus needs much effort, without being ensured that the inferred
networks are representative of the studied communities (interactions can
indeed depend on space and time, but also on the abundances of the species
which makes them harder or not to observe). From a theoretical point of view,
the study of the dynamics and stability of communities most often relies
either on the numerical analysis of large systems of ordinary differential
equations (ODEs), or on the study of the system near equilibrium in simplified
communities (see e.g. [40, 41, 7, 14, 17, 26, 44, 45, 46, 38] and [2] for a
review). This makes difficult the identification of general properties of
ecological communities. Overall, despite being more than 50 years old, the
question of how correctly modelling ecological networks and communities is
still challenging. In this paper, we develop a new theoretical framework that
could help in addressing the previous questions, based on the fact that
networks are hierarchically structured. Interactions can indeed be considered
at different scales: between individuals, between species or at the level of
the whole community [29]. At all scales however, several studies showed
evidence that the structure of the interaction networks is due to traits
variation between or within species’ such as size, symmetry or time of
activity, e.g. [15, 49, 47, 4]. Starting from a stochastic individual-based
model, our goal is to provide simplifications of an ecological network
deriving the continuous limits of both the population sizes and the
interaction graph structured by a trait.

We focus on a particular ecological community: plant and pollinator species.
The interactions between individuals of each species are modelled by a
bipartite random network. Interactions affect the demographic rates of plant
and pollinator individuals because of resources exchanges. On the one hand,
plants benefit from pollinators visitations to increase their reproduction
rate. On the other hand, pollinators consume nectar, leaves, pollen, etc.
which increases their reproduction or survival rates. However, producing such
resources is costly for the plants, which can negatively affect plant
demographic rates, proportionally to pollinators density [30, 32].

The interaction rates between plant and pollinator individuals are modelled by
a random graph where the species are the nodes of the graph. The existence of
an edge between a plant and a pollinator species indicates that individuals of
these two species can interact. The resulting graph is bipartite, since there
is by definition no direct edges between insect species or between plant
species. The topology of the graph depends on a trait that represents the
degree of generalism of the species. Any pair of plant-pollinator species is
connected independently from the other pairs with a probability that depends
on the traits of the two considered species. This simple model corresponds to
Erdös-Rényi graphs, where the probability of connection is the same for every
pairs. Such a model can also generate particular structures such as nested or
modular graphs which are commonly found in ecological networks [46, 29].

In Section 1, we present the stochastic individual based model, how the random
interactions graph is modelled, and how interactions affect demographic rates.
Considering large populations but with a fixed number of species, we show in
Section 2 how the dynamics can be approximated by a system of ODEs (close to
the ones commonly studied in the ecological literature). The fluctuations
between the approximated deterministic limit and the stochastic individual
based process are established.
In Section 3, we derive a continuous approximation of the random graph when
the numbers of plant and pollinator species, say $n$ and $m$, are also large.
When the graph is dense (i.e. the order of the number of edges is in
$O(n\times m)$), the complex random network can be replaced by a continuous
object, namely a graphon [13, 39]. In this case, the high-dimensional system
of ODEs can be replaced by two partial integro-differential equations: one for
the plants and the other for the pollinators. The latter equations fall into
the broader class of kinetic equations, the most well known and studied being
the seminal Boltzmann and Vlasov equations. The use of such mesoscopic scale,
namely intermediate between a microscopic agent based approach [25], and an
averaged macroscopic one (over space) such as [42], is not a novelty for
modeling competitive interactions between species [20, 35, 9, 10], i.e. with a
large but not infinite number of species with trait structure. Indeed, such
approach can be traced back as early as the book of [43] for the mathematical
modeling of evolutionary ecology. However, the kinetic treatment of large
random networks and their approximations by integro-differential equations is
new to our knowledge. Finally, in Section 4, the large-time behaviors of the
ODE system or of the integro-differential equations are studied and
simulations are produced to explore several situations.

##  1 Stochastic individual based model with species interactions

In this section, we introduce a stochastic individual-based model of plant and
pollinator species. In this model, the number of species and the number of
individuals in each species are finite. Each species is defined by a trait,
for instance a morphological or functional trait, which determines its degree
of generalism. This trait shapes the interaction network: generalist plant
species can be visited by a large number of insect species and generalist
pollinator species can visit a large number of plant species. Even though
several coevolved traits can be considered in the model (e.g. orchids have
flower with a particular morphology, color, phenology, etc.), we will focus on
a single trait for the sake of simplicity. The dynamics is driven by the
births and deaths of individuals at random times that depend 1) on their
position in the plant-pollinator network, determined by their species’ trait,
2) by the weights of the interactions between pairs of species, and 3) by the
population sizes of the species. The evolution of these interacting
populations is described by a stochastic differential equation (SDE) involving
Poisson point measures and an acceptance-rejection algorithm as in the
Gillespie algorithm [27] (see Fournier and Méléard [25] or Ferrière and Tran
[23] for the mathematical formulation of the Gillespie algorithm by mean of
SDEs in cases with interactions).

###  1.1 Description of the plant-pollinator community

We consider $n$ and $m$ plant and pollinator species, respectively. Each plant
specie is characterized by a trait $x\in[0,1]$ and each pollinator specie by a
trait $y\in[0,1]$. For
$i\in\\{1,\dots,n\\}=\left[\\!\left[1,n\right]\\!\right]$ and
$j\in\\{1,\dots,m\\}=\left[\\!\left[1,m\right]\\!\right]$, we respectively
denote $x^{i}$ and $y^{j}$ the trait of the plant species $i$ and of the
pollinator species $j$. These traits can represent for instance their degree
of generalism, i.e. their tendency to interact with a large number of other
species: a species is considered as generalist when its trait ($x^{i}$ or
$y^{j}$) is close to $1$, or specialist when it is close to $0$.

###  1.2 The plant-pollinator interaction network as a bipartite random graph

Plants and pollinators interact through a bipartite network where each species
is a vertex, and an interaction is an edge. We denote $i\sim j$ or $j\sim i$
when individuals of the plant species $i$ can interact with individuals of the
pollinator species $j$. There are no edges between two plant or two pollinator
species because the network only represents interactions between plants and
pollinators. Yet, when specifying the birth and death rates, competition
kernels within and between plant species, and within and between pollinator
species, will be introduced in the model.

The bipartite graph can be represented by an $n\times m$ adjacency matrix
$G^{n,m}$, with

$G^{n,m}_{ij}=\left\\{\begin{aligned} 1&\text{ if the pollinator species $j$
can interact with the plant species $i$},\\\ 0&\text{
otherwise}.\end{aligned}\right.$

Regarding the structure of the network, we consider a stochastic bipartite
graph that generalizes the Erdös-Rényi random graph ([12, 21, 19]). The
probability $\phi(x^{i},y^{j})$ that there is an edge $i\sim j$ between plant
species $i$ and pollinator species $j$ is assumed to depend only on the traits
$x^{i}$ and $y^{j}$ of these two species. Different pairs are supposed to be
connected independently, the interaction graph is thus built such that:

######  Assumption 1.1.

We assume that there exists a continuous function $\phi:[0,1]^{2}\to[0,1]$,
such that conditionally on
$(x^{i},y^{j})_{i,j\in\left[\\!\left[1,n\right]\\!\right]\times\left[\\!\left[1,m\right]\\!\right]}$,
the random variables
$(G^{n,m}_{ij})_{(i,j)\in\left[\\!\left[1,n\right]\\!\right]\times\left[\\!\left[1,m\right]\\!\right]}$
are independent and distributed as Bernoulli random variables with parameters
$(\phi(x^{i},y^{j}))_{(i,j)\in\left[\\!\left[1,n\right]\\!\right]\times\left[\\!\left[1,m\right]\\!\right]}$.

|
---|---
(a) | (b)
Figure 1: (a): Nested (left) or modular (right) bipartite networks, from
Fontaine et al. [24]. (b) Pollination network for diurnal and nocturnal insect
species, from Knop et al. [37]

######  Example 1.2.

Case 1. If there exists $\phi_{0}\in(0,1)$ such that for all $i$ and $j$,
$\phi(x^{i},y^{j})=\phi_{0}$, then the stochastic network is a bipartite
Erdös-Rényi graph. The degrees of nodes representing pollinators (resp.
plants), i.e. the numbers of edges, do not depend on the species traits.
Degrees are thus independent binomial r.v. with parameters $(m,\phi_{0})$
(resp. $(n,\phi_{0})$) and the total number of edges follows a binomial
distribution with parameters $(nm,\phi_{0})$.
Case 2. The probability of an edge $i\sim j$ between two species $i$ and $j$
depends on their degree of generalism, for instance by assuming:

$\phi(x^{i},y^{j})=x^{i}y^{j}.$ (1.1)

Under this assumption, two generalist species ($x^{i}$ and $y^{j}$ close to
$1$) have a higher probability to be connected than two specialist species
($x^{i}$ and $y^{j}$ close to $0$). Nested graphs result from such an
assumption (see Fig. 1 a-left for an illustration of nested graphs
definition).
Case 3. When plant and pollinator species can be structured in groups, for
instance because of pollination syndrome or spatio-temporal segregation, the
probability of an edge $i\sim j$ $\phi(x^{i},y^{j})=\phi_{IJ}$ depends only on
classes $I\ni i$ and $J\ni j$. The resulting random graph is a Stochastic
Block Model (SBM [33], see e.g. [1] for a review), often called modular
networks in community ecology (see Fig. 1 a-right for an illustration of
modular graphs definition).

If species $i$ and $j$ interact, i.e. if $G^{n,m}_{ij}=1$, we denote by
$c_{ij}^{n,m}$ the weight of the interaction $i\sim j$. The quantity
$c_{ij}^{n,m}$ describes the intensity and frequency of the relation between
the pollinator species $j$ and the plant species $i$. From the point of view
of the plant, $c_{ij}^{n,m}$ can be interpreted as a measure of the
pollination services received from the pollinators. From the point of view of
the pollinators, $c_{ij}^{n,m}$ measures the quantity and quality of nutrients
collected from the plants.

######  Assumption 1.3.

For all $n,m\in\mathbb{N}^{2}$, conditionally on
$(x^{i},y^{j})_{(i,j)\in\left[\\!\left[1,n\right]\\!\right]\times\left[\\!\left[1,m\right]\\!\right]}$,
$(c_{ij}^{n,m})_{(i,j)\in\left[\\!\left[1,n\right]\\!\right]\times\left[\\!\left[1,m\right]\\!\right]}$
is assumed to be a sequence of independent random variables such that

  * •

the expected values only depend on the traits of the plants and pollinators
through a function $c^{n,m}\ :[0,1]^{2}\rightarrow\mathbb{R}$: for all
$i\in\left[\\!\left[1,n\right]\\!\right],j\in\left[\\!\left[1,m\right]\\!\right]$:

$\mathbb{E}_{x,y}\left[c^{n,m}_{ij}\right]=c^{n,m}\big{(}x^{i},y^{j}\big{)},$
(1.2)

with
$\mathbb{E}_{x,y}[A]=\mathbb{E}[A|(x^{i},y^{j})_{(i,j)\in\left[\\!\left[1,n\right]\\!\right]\times\left[\\!\left[1,m\right]\\!\right]}]$,

  * •

and the variances decrease with $n$ and $m$ such that:

$V_{max}:=\sup_{n\in\mathbb{N},m\in\mathbb{N}}\sup_{(i,j)\in\left[\\!\left[1,n\right]\\!\right]\times\left[\\!\left[1,m\right]\\!\right]}\text{Var}_{x,y}((m+n)c^{n,m}_{i,j})<\infty.$
(1.3)

######  Example 1.4.

Depending on plants and pollinators’ traits considered, different function
$c^{n,m}_{i,j}$ be considered:

  * •

The intensity and frequency of the interaction can be independent of the
species traits, i.e. $c^{n,m}_{i,j}=\frac{c}{n+m}$, which would mean that
given a plant and a pollinator species interact, the effect of the interaction
on their demographic rates is the same for all species;

  * •

$x^{i}$ and $y^{j}$ can represent the time of peak activity during a day or a
season for plants $i$ and pollinators $j$, and $c^{n,m}_{i,j}$ a function
measuring the overlap of activity for both species; Similarly, $c^{n,m}_{i,j}$
can represent geographical rather than time overlap between both species;

  * •

$x^{i}$ and $y^{j}$ can virtually represent some aspects of the morphology of
plants’ flowers and pollinators’ body and the function $c^{n,m}_{i,j}$ some
morphological preference; For instance, assuming $(x^{i},y^{j})\in[0,1]^{2}$,
$c^{n,m}_{i,j}=\frac{x^{i}y^{j}}{n+m}$ would mean that pollinators and plants
tend to interact more when the traits are both large, while if
$c^{n,m}_{i,j}=\frac{x^{i}(1-y^{j})}{n+m}$ would mean that they will interact
more when the plant has a large and the pollinator a low trait values.

###  1.3 Stochastic dynamics of the plant-pollinator community

The dynamics of the plant and pollinator populations within the community is
ruled by random point birth and death events. We assume that the size of the
population of plants and pollinators is scaled by a factor $K>0$, called the
carrying capacity in ecology. This carrying capacity $K$ is a measure of the
size of the system, in other words, it controls the total abundance of the
whole community that can be sustained by the environment. A continuous limit
of the stochastic dynamics will be obtained when the species abundances tend
to infinity, in other words when $K\rightarrow\infty$.

Given the scaling factor $K$, we denote $P^{K,i}_{t}$ and $A^{K,j}_{t}$ the
size of the plant and pollinator species $i$ and $j$ at time $t$. The plant
and pollinator populations at time $t$ can be represented by the following
point measures:

$\mathbf{P}^{K,n,m}_{t}(dx)=\frac{1}{nK}\sum_{i=1}^{n}P^{K,i}_{t}\delta_{x^{i}}(dx),\qquad\mathbf{A}^{K,n,m}_{t}(dy)=\frac{1}{mK}\sum_{j=1}^{m}A^{K,j}_{t}\delta_{y^{j}}(dy).$
(1.4)

The interactions between plants and pollinators are characterized by an
exchange of resources [32]. The quantity of resources gained by plants or
pollinators $R$ is modeled through a "mass-action model". At time $t$, in the
population scaled by $K$, a single individual of the plant species $i$
interacting with pollinator species $j$ is supposed to gain a quantity of
resources (here the pollination service) proportional to the abundance of
pollinators $A^{K,j}_{t}/K$ weighted by the interaction efficiency
$c_{ij}^{n,m}$, such that the total resource gained by a plant individual of
species $i$ through the pollination interactions is

$R^{A,K,i}_{t}:=\sum_{j\sim i}\frac{c^{n,m}_{ij}A^{K,j}_{t}}{K}.$ (1.5)

Similarly, for a given pollinator of the species $j$, the resources gained
from the interaction with the plant species at time $t$ in the population
parameterized by $K$ is assumed to be:

$R^{P,K,j}_{t}=\sum_{i\sim j}\frac{c^{n,m}_{ij}P^{K,i}_{t}}{K}.$ (1.6)

The dynamics of the community is supposed to be governed by the birth and
death rates of the plant and pollinator populations. We denote $b^{P}(R)$ and
$b^{A}(R)$ the individual birth rate of plant and pollinator species,
respectively, each depending on the quantity of resources exchanged $R^{P}$ or
$R^{A}$. Similarly, we denote $d^{P}(R)$ and $d^{A}(R)$ the individual death
rates. Finally,

$g^{P}(R):=b^{P}(R)-d^{P}(R)\quad\text{and}\quad g^{A}(R):=b^{A}(R)-d^{A}(R),$
(1.7)

are the component of the growth rate due to the interactions between plants
and pollinators.

The plants and pollinators dynamics are also assumed to be affected by
logistic competition among plants and among pollinators (within and between
species competition). We suppose that competition strength depends on the
traits $x$ and $y$. This can represent the fact that plants or pollinators
with similar traits tend to share similar ecological niche, or phenology, etc.
A plant with trait $x\in[0,1]$ suffers an additional death rate term due to
competition such that

$k\star\mathbf{P}^{K,n,m}_{t}(x):=\int_{[0,1]}k(x,x^{\prime})d\mathbf{P}^{K,n,m}_{t}(dx^{\prime})=\frac{1}{nK}\sum_{i=1}^{n}k(x,x^{i})P^{K,i}_{t},$
(1.8)

where $k(x,x^{\prime})$ quantifies the competition pressure exerted by another
plant of trait $x^{\prime}$. Similarly, a pollinator with trait
$y^{\prime}\in[0,1]$ suffers an additional death rate due to competition with
pollinators of trait $y$ is

$h\star\mathbf{A}^{K,n,m}_{t}(y):=\int_{[0,1]}h(y,y^{\prime})d\mathbf{A}^{K,n,m}_{t}(dy^{\prime})=\frac{1}{mK}\sum_{j=1}^{m}h(y,y^{j})A^{K,j}_{t},$
(1.9)

where $h(y,y^{\prime})$ quantifies the competition of individuals with trait
$y^{\prime}$ on individuals with trait $y$.

The following assumptions on the functions involved in the model will be
needed, both for modeling and mathematical purposes.

######  Assumption 1.5.

(i) The birth rates $b^{P}$ and $b^{A}$ are assumed to be bounded on
$\mathbb{R}^{+}$ by constants $M^{P}>0$ and $M^{A}>0$ respectively. Moreover,
all the rate functions $b^{P}$, $b^{A}$, $d^{P}$ and $d^{A}$ are assumed to be
locally Lipschitz continuous on $[0,\infty)$.
(ii) The competition kernels $k$ and $h$ are assumed to be continuous on
$[0,1]^{2}$.

###  1.4 Stochastic differential equations

Following works by Méléard and co-authors [25, 16], it is possible to describe
the evolution of the population measures
$(\mathbf{P}^{K,n,m}_{t},\mathbf{A}^{K,n,m}_{t})_{t\in\mathbb{R}_{+}}$ defined
in (1.4) by stochastic differential equations (SDEs) driven by Poisson point
measures. Let us first present these SDEs and then explain their heuristic
meaning.

######  Definition 1.6.

Suppose that assumptions 1.1, 1.3 and 1.5 hold. Let
$Q_{B}^{P}(ds,dk,d\theta)$, $Q_{D}^{P}(ds,dk,d\theta)$, (resp.
$Q_{B}^{A}(ds,dk,d\theta)$ and $Q^{A}_{D}(ds,dk,d\theta)$) be Poisson point
measures on $\mathbb{R}_{+}\times
E:=\mathbb{R}_{+}\times\left[\\!\left[1,n\right]\\!\right]\times\mathbb{R}_{+}$
(resp. $\mathbb{R}_{+}\times
F:=\mathbb{R}_{+}\times\left[\\!\left[1,m\right]\\!\right]\times\mathbb{R}_{+}$)
with intensity measure
$q(ds,dk,d\theta)=ds\,{\color[rgb]{0,0,0}\mathfrak{n}}(dk)\,d\theta$ where
$ds$ and $d\theta$ are Lebesgue measures on $\mathbb{R}_{+}$ and where
${\color[rgb]{0,0,0}\mathfrak{n}}(dk)$ is the counting measure on
$\mathbb{N}^{*}=\\{1,2,\dots\\}$. For $i\in\\{1,\dots n\\}$ and
$j\in\\{1,\dots m\\}$, one has

$\displaystyle P^{K,i}_{t}=$ $\displaystyle
P^{K,i}_{0}+\int_{0}^{t}\int_{E}{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}_{i=k}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\theta\leq
b^{P}\big{(}R^{A,K,i}_{s-}\big{)}P^{K,i}_{s-}}Q^{P}_{B}(ds,dk,d\theta)$
$\displaystyle\hskip 28.45274pt-\int_{0}^{t}\int_{E}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{i=k}{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\theta\leq\big{[}d^{P}\big{(}R^{A,K,i}_{s-}\big{)}+k\star\mathbf{P}^{K,n}_{s-}(x^{i})\big{]}P^{K,i}_{s-}}Q^{P}_{D}(ds,dk,d\theta)$
(1.10) $\displaystyle A^{K,j}_{t}=$ $\displaystyle
A^{K,j}_{0}+\int_{0}^{t}\int_{F}{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}_{j=k}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\theta\leq
b^{A}\big{(}R^{P,K,j}_{s-}\big{)}A^{K,j}_{s-}}Q^{A}_{B}(ds,dk,d\theta)$
$\displaystyle\hskip 28.45274pt-\int_{0}^{t}\int_{F}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{j=k}{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\theta\leq\big{[}d^{A}\big{(}R^{P,K,j}_{s-}\big{)}+h\star\mathbf{A}^{K,m}_{s-}(y^{j})\big{]}A^{K,j}_{s-}}Q^{A}_{D}(ds,dk,d\theta).$
(1.11)

The above SDEs correspond to the mathematical formulation of the individual-
based simulations classically used in ecology and originated with Gillespie’s
algorithm [27, 28]. The Poisson point measures $Q_{B}^{P}$ and $Q_{B}^{A}$
(resp. $Q_{D}^{P}$ and $Q_{D}^{A}$) give the random times of possible birth
(resp. death) events. The indicators ${\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}_{i=k}$ or
${\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{j=k}$ in the integral ensure that the equations for
$P^{K,i}$ or $A^{K,j}$ are indeed modified when the birth or death events
affect the plant population $i$ or pollinator population $j$. The indicators
in $\theta$ correspond to an acceptance-rejection algorithm so that the birth
and death times occur with the correct rates. Note that in these SDEs, the
network is hidden in the definitions of the terms $R^{A,K,i}_{s}$ and
$R^{P,K,j}_{s}$ (see (1.5) and (1.6)).

In order to approximate the individual-based model with ODEs or integro-
differential equations, we can reformulate the SDEs (1.10) and (1.11) by
expressing the processes $P^{K,i}$ and $A^{K,j}$ as semi-martingales. The
following proposition also states the existence and uniqueness of a solution
to (1.10)-(1.11).

######  Proposition 1.7.

Consider here $n$, $m$ and $K$ fixed. Assume 1.5 and that

$\mathbb{E}\left[\langle\mathbf{P}^{K,n,m}_{0},\mathds{1}\rangle^{2}+\langle\mathbf{A}^{K,n,m}_{0},\mathds{1}\rangle^{2}\right]<+\infty,$

then the processes defined by (1.10) and (1.11) are well defined on
$\mathbb{R}_{+}$. Moreover, for all $T\geq 0$,

$\mathbb{E}\left[\sup_{t\in[0,T]}\langle\mathbf{P}^{K,n,m}_{t},\mathds{1}\rangle^{2}+\langle\mathbf{A}^{K,n,m}_{t},\mathds{1}\rangle^{2}\right]<+\infty.$
(1.12)

For any $f:[0,1]\rightarrow\mathbb{R}$ measurable test function, we have:

$\displaystyle\langle\mathbf{P}^{K,n,m}_{t},f\rangle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{1}{K}P^{K,i}_{t}f(x^{i})$
$\displaystyle=$
$\displaystyle\langle\mathbf{P}^{K,n,m}_{0},f\rangle+\int_{0}^{t}\frac{1}{nK}\sum_{i=1}^{n}f(x^{i})\Big{[}g^{P}\big{(}R^{A,K,i}_{s_{-}}\big{)}-k\star\mathbf{P}^{K,n,m}_{s}(x^{i})\Big{]}P^{K,i}_{s}\
ds$ $\displaystyle+\frac{1}{nK}\sum_{i=1}^{n}f(x^{i})M^{K,i}_{t},$ (1.13)

where $(M^{K,i})_{i\in\\{1,..,n\\}}$ are square integrable martingales with
predictable quadratic variation processes:

$\langle
M^{K,i}\rangle_{t}=\int_{0}^{t}\left(b^{P}\big{(}R^{A,K,i}_{s}\big{)}+d^{P}\big{(}R^{A,K,i}_{s}\big{)}+k\star\mathbf{P}^{K,n,m}_{s}(x^{i})\right)P^{K,i}_{s}\
ds.$ (1.14)

A similar expression holds for the pollinator populations. For any
$f:[0,1]\rightarrow\mathbb{R}$ measurable test function, we have:

$\displaystyle\langle\mathbf{A}^{K,n,m}_{t},f\rangle=$
$\displaystyle\frac{1}{m}\sum_{j=1}^{m}\frac{1}{K}A^{K,j}_{t}f(y^{j})$
$\displaystyle=$
$\displaystyle\langle\mathbf{A}^{K,n,m}_{0},f\rangle+\int_{0}^{t}\frac{1}{mK}\sum_{j=1}^{m}f(y^{j})\Big{[}g^{A}\big{(}R^{P,K,j}_{s_{-}}\big{)}-h\star\mathbf{A}^{K,n,m}_{s}(y^{j})\Big{]}A^{K,j}_{s}\
ds$ $\displaystyle+\frac{1}{mK}\sum_{j=1}^{m}f(y^{j})M^{K,j}_{t},$ (1.15)

where $(M^{K,j})_{j\in\\{1,..,m\\}}$ are square integrable martingales with
predictable quadratic variation processes:

$\langle
M^{K,j}\rangle_{t}=\int_{0}^{t}\left(b^{A}\big{(}R^{P,K,j}_{s}\big{)}+d^{A}\big{(}R^{P,K,j}_{s}\big{)}+h\star\mathbf{A}^{K,n,m}_{s}(y^{j})\right)A^{K,j}_{s}\
ds.$

The proof follows from usual stochastic calculus with Poisson point processes
(e.g. [34]), as developed in [6, 25] for example. We however give a sketch of
proof in Appendix A.1.

Thus the processes $(\mathbf{P}^{K,n,m}_{t})_{t\in\mathbb{R}_{+}}$ and
$(\mathbf{A}^{K,n,m}_{t})_{t\in\mathbb{R}_{+}}$ are well-defined in the set
$\mathbb{D}(\mathbb{R}_{+},\mathcal{M}_{F}([0,1]))$ of right-continuous left-
limited (càdlàg) processes with values in the set of finite measures on
$[0,1]$. The space $\mathcal{M}_{F}([0,1])$ embedded with the weak topology is
a Polish space and the set of càdlàg functions is embedded with the Skorokhod
topology which makes it Polish as well (e.g. [11, Th. 6.8 and 12.2]).

##  2 Community dynamics limit when abundances are large but the number of
species is fixed

In this section, we consider the numbers of species $n$ and $m$ as fixed while
the plant and pollinator populations size, represented by $K$, tend to
$+\infty$. In this section, since $n$ and $m$ are fixed, these notations will
be omitted in order to avoid cumbersomeness when there is no possible
confusion.

###  2.1 Law of large numbers

######  Proposition 2.1.

We consider a sequence
$(\mathbf{P}^{K,n,m},\mathbf{A}^{K,n,m})_{K\in\mathbb{N}}$ of processes as in
Definition 1.6, with initial conditions such that

$\sup_{K\in\mathbb{N}}\mathbb{E}\left[\langle\mathbf{P}^{K,n,m}_{0},\mathds{1}\rangle^{3}+\langle\mathbf{A}^{K,n,m}_{0},\mathds{1}\rangle^{3}\right]<+\infty,$
(2.1)

and such that there exist
$(\widetilde{P}^{1}_{0},\dots,\widetilde{P}^{n}_{0})$ and
$(\widetilde{A}_{0}^{1},\dots\widetilde{A}_{0}^{m})$ satisfying the following
convergences almost surely:

$\lim_{K\rightarrow+\infty}\frac{P^{K,i}_{0}}{K}=\widetilde{P}^{i}_{0},\qquad\lim_{K\rightarrow+\infty}\frac{A^{K,j}_{0}}{K}=\widetilde{A}^{j}_{0}.$

Then, for all $T\geq 0$, the following convergence holds almost surely for all
$i\in\left[\\!\left[1,n\right]\\!\right]$ and
$j\in\left[\\!\left[1,m\right]\\!\right]$:

$\lim_{K\to\infty}\sup_{t\leq
T}\sup_{i,j}\left\\{\left|\frac{P^{K,i}_{t}}{K}-\widetilde{P}^{i}_{t}\right|,\left|\frac{A^{K,j}_{t}}{K}-\widetilde{A}^{j}_{t}\right|\right\\}=0,$
(2.2)

where $(\widetilde{P}^{1}_{t},\dots,\widetilde{P}^{n}_{t})_{t\geq 0}$ and
$(\widetilde{A}_{t}^{1},\dots\widetilde{A}_{t}^{m})_{t\geq 0}$ are the unique
solution of the system, which exists on $\mathbb{R}^{n+m}$:

$\displaystyle\forall 1\leq i\leq n,$
$\displaystyle\frac{d\widetilde{P}^{i}_{t}}{dt}=\left(g^{P}\big{(}\sum_{j\sim
i}\widetilde{R}^{A,ij}_{t}\big{)}-\frac{1}{n}\sum_{\ell=1}^{n}k(x^{i},x^{\ell}){\widetilde{P}}^{\ell}_{t}\right)\widetilde{P}^{i}_{t}$
(2.3) $\displaystyle\forall 1\leq j\leq m,$
$\displaystyle\frac{d\widetilde{A}^{j}_{t}}{dt}=\left(g^{A}\big{(}\sum_{i\sim
j}\widetilde{R}^{P,ij}_{t}\big{)}-\frac{1}{m}\sum_{\ell=1}^{m}h(y^{j},y^{\ell})\widetilde{A}^{\ell}_{t}\
\right)\widetilde{A}^{j}_{t},$

where, for all $t\in\mathbb{R}^{+}$,
$\widetilde{R}^{P,ij}_{t}=c_{ij}^{n,m}\widetilde{P}^{i}_{t}$ and
$\widetilde{R}^{A,ij}_{t}=c_{ij}^{n,m}\widetilde{A}^{j}_{t}$.

Eq. (2.3) is similar to a classical Lotka-Volterra system applied to
mutualistic interactions with competition. In our case, the species community
is structured by the plants traits $x^{i}$ and the pollinators traits $y^{j}$
which determine the probability and strength of the interactions. In addition,
how interactions translate into births and deaths (the so-called numerical
response in ecological terms) is embedded in functions $g^{A}$ and $g^{P}$
(demographic growth) and $k$ and $h$ (competitive kernels). Those functions
can take any form (see section 4 for some examples). As a consequence, Eqs.
(2.3) can capture a large variety of ecological situations. In particular,
many ODE models published in the ecological literature are special cases of
Eqs. (2.3) (e.g. [46, 8, 38]). Our analysis thus shows that ODE models can be
commonly seen as a limits of stochastic individual-based models as defined in
Section 1. Notice that coefficients $c_{ij}^{n,m}$ are random coefficients as
explained in the previous section. However, these coefficients remain
identical with time.

###### Proof.

Since the numbers $n$ and $m$ of species are supposed constant, working with
the vector processes
$(P^{K,i}_{.},A^{K,j}_{.})_{i\in\left[\\!\left[1,n\right]\\!\right],j\in\left[\\!\left[1,m\right]\\!\right]}$
with values in $\mathbb{R}_{+}^{n+m}$ (and not the measure-valued processes)
is here sufficient. Notice that under the assumption (2.1), we can obtain by
computations similar to the Step 1 of the proof of Proposition 1.7 that:

$\sup_{K\in\mathbb{N}}\mathbb{E}\Big{[}\sup_{t\in[0,T]}\langle\mathbf{P}^{K,n,m}_{t},\mathds{1}\rangle^{3}+\langle\mathbf{A}^{K,n,m}_{t},\mathds{1}\rangle^{3}\Big{]}<+\infty.$
(2.4)

The global existence and uniqueness of solution to (2.3) can be proved by
classical results for ordinary differential equations: it follows from the
local boundedness and Lipschitz property of the functions on the right-hand
side. Then this proposition is a direct application of Theorem 2.1 p.456 in
the book by Ethier and Kurtz [22]. ∎

As a consequence:

######  Corollary 2.2.

Under the assumptions of Proposition 2.1, and for any $T>0$, the sequence of
measure-valued processes
$(\mathbf{P}^{K,n,m},\mathbf{A}^{K,n,m})_{K\in\mathbb{N}}$ converges a.s. and
uniformly in the Skorohod space
$\mathbb{D}\left([0,T],\mathcal{M}_{F}([0,1])^{2}\right)$ to the process
$({\mathbf{\widetilde{P}}}^{n},{\mathbf{\widetilde{A}}}^{m})$ defined by:

${\mathbf{\widetilde{P}}}^{n}_{t}(dx)=\frac{1}{n}\sum_{i=1}^{n}\widetilde{P}^{i}_{t}\delta_{x^{i}}(dx),\qquad{\mathbf{\widetilde{A}}}^{m}_{t}(dy)=\frac{1}{m}\sum_{j=1}^{m}\widetilde{A}^{j}_{t}\delta_{y^{j}}(dy).$

In this corollary, the space $\mathcal{M}_{F}([0,1])^{2}$ can be embedded with
the total variation topology which is stronger than the weak topology: as $n$,
$m$, the sequences $(x^{1},\dots x^{n})$ and $(y^{1},\dots y^{m})$ remain
unchanged, all the measures are absolutely continuous with respect to the same
counting measures $\sum_{i=1}^{n}\delta_{x^{i}}$ or
$\sum_{j=1}^{m}\delta_{y^{j}}$. Also, the uniform convergence in the Skorohod
space for any $T>0$ yields the convergence for the Skorohod topology on
$\mathbb{R}_{+}$ (see Theorem 16.2 in [11]).

For the remainder, let us denote by
$\Phi^{P,n,m}=(\Phi^{P,n,m}_{1},\cdots\Phi^{P,n,m}_{n})$ and
$\Phi^{A,n,m}=(\Phi^{A,n,m}_{1},\cdots\Phi^{A,n,m}_{m})$ the applications from
$\mathbb{R}_{+}^{n+m}$ into $\mathbb{R}^{n}$ and $\mathbb{R}^{m}$ respectively
such that (2.3) rewrites for all $i\in\\{1,\dots n\\}$ and $j\in\\{1,\dots
m\\}$:

$\displaystyle\frac{d\widetilde{P}^{i}_{t}}{dt}=\Phi^{P,n,m}_{i}(\widetilde{P}_{t},\widetilde{A}_{t})\quad\mbox{
and
}\quad\frac{d\widetilde{A}^{j}_{t}}{dt}=\Phi^{A,n,m}_{j}(\widetilde{P}_{t},\widetilde{A}_{t}).$
(2.5)

###  2.2 Central limit theorem

The individual based model (mathematically described by the SDE of Definition
1.6) is very commonly encountered in simulations in Biology for simulation,
while the ODEs (2.3) are often used in the modelling papers. The law of large
number (Proposition 2.1) shows that both point of view are related in the
limit $K\rightarrow+\infty$. In this section, we quantify the speed at which
the convergence holds by establishing a central limit theorem. This can be
useful for instance to compute confidence intervals or study the statistical
properties of parameter estimators.

Let us introduce the fluctuation processes:

$\eta^{K,P}_{t}=\sqrt{K}\left(\begin{array}[]{c}P^{K,1}_{t}/K-\widetilde{P}^{1}_{t}\\\
\vdots\\\ P^{K,n}_{t}/K-\widetilde{P}^{n}_{t}\end{array}\right),\qquad\mbox{
and
}\eta^{K,A}_{t}=\sqrt{K}\left(\begin{array}[]{c}A^{K,1}_{t}/K-\widetilde{A}^{1}_{t}\\\
\vdots\\\ A^{K,m}_{t}/K-\widetilde{A}^{m}_{t}\end{array}\right).$ (2.6)

######  Proposition 2.3.

Consider the process of Definition 1.6, assume that the functions $b^{A}$,
$b^{P}$, $d^{A}$ and $d^{P}$ are of class $\mathcal{C}^{1}$, and that the
initial conditions $(\eta^{K,P}_{0},\eta^{K,A}_{0})_{K\geq 0}$ converge in
distribution towards a deterministic vector
$(\widetilde{\eta}^{P}_{0},\widetilde{\eta}^{A}_{0})$ when
$K\rightarrow+\infty$. Under the same assumptions as Proposition 2.1, one has

$(\eta^{K,P}_{t},\eta^{K,A}_{t})_{t\geq
0}\underset{K\to\infty}{\Rightarrow}(\widetilde{\eta}^{P}_{t},\widetilde{\eta}^{A}_{t})_{t\geq
0},$

where the converge holds in law in the Skorohod space
$\mathbb{D}(\mathbb{R}_{+},\mathbb{R}^{n+m})$ and where the processes
$(\widetilde{\eta}^{P}_{.},\widetilde{\eta}^{A}_{.})=(\widetilde{\eta}^{P,i}_{.},\widetilde{\eta}^{A,j}_{.})_{i\in\\{1,\dots
n\\},j\in\\{1,\dots m\\}}$ are solutions of the following SDEs driven by $n+m$
independent standard Brownian motions $(W^{P,i},W^{A,j})_{1\leq i\leq n,1\leq
j\leq m}$ in $\mathbb{R}$:

for all $1\leq i\leq n$ and $1\leq j\leq m$,

$\displaystyle\widetilde{\eta}^{P,i}_{t}=$
$\displaystyle\widetilde{\eta}^{P,i}_{0}+\int_{0}^{t}\sqrt{\bigg{(}b^{P}\big{(}\widetilde{R}^{A,i}_{s}\big{)}+d^{P}\big{(}\widetilde{R}^{A,i}_{s}\big{)}+\frac{1}{n}\sum_{\ell=1}^{n}k(x^{i},x^{\ell})\widetilde{P}^{\ell}_{s}\bigg{)}\widetilde{P}^{i}_{s}}\
dW^{P,i}_{s}$ (2.7)
$\displaystyle+\int_{0}^{t}\left[\sum_{\ell=1}^{n}\frac{\partial\Phi^{P,n,m}_{i}}{\partial
p^{\ell}}(\widetilde{P}_{s},\widetilde{A}_{s})\widetilde{\eta}^{P,\ell}_{s}+\sum_{\ell=1}^{m}\frac{\partial\Phi^{P,n,m}_{i}}{\partial
a^{\ell}}(\widetilde{P}_{s},\widetilde{A}_{s})\widetilde{\eta}^{A,\ell}_{s}\right]ds,$
$\displaystyle\widetilde{\eta}^{A,j}_{t}=$
$\displaystyle\widetilde{\eta}^{A,j}_{0}+\int_{0}^{t}\sqrt{\bigg{(}b^{A}\big{(}\widetilde{R}^{P,j}_{s}\big{)}+d^{A}\big{(}\widetilde{R}^{P,j}_{s}\big{)}+\frac{1}{m}\sum_{\ell=1}^{m}h(y^{j},y^{\ell})\widetilde{A}^{\ell}_{s}\bigg{)}\widetilde{A}^{j}_{s}}\
dW^{A,j}_{s}$
$\displaystyle+\int_{0}^{t}\left[\sum_{\ell=1}^{n}\frac{\partial\Phi^{A,n,m}_{j}}{\partial
p^{\ell}}(\widetilde{P}_{s},\widetilde{A}_{s})\widetilde{\eta}^{P,\ell}_{s}+\sum_{\ell=1}^{m}\frac{\partial\Phi^{A,n,m}_{j}}{\partial
a^{\ell}}(\widetilde{P}_{s},\widetilde{A}_{s})\widetilde{\eta}^{A,\ell}_{s}\right]ds.$

Notice that the limiting fluctuation process is an Ornstein-Uhlenbeck process
that is centered. Systems of SDEs have already been introduced in the
literature to describe the evolution of communities, but they are to our
knowledge of a different nature. In (2.7) the noise relates to the fluctuation
of the stochastic individual-based model around its deterministic limit (2.5)
for $K\rightarrow+\infty$. In other works, such as in [14] for instance, the
white noise corresponds to a diffusive limit obtained when considering a
different longer time-scale, as in the Donsker theorem (see e.g. [16, Section
4.2]): the random noise comes from the rapid successions of birth and death
events in this accelerated time-scale. In recent works, following the steps of
May [40], [3, 26]introduce a system of equations, coupled via a smooth random
vector field, which describes the approximated dynamical system around an
equilibrium state: in this case, the noise models the complexity and
nonlinearity of interactions.

###### Proof of Proposition 2.3.

It is a direct application of Theorem 2.3 from Chapter 11 of [22]. The proof
can also be carried from the semi-martingale expressions of $P^{K,i}_{t}$ and
$A^{K,j}_{t}$ (1.13)-(1.15) and (2.5). To understand (2.7), consider for
example the $i$th component of $\eta^{K,P}_{t}$, $\eta^{K,P,i}_{t}$. The other
terms can be treated similarly. We sketch here this alternative proof:

$\displaystyle\eta^{K,P,i}_{t}=$
$\displaystyle\eta^{K,P,i}_{0}+\sqrt{K}\int_{0}^{t}\frac{1}{n}\sum_{\ell=1}^{n}k(x^{\ell},x^{i})\Big{(}{\color[rgb]{0,0,0}\widetilde{P}^{\ell}_{s}\widetilde{P}^{i}_{s}-\frac{P^{K,\ell}_{s}}{K}\frac{P^{K,i}_{s}}{K}}\Big{)}\
ds+\sqrt{K}\frac{M^{K,i}_{t}}{K}$ $\displaystyle+$
$\displaystyle\sqrt{K}\int_{0}^{t}\Big{[}g^{P}\Big{(}\sum_{j\sim
i}c_{ij}^{n,m}\frac{A^{K,j}_{s}}{K}\Big{)}{\color[rgb]{0,0,0}\frac{P^{K,i}_{s}}{K}}-g^{P}\Big{(}\sum_{j\sim
i}c_{ij}^{n,m}\widetilde{A}^{j}_{s}\Big{)}\widetilde{P}^{i}_{s}\Big{]}ds$
$\displaystyle=$
$\displaystyle\eta^{K,P,i}_{0}+\frac{1}{n}\sum_{\ell=1}^{n}k(x^{\ell},x^{i})\int_{0}^{t}{\color[rgb]{0,0,0}\Big{[}\eta^{K,P,\ell}_{s}\frac{P^{K,i}_{s}}{K}+\widetilde{P}^{\ell}_{s}\eta^{K,P,i}_{s}\Big{]}}ds+\sqrt{K}\frac{M^{K,i}_{t}}{K}$
(2.8) $\displaystyle+$
$\displaystyle\int_{0}^{t}\Big{[}g^{P}\Big{(}\sum_{j\sim
i}c_{ij}^{n,m}\frac{A^{K,j}_{s}}{K}\Big{)}{\color[rgb]{0,0,0}\eta^{K,P,i}_{s}}+\sqrt{K}\widetilde{P}^{i}_{s}\Big{(}g^{P}\Big{(}\sum_{j\sim
i}c_{ij}^{n,m}\frac{A^{K,j}_{s}}{K}\Big{)}-g^{P}\Big{(}\sum_{j\sim
i}c_{ij}^{n,m}\widetilde{A}^{j}_{s}\Big{)}\Big{)}\Big{]}ds,$

for the martingale $M^{K,i}$ appearing in (1.13) with quadratic variation
(1.14). Because the birth and death rates are assumed to be of class
$\mathcal{C}^{1}$, so is $g^{P}$ and a Taylor expansion can be used for the
last term:

$\displaystyle\sqrt{K}\Big{(}g^{P}\Big{(}\sum_{j\sim
i}c_{ij}^{n,m}\frac{A^{K,j}_{s}}{K}\Big{)}-g^{P}\Big{(}\sum_{j\sim
i}c_{ij}^{n,m}\widetilde{A}^{j}_{s}\Big{)}\Big{)}=$
$\displaystyle(g^{P})^{\prime}\Big{(}\sum_{j\sim
i}c_{ij}^{n,m}\widetilde{A}^{j}_{s}\Big{)}\times\sum_{j\sim
i}c_{ij}^{n,m}\sqrt{K}\Big{(}\frac{A^{K,j}_{s}}{K}-\widetilde{A}^{j}_{s}\Big{)}+\varepsilon_{K}$
$\displaystyle=$ $\displaystyle(g^{P})^{\prime}\Big{(}\sum_{j\sim
i}c_{ij}^{n,m}\widetilde{A}^{j}_{s}\Big{)}\times\sum_{j\sim
i}c_{ij}^{n,m}{\color[rgb]{0,0,0}\eta^{K,A,j}_{s}}+\varepsilon_{K},$ (2.9)

where $\varepsilon_{K}$ is a remainder term. From (2.8) and (2.9), we
recognize that:

$\eta^{K,P,i}_{t}=\eta^{K,P,i}_{0}+\sqrt{K}\frac{M^{K,i}_{t}}{K}+\int_{0}^{t}\Big{[}\sum_{\ell=1}^{n}\frac{\partial\Phi^{P,n,m}_{i}}{\partial
p^{\ell}}(\widetilde{P}_{s},\widetilde{A}_{s})\widetilde{\eta}^{P,\ell}_{s}+\sum_{\ell=1}^{m}\frac{\partial\Phi^{P,n,m}_{i}}{\partial
a^{\ell}}(\widetilde{P}_{s},\widetilde{A}_{s})\widetilde{\eta}^{A,\ell}_{s}\Big{]}ds+\varepsilon_{K}.$

Using the Aldous-Rebolledo criterion (e.g. [36]) it is possible to prove that
the distributions of the processes $(\eta^{K,P},\eta^{K,A})_{K}$ form a tight
family with a unique limiting value that solves the SDEs (2.7). ∎

##  3 Continuous limits when abundances and the number of species are large

We now consider that the numbers of plant and pollinator species in the
network tend to infinity. Taking the limit $n,m\rightarrow+\infty$, we obtain
equations describing the evolution of a population consisting in a continuum
of species.

Let the traits of plants and pollinators be chosen according to i.i.d. random
variables with cumulative distribution functions $F_{P}$ and $F_{A}$
respectively. To order species according to their respective trait, we proceed
as follows: let $(\tilde{u}_{i})_{i\geq 1}$ and $(\tilde{v}_{j})_{j\geq 1}$ be
two sequences of i.i.d. random variables with uniform distribution in $[0,1]$;
for any $n,m\in\mathbb{N}^{2}$, let $(u^{i,n})_{i=1..n}$ be the ordered
$n^{th}$ first values of $(\tilde{u}_{i})_{i\geq 1}$ and $(v^{j,m})_{j=1...m}$
be the ordered $m^{th}$ first values of $(\tilde{v}_{j})_{j\geq 1}$, then for
any $i\in\\{1,..,n\\}$ and $j\in\\{1,..,m\\}$,

$x^{i,n}=F_{P}^{-1}(u^{i,n})\qquad\mbox{ and }\qquad
y^{j,m}=F_{A}^{-1}(v^{j,m}).$

The indices $n$ and $m$ will be dropped when no confusion is possible.

Recall that the plant-pollinator network is defined by its adjacency matrix
$(G^{n,m}_{ij})_{(i,j)\in\left[\\!\left[1,n\right]\\!\right]\times\left[\\!\left[1,m\right]\\!\right]}$
and the weight of $c^{n,m}_{ij}$ of the interaction $i\sim j$ whenever
$G^{n,m}_{ij}=1$. The random variables $G^{n,m}_{ij}$ are supposed to satisfy
Assumption 1.1 for all $n$ and $m\in\mathbb{N}$: the entries of the adjacency
matrix are here independent Bernoulli with parameters $\phi(x^{i},y^{j})$ for
the term
$(i,j)\in\left[\\!\left[1,n\right]\\!\right]\times\left[\\!\left[1,m\right]\\!\right]$.
For the harvesting coefficients $c^{n,m}_{ij}$, the following assumption is
made:

######  Assumption 3.1.

We assume that there exists for any $n,m\in\mathbb{N}^{2}$ a continuous
function $c^{n,m}$ satisfying the condition (1.2) of Assumption 1.3.
Additionally, we assume that there exists a function
$c\,:\,[0,1]^{2}\mapsto\mathbb{R}$ to which the sequences of functions
$nc^{n,m}(.,.)$ and $mc^{n,m}(.,.)$ converge uniformly, in
$L^{\infty}([0,1]^{2})$, when $n$ and $m$ tend to $+\infty$.

The idea in this assumption is that the function $c(.,.)$ is a “harvesting
function” underlying the matrix
$(c_{ij}^{n,m})_{i\in\left[\\!\left[1,n\right]\\!\right],j\in\left[\\!\left[1,m\right]\\!\right]}$
of harvesting coefficients. For instance, Examples 1.4 satisfy the previous
assumptions.

Notice also that Assumption 3.1 implies that $n$ and $m$ grow to infinity with
a similar speed. To state this idea in specific terms, we assume that there
exists a sequence $(\alpha_{n})_{n\geq 1}$ such that

$m=\alpha_{n}n\quad\text{and}\quad\lim_{n\to\infty}\alpha_{n}=1.$ (3.1)

In what follows, we write $n\to\infty$, $m\to\infty$ or $n,m\to\infty$
indistinctly.

######  Proposition 3.2.

Let us assume Assumptions 1.1, 1.5 and 3.1 hold and that there exist
deterministic continuous bounded densities $\bar{p}_{0}$ and $\bar{a}_{0}$
such that the following weak convergences hold:

$\lim_{n\rightarrow+\infty}\frac{1}{n}\sum_{i=1}^{n}\widetilde{P}^{i}_{0}\delta_{x^{i}}\stackrel{{\scriptstyle
w}}{{=}}\bar{p}_{0}(x)dx,\quad\lim_{m\rightarrow+\infty}\frac{1}{m}\sum_{j=1}^{m}\widetilde{A}^{j}_{0}\delta_{y^{j}}\stackrel{{\scriptstyle
w}}{{=}}\bar{a}_{0}(y)dy\quad\text{a.s.}$ (3.2)

For any $T\geq 0$, and for $n,m\rightarrow+\infty$, the sequence of measure-
valued processes

$\left(\widetilde{\mathbf{P}}^{n,m}_{t}(dx),\widetilde{\mathbf{A}}_{t}^{n,m}(dy)\right)_{t\geq
0}:=\Big{(}\frac{1}{n}\sum_{i=1}^{n}\widetilde{P}^{i}_{t}\delta_{x^{i}},\
\frac{1}{m}\sum_{j=1}^{m}\widetilde{A}^{j}_{t}\delta_{y^{j}}\Big{)}_{t\geq
0},\quad n,m\geq 1$

converge in law in $\mathcal{C}([0,T],\mathcal{M}_{F}([0,1])^{2})$ to a
deterministic process $(\bar{\mathbf{P}},\bar{\mathbf{A}})$ in
$\mathcal{C}([0,T],\mathcal{M}_{F}([0,1])^{2})$, such that:
(i) for all $t\geq 0$, $\bar{\mathbf{P}}_{t}$ and $\bar{\mathbf{A}}_{t}$ admit
densities $\bar{p}_{t}$ and $\bar{a}_{t}$ with respect to the Lebesgue measure
on $[0,1]$,
(ii) $(\bar{\mathbf{P}},\bar{\mathbf{A}})$ is the unique solution, for
$f\in\mathcal{C}([0,1],\mathbb{R})$, of

$\displaystyle\int_{0}^{1}f(x)d\bar{\mathbf{P}}_{t}(x)=$
$\displaystyle\int_{0}^{1}f(x)\bar{p}_{0}(x)dx$ (3.3)
$\displaystyle+\int_{0}^{t}\int_{0}^{1}f(x)\Big{[}g^{P}\Big{(}\int_{0}^{1}c(x,y)\phi(x,y)\bar{a}_{s}(y)dy\Big{)}-k\star\bar{p}_{s}(x)\Big{]}\bar{p}_{s}(x)dx\
ds,$ $\displaystyle\int_{0}^{1}f(y)d\bar{\mathbf{A}}_{t}(y)=$
$\displaystyle\int_{0}^{1}f(y)\bar{a}_{0}(y)dy$
$\displaystyle+\int_{0}^{t}\int_{0}^{1}f(y)\Big{[}g^{A}\Big{(}\int_{0}^{1}c(x,y)\phi(x,y)\bar{p}_{s}(x)dx\Big{)}-h\star\bar{a}_{s}(y)\Big{]}\bar{a}_{s}(y)dy\
ds,$

where $k\star\nu$ denotes in this case the convolution on two functions and
$\phi$ was introduced in Assumption 1.1. In all this statement, the space
$\mathcal{M}_{F}^{2}([0,1])$ is endowed with its weak topology.

The proof of this Proposition is given in Appendix A.2

Equation (3.3) is analogous to the ODE system given in Equation (2.3).
However, here species are not considered as discrete but continuously
distributed along a continuous trait. Modelling the dynamics of a plant-
pollinator community with the PDEs given in Equation (3.3) is thus a
functional rather than a species representation of the same system. The
connexions are modelled in (3.3) by $\phi$. This function
$\phi\,:\,[0,1]^{2}\mapsto[0,1]$ is a graphon (see [39]): it can be understood
(in this case) as a graph on node sets $[0,1]$ for the plants and $[0,1]$ for
the pollinators, where $\phi(x,y)$ describes the density of connexions between
plants $x$ and pollinator $y$. The term $c\phi$ reflects both the topology
($\phi$) and the intensity ($c$) of plant-pollinator interactions throughout
the community depending on the traits values $x$ and $y$ involved.

######  Remark 3.3.

Note that one can formally write a strong form of (3.3). This leads to the
following system of integral equations on the densities $\bar{p}_{t}$ and
$\bar{a}_{t}$:

$\displaystyle\partial_{t}\bar{p}_{t}(x)=$
$\displaystyle\Big{[}g^{P}\Big{(}\int_{0}^{1}{c(x,y)}\phi(x,y)\bar{a}_{t}(y)dy\Big{)}-k\star\bar{p}_{t}(x)\Big{]}\bar{p}_{t}(x),$
(3.4) $\displaystyle\partial_{t}\bar{a}_{t}(y)=$
$\displaystyle\Big{[}g^{A}\Big{(}\int_{0}^{1}{c(x,y)}\phi(x,y)\bar{p}_{t}(x)dx\Big{)}-h\star\bar{a}_{t}(y)\Big{]}\bar{a}_{t}(y),$

with initial conditions $\bar{p}_{0}$ and $\bar{a}_{0}$.

##  4 Study of the limiting dynamical systems

Let us now study the behaviour of the limiting dynamical systems that have
been obtained: the ODE system (2.3) when $K\rightarrow+\infty$, and the
kinetic PDEs (3.3) when additionally $n,m\rightarrow+\infty$.

###  4.1 Stationary states of the ODE system (2.3)

In this section, we give some results about the dynamics of solutions to
Systems (2.3), for general forms of $g^{P}$ and $g^{A}$. Notice that any
stationary points
$(\widetilde{P}^{1}_{\infty},\dots\widetilde{P}^{n}_{\infty},\widetilde{A}^{1}_{\infty},\dots\widetilde{A}^{m}_{\infty})$
of (2.3) are solutions to the following system:

$\displaystyle\forall 1\leq i\leq n,$
$\displaystyle\widetilde{P}^{i}_{\infty}=0\mbox{ or }g^{P}\big{(}\sum_{j\sim
i}c_{ij}^{n,m}\widetilde{A}^{j}_{\infty}\big{)}=\frac{1}{n}\sum_{\ell=1}^{n}k(x^{i},x^{\ell})\widetilde{P}^{\ell}_{\infty},$
(4.1) $\displaystyle\forall 1\leq j\leq m,$
$\displaystyle\widetilde{A}^{j}_{\infty}=0\mbox{ or }g^{A}\big{(}\sum_{i\sim
j}c_{ij}^{n,m}\widetilde{P}^{i}_{\infty}\big{)}=\frac{1}{m}\sum_{\ell=1}^{m}h(y^{j},y^{\ell})\widetilde{A}^{\ell}_{\infty}.$

The community with no individuals is an obvious stationary point called the
null equilibrium in the rest of the paper:

$\widetilde{P}^{i}_{\infty}=0\text{ for all
}i\in\\{1,..n\\}\quad\text{and}\quad\widetilde{A}^{j}_{\infty}=0\text{ for all
}j\in\\{1,..m\\},$

is a stationary state of (2.3). The local stability of this stationary state
depends only on the sign of $g^{P}(0)$ and $g^{A}(0)$.

######  Remark 4.1.

According to the expression of the non trivial equilibria, a plant-pollinator
community can exist, i.e. $\widetilde{P}^{i}_{\infty}>0$ and
$\widetilde{A}^{j}_{\infty}>0$, as soon as the interactions between plants and
pollinators translate into positive growth rate for both plants and
pollinators. Despite this condition is little restrictive, a non-trivial
equilibrium is possible only if trajectories have to converge to it, (i.e. a
non-trivial equilibrium has to be stable). The stability of an equilibrium can
be more tricky to study and obtained in general.

######  Lemma 4.2.

If $g^{P}(0)\vee g^{A}(0)<0$, the null equilibrium is locally stable.
If $g^{P}(0)\vee g^{A}(0)>0$, the null equilibrium is locally unstable.

###### Proof.

From the Hartman–Grobman theorem, we know that local stability is given by the
signs of the Jacobian matrix eigenvalues. The Jacobian matrix around the null
equilibrium can be directly computed:

$\begin{pmatrix}g^{P}(0)&0&0&...&0&0\\\ 0&g^{P}(0)&0&...&0&0\\\
...&...&...&...&...&...\\\ 0&0&...&0&g^{A}(0)&0\\\ 0&0&...&0&0&g^{A}(0)\\\
\end{pmatrix},$

whose eigenvalues are $g^{P}(0)$ with multiplicity $n$ and $g^{A}(0)$ with
multiplicity $m$. This ends the proof. ∎

###  4.2 Study of a plant-pollinator interaction with a trade-off for plants

As explained in Remark 4.1, more precise functional forms for the birth and
death rates are needed to carry further computation.

In this Section, we give results for the following particular forms for the
individual growth rate functions:

$\displaystyle g^{P}(R^{A})=$ $\displaystyle
b^{P}(R^{A})-d^{P}(R^{A})=\frac{\alpha_{P}R^{A}}{\beta_{P}+\gamma_{P}R^{A}}-(d_{P}+\delta_{P}R^{A}),$
(4.2) $\displaystyle g^{A}(R^{P})=$ $\displaystyle
b^{A}(R^{P})-d^{A}(R^{P})=\frac{\alpha_{A}R^{P}}{\beta_{A}+\gamma_{A}R^{P}}-d_{A},$

where $R^{A}$ and $R^{P}$ are the total resources collected by respectively
plants and pollinators (Eqs. (1.5) and (1.6)). Parameters $\alpha_{P}$,
$\alpha_{A}$, $\beta_{P}$, $\beta_{A}$, $\gamma_{P}$, $\gamma_{A}$, $d_{A}$,
$d_{P}$ and $\delta_{P}$ are assumed positive. Note that according to Lemma
4.2, the null equilibrium is stable in this particular case.

On the right-hand side of Equation (4.2), the death rate is an increasing
function of the resources exchanged $R$ for plants. This reflects an
interaction trade-off for the plant, i.e. it is supposed that there is a cost
for interacting with pollinators due to nectar production, leaves consumption,
etc. Graphical representations of $g^{P}$ and $g^{A}$ are given in Figure 2.

Figure 2:  Graphical representation of $g_{A}$ and $g_{P}$

The form given here to the growth rates (or numerical responses) $g^{P}$ and
$g^{A}$ can be found in [30]. It assumes that there is no competition between
plants for collecting resources from pollinators, and that there is no
competition between pollinators for collecting resources from plants.

We state an additional assumption on the parameters of both functions to avoid
the case where one growth rate is never positive.

######  Assumption 4.3.

Assume that there exist $r^{P}>0$ and $r^{A}>0$ such that

$g^{P}(r^{P})>0\quad\text{and}\quad g^{A}(r^{A})>0.$

####  4.2.1 Case $n=m=1$

Let us first consider the particular case of a single plant species
interacting with a single pollinator species. In other words, we consider the
system (2.3) with $n=m=1$ and omit the indices and exponents $i=1$ and $j=1$
for the sake of simplicity (so that $\widetilde{P}^{1}_{t}$ becomes ${P}_{t}$
and $c_{11}^{1,1}$ becomes $c$ for instance). We also assume that the two
species interact, i.e. $G_{11}=1$. Otherwise, computations are trivial: the
community goes to extinction as only death events occur. We have then:

$\displaystyle\frac{d{P}_{t}}{dt}$
$\displaystyle=\left(g^{P}\big{(}c{A}_{t}\big{)}-k{P}_{t}\right){P}_{t}$ (4.3)
$\displaystyle\frac{d{A}_{t}}{dt}$
$\displaystyle=\left(g^{A}\big{(}c{P}_{t}\big{)}-h{A}_{t}\right){A}_{t}$

(see [32] for a similar ODE model in the ecology literature).

As discussed in Lemma 4.2, the null equilibrium is a locally stable
equilibrium of the system (4.3). Let us search for positive equilibrium. To
help the study, we may draw some phase plan of the system, for the choice of
parameters $\alpha_{A}=25$, $\alpha_{P}=9$,
$\beta_{A}=\beta_{P}=\gamma_{A}=\gamma_{P}=1$ and $\delta_{P}=3$, and some
various $d_{A}$ and $d_{P}$. (see Figure 3).

Figure 3: Phase plan and nullclines of the system of ODE (4.3): nullclines for
the pollinators dynamics in blue; nullclines for the plants dynamics in cyan.
Right: $d_{A}=2$ and $d_{P}=1$, the system has $3$ stationary states: the null
equilibrium, $1$ stable positive equilibrium and $1$ unstable positive
equilibrium. Left: $d_{A}=3$ and $d_{P}=1.2$, the unique equilibrium of the
system is the null equilibrium.

Precisely, we prove the following lemma that gives the number of equilibria
and their stability.

######  Lemma 4.4.

In addition to the null-equilibrium, System (4.3) has $0$, $1$ or $2$
equilibria, depending on the values of the parameters.
Precisely, denote by $C_{0}^{-}$ and $C_{0}^{+}$ the two zeros of $g^{P}$.
System (4.5) has $2$ positive equilibrium if and only if

$\max_{x\in(C^{-}_{0}/c,C^{+}_{0}/c)}f(x)>0.$

In this case, the equilibria are $({P}^{-}_{\infty},{A}^{-}_{\infty})$ and
$({P}^{+}_{\infty},{A}^{+}_{\infty})$ with ${A}^{-}_{\infty}<{A}^{+}_{\infty}$
and are the two positive solutions to

$\left\\{\begin{aligned} &{P}_{\infty}=\frac{1}{k}g^{P}(c{A}_{\infty})\\\
&{A}_{\infty}=\frac{1}{h}g^{A}\left(\frac{c}{k}g^{P}(c{A}_{\infty})\right),\text{
and }c{A}_{\infty}\in(C_{0}^{-},C_{0}^{+}).\end{aligned}\right.$ (4.4)

$({P}^{-}_{\infty},{A}^{-}_{\infty})$ is an unstable positive equilibrium and
$({P}^{+}_{\infty},{A}^{+}_{\infty})$ is a stable positive equilibrium.

To conclude on the dynamics of the trajectories of System (4.3), we performed
simulations. In summary, it appeared that, depending on the parameters values,
either all trajectories are attracted by the null equilibrium (Fig 3(left)),
or trajectories converge to $0$ or some positive equilibrium, depending on the
initial conditions (Fig 3(right)). We never observed cycles. Moreover, the
competitive terms ensure that trajectories remain bounded.

###### Proof of Lemma 4.4.

Any stationary point $({P}_{\infty},{A}_{\infty})$ is solution to

$\displaystyle g^{P}\big{(}c{A}_{\infty}\big{)}=k{P}_{\infty}$
$\displaystyle\quad\text{ or }\quad{P}_{\infty}=0$ (4.5) $\displaystyle
g^{A}\big{(}c{P}_{\infty}\big{)}=h{A}_{\infty}$ $\displaystyle\quad\text{ or
}\quad{A}_{\infty}=0.$

Recall that either both or no coordinates are null. In view of (4.5), the
existence of a positive equilibrium requires at least that
$\sup_{x\in\mathbb{R}^{+}}g^{P}(x)>0$ and
$\sup_{y\in\mathbb{R}^{+}}g^{A}(y)>0$. This is true under Assumption 4.3.
Recall that $C_{0}^{-}$ and $C_{0}^{+}$ are the two zeros of $g^{P}$. From
(4.5), we deduce that, if there exists a positive stationary state, it
satisfies System (4.4).

It remains to find the number of zeros of
$f(x)=\frac{1}{h}g^{A}\left(\frac{c}{k}g^{P}(cx)\right)-x$ on
$\left(\frac{C^{-}_{0}}{c},\frac{C^{+}_{0}}{c}\right)$. A rapid
differentiation gives

$f^{{}^{\prime\prime}}(x)=\frac{c^{3}}{hk}(g^{P})^{{}^{\prime\prime}}(cx)(g^{A})^{\prime}\left(\frac{c}{k}g^{P}(cx)\right)+\frac{c^{3}}{(hk)^{2}}(g^{P})^{{}^{\prime}}(cx)^{2}(g^{A})^{{}^{\prime\prime}}\left(\frac{c}{k}g^{P}(cx)\right)<0,$

since $g^{A}$ is increasing and $g^{P}$ and $g^{A}$ are concave functions.
Thus, $f^{\prime}$ is decreasing and

$f^{{}^{\prime}}(x)=\frac{c^{2}}{hk}(g^{P})^{{}^{\prime}}(cx)(g^{A})^{\prime}\left(\frac{c}{k}g^{P}(cx)\right)-1,$

converges to $+\infty$ (resp. $-\infty$) when $x$ goes to $C^{-}_{0}/c$ (resp.
$C^{+}_{0}/c$). Finally, ${f}$ is increasing then decreasing and admits a
unique maximum. Since ${f}$ converges to $-\infty$ when $x$ goes to
$C^{-}_{0}/c$ and $C^{+}_{0}/c$, it admits

  * •

$2$ zeros if $\max_{x\in(C^{-}_{0}/c,C^{+}_{0}/c)}f(x)>0$, and System (4.5)
has $2$ positive equilibrium,

  * •

$1$ zero if $\max_{x\in(C^{-}_{0}/c,C^{+}_{0}/c)}f(x)=0$, and System (4.5) has
$1$ positive equilibrium,

  * •

no zero otherwise, and System (4.5) has no positive equilibrium.

Finally, in the case of $2$ positive equilibrium, which will be denoted by
$({P}^{-}_{\infty},{A}^{-}_{\infty})$ and
$({P}^{+}_{\infty},{A}^{+}_{\infty})$ with
${A}^{-}_{\infty}<{A}^{+}_{\infty}$, the Jacobian matrix can be computed in
order to study their stability (see the Hartman–Grobman theorem):

$J^{\pm}=\begin{pmatrix}-k{P}^{\pm}_{\infty}&c{P}^{\pm}_{\infty}(g^{P})^{\prime}(c{A}^{\pm}_{\infty})\\\
c{A}^{\pm}_{\infty}(g^{A})^{\prime}(c{P}^{\pm}_{\infty})&-h{A}^{\pm}_{\infty}\end{pmatrix}.$

The trace is negative, thus the stability of the stationary states depend on
the sign of the determinant.

$\det(J^{\pm})={P}^{\pm}_{\infty}{A}^{\pm}_{\infty}\Big{(}hk-c^{2}(g^{P})^{\prime}(c{A}^{\pm}_{\infty})(g^{A})^{\prime}(c{P}^{\pm}_{\infty})\Big{)}=-{P}^{\pm}_{\infty}{A}^{\pm}_{\infty}hkf^{\prime}({A}^{\pm}_{\infty}),$

and its sign depends only on the sign of $f^{\prime}({A}^{\pm}_{\infty})$.
According to the previous study on function $f$, we deduce that $\det(J^{+})$
is negative, and $\det(J^{-})$ is positive. In other words, System (4.5)
admits $1$ stable positive equilibrium and $1$ unstable positive equilibrium.
In the case when System (4.5) admits a unique positive equilibrium, the same
study implies that this equilibrium is a non-hyperbolic equilibrium. ∎

####  4.2.2 Behaviour of the kinetic equations (3.3)

Recall that we are still working with the growth rates defined in (4.2).

The stationary solutions of (3.3) are couples of measures
$\bar{P}^{\infty}(dx)$ and $\bar{A}^{\infty}(dy)$ in $\mathcal{M}_{F}([0,1])$
such that, for all positive, bounded and continuous function $f$ on $[0,1]$:

$\displaystyle 0=$
$\displaystyle\int_{0}^{1}f(x)\Big{[}g^{P}\Big{(}\int_{0}^{1}\psi(x,y)\bar{A}^{\infty}(dy)\Big{)}-k\star\bar{P}^{\infty}(x)\Big{]}\bar{P}^{\infty}(dx),$
(4.6) $\displaystyle 0=$
$\displaystyle\int_{0}^{1}f(y)\Big{[}g^{A}\Big{(}\int_{0}^{1}\psi(x,y)\bar{P}^{\infty}(dx)\Big{)}-h\star\bar{A}^{\infty}(y)\Big{]}\bar{A}^{\infty}(dy)$
(4.7)

where

$\psi(x,y):=c(x,y)\phi(x,y).$ (4.8)

The null measures constitute a trivial solution to (4.6)-(4.7), which is
stable in our particular case (4.2). Let us discuss non-zero solutions.

######  Proposition 4.5.

Assume that

  * •

the competitive kernels $k$ and $h$ are constant functions;

  * •

for all $x_{0},y_{0}\in[0,1]$, $y\mapsto\psi(x_{0},y)$ and
$x\mapsto\psi(x,y_{0})$ are increasing and continuous functions;

then, System (3.3) does not admit non-null stationary state with densities
w.r.t Lebesgue measure.
Moreover, any non-null stationary state in $L^{1}([0,1]^{2})$ is a couple of
measures $(\bar{P}^{\infty},\bar{A}^{\infty})$ such that

$\displaystyle\exists\bar{a}_{0},$
$\displaystyle\bar{p}_{1}\in\mathbb{R}^{*}_{+},\bar{p}_{2}\in\mathbb{R}_{+},\bar{x}_{1},\bar{x}_{2},\bar{y}_{0}\in[0,1],\
\left\\{\begin{aligned}
&\bar{P}^{\infty}=\bar{p}_{1}\delta_{\bar{x}_{1}}+\bar{p}_{2}\delta_{\bar{x}_{2}}\\\
&\bar{A}^{\infty}=\bar{a}_{0}\delta_{\bar{y}_{0}}\end{aligned}\right.$ (4.9)
$\displaystyle\text{ with }\left\\{\begin{aligned}
&g^{P}\Big{(}\bar{a}_{0}\psi(\bar{x}_{1},\bar{y}_{0})\Big{)}=g^{P}\Big{(}\bar{a}_{0}\psi(\bar{x}_{2},\bar{y}_{0})\Big{)}=k(\bar{p}_{1}+\bar{p}_{2})\\\
&g^{A}\Big{(}\bar{p}_{1}\psi(\bar{x}_{1},\bar{y}_{0})+\bar{p}_{2}\psi(\bar{x}_{2},\bar{y}_{0})\Big{)}=h\bar{a}_{0}.\end{aligned}\right.$

All these stationary states are unstable, except the state

$\left\\{\begin{aligned}
&\bar{P}^{\infty}=\frac{\max_{\mathbb{R}^{+}}g^{P}}{k}\ \delta_{x_{0}}\\\
&\bar{A}^{\infty}=\frac{\arg\max_{\mathbb{R}^{+}}g^{P}}{\psi(x_{0},1)}\
\delta_{1}\end{aligned}\right.\\\ $ (4.10)

if $x_{0}$, solution to

$g^{A}\left(\frac{\max_{\mathbb{R}^{+}}g^{P}}{k}\psi(x_{0},1)\right)\psi(x_{0},1)=h\cdot\arg\max_{\mathbb{R}^{+}}g^{P},$
(4.11)

exists and is unique.
Finally, assuming that, for all initial conditions with positive densities
w.r.t Lebesgue measure, the quantities
$\int_{0}^{1}\psi(x,.)\bar{p}_{t}(x)dx$, $\int_{0}^{1}\bar{p}_{t}(x)dx$,
$\int_{0}^{1}\psi(.,y)\bar{a}_{t}(y)dy$ and $\int_{0}^{1}\bar{a}_{t}(y)dy$
converge when $t$ grow to infinity, then the trajectory converges to
equilibrium (4.10).

This proposition ensures the fact that the only possible stable equilibrium is
composed of only one plant species and one pollinator species. Simulations
below will illustrate this Proposition.

######  Remark 4.6.

The proof of Proposition 4.5 is not restricted to the specified forms of
$g^{P}$ and $g^{A}$, (4.2). The proposition is still true as soon as shapes
(successions of increases and decreases) of $g^{P}$ and $g^{A}$ are the same
as the ones of specified functions (4.2).
Otherwise, in any case, System (3.3) has no stationary state with densities,
all stationary states will be composed of Dirac measures. Moreover, the
maximal number of Dirac measures found in such a stationary state corresponds
to the maximal number of points that can be found in an inverse image of a
positive real for functions $g^{P}$ and $g^{A}$ respectively. The stability of
these stationary states can then be deduced using the same kind of arguments
as these of the following proof.
From a biological point of view, it means that that the number of pollinator
species and plant species is reduced to $1$ when considering all type of
growths given in Figure 2 of [31].
This result shows that when the plant-pollinator network is nested and the
competition among plants and among pollinators is constant, then the plant-
pollinator community collapses to a single plant-pollinator species pair. At
the stationary state, the value of the trait for the pollinators is 1, where
the probability of interactions with plant is maximized. The trait for plants
is <1 ((4.11)), because of the trade-offs for plants due to the cost of the
interaction. This is in line with the numerical analysis of a system of ODEs
by [38]. Our results is a formal demonstration of this necessary collapse. It
is also more general since we show that it does not depend on the specific
form of $g^{P}$ and $g^{A}$. This raises the question whether it is possible,
in this type of ecological model, to maintain a stable coexistence of many
plant and pollinator species in a single community by modifying the structure
of the interaction graph (i.e. with other assumptions for the function
$\psi$), or the structure of the competition graph (i.e. with non-constant
competition functions $k$ and $h$).

###### Proof.

Any stationary state in $L^{1}([0,1]^{2})$ satisfies that
$(\bar{P}^{\infty},\bar{A}^{\infty})$

$\left\\{\begin{aligned}
g^{P}\Big{(}\int_{0}^{1}\psi(x,y)\bar{A}^{\infty}(dy)\Big{)}=\int_{0}^{1}k\bar{P}^{\infty}(dx^{\prime}),\text{
for all }{\color[rgb]{0,0,0}x\in}\text{supp}\bar{P}^{\infty},\\\
g^{A}\Big{(}\int_{0}^{1}\psi(x,y)\bar{P}^{\infty}(dx)\Big{)}=\int_{0}^{1}h\bar{A}^{\infty}(dy^{\prime}),\text{
for all
}{\color[rgb]{0,0,0}y\in}\text{supp}\bar{A}^{\infty}\end{aligned}\right.$
(4.12)

Since $(\bar{P}^{\infty},\bar{A}^{\infty})\in L^{1}([0,1]^{2})$,
$\int_{0}^{1}k\bar{P}^{\infty}(dx^{\prime})$ and
$\int_{0}^{1}h\bar{A}^{\infty}(dy^{\prime})$ are finite. Then, according to
the continuity and the variations of the functions ($\psi$ is increasing w.r.t
each variable, $g^{A}$ is increasing and $g^{P}$ is increasing then
decreasing), we deduce (4.9).
Let us denote by $\mathcal{S}$ the set of stationary states given by (4.9)
excluding equilibrium (4.10) Let us prove that all states in $\mathcal{S}$ are
unstable. To this aim, we develop an argument by contradiction, assuming that
it is not true and that there exists at least a stationary state
$(\bar{P}^{\infty},\bar{A}^{\infty})\in\mathcal{S}$. We deal with the case
where $\bar{x}_{1}<x_{0}<\bar{x}_{2}$, $\bar{y}_{0}$ being any value in
$[0,1]$. Other cases ((2) $\bar{x}_{2}<x_{0}<\bar{x}_{1}$,
$\bar{y}_{0}\in[0,1]$ or (3) $\bar{x}_{1}=\bar{x}_{2}=x_{0}$, $\bar{y}_{0}<1$)
can be treated similarly.
Let consider an initial state that is close (in Wasserstein distance) to
$(\bar{P}^{\infty},\bar{A}^{\infty})$ and that has positive densities w.r.t
Lebesgue measures. Thus, $p_{t}$ exists and is positive for all $t\geq 0$ and
we can study $\log(\bar{p}_{t}(x))$ for all $t\in\mathbb{R}_{+}$ and all
$x\in[0,1]^{2}$:

$\frac{d}{dt}\log({\color[rgb]{0,0,0}\bar{p}}_{t}(x))=g^{P}\Big{(}\int_{0}^{1}\psi(x,y){\color[rgb]{0,0,0}\bar{a}}_{t}(y)dy\Big{)}-\int_{0}^{1}k,hel{\bar{p}}_{0}\,t(x^{\prime})dx^{\prime}.$

Under our assumption, $\mathcal{S}$ admits at least a stable state and since
the initial state is close to this set, we deduce that
$({\color[rgb]{0,0,0}\bar{p}}_{t},{\color[rgb]{0,0,0}\bar{a}}_{t})$ will
converge to some $(\hat{P}^{\infty},\hat{A}^{\infty})\in\mathcal{S}$, which is
close to $(\bar{P}^{\infty},\bar{A}^{\infty})$ (if not it), sufficiently close
to have a similar form. Then for all $x\in[0,1]$,

$\displaystyle\frac{d}{dt}\log({\color[rgb]{0,0,0}\bar{p}}_{t}(x))\underset{t\to\infty}{\longrightarrow}$
$\displaystyle
g^{P}\Big{(}\int_{0}^{1}\psi(x,y)\hat{a}^{\infty}(y)dy\Big{)}-\int_{0}^{1}k\hat{p}^{\infty}(x^{\prime})dx^{\prime}$
$\displaystyle=g^{P}\Big{(}\int_{0}^{1}\psi(x,y)\hat{a}^{\infty}(y)dy\Big{)}-g^{P}\Big{(}\int_{0}^{1}\psi(\hat{x}_{1},y)\hat{a}^{\infty}(y)dy\Big{)}.$

The latter quantity is positive as soon as $x\in(\hat{x}_{1},x_{0})$, which
contradicts the fact that ${\color[rgb]{0,0,0}\bar{p}}_{t}(x)$ converges to
$0$ for all $x\not\in\\{\hat{x}_{1},\hat{x}_{2}\\}$.
The last point of Proposition (4.5) can be obtained using the same argument,
once the convergences of the four quantities detailed in the statement are
assumed. ∎

Let us illustrate numerically Proposition 4.5. Let $(p_{i}(t))_{i}$, resp.
$(a_{j}(t))_{j}$ be an approximation of $(\bar{p}_{t}(x_{i}))_{i}$, resp.
$(\bar{a}(y_{j}))_{j}$, for $x_{i}=i/N$, $y_{j}=j/N$, $0\leq i,j\leq N$. Then
the continuous model (3.4) can be approximated using the rectangular rule by
the following system of coupled ODEs:

$\displaystyle\frac{dp_{i}}{dt}=\left[g^{P}\left(\frac{1}{N}\sum_{j=0}^{N}c_{ij}a_{j}\right)-\frac{1}{N}\sum_{j=0}^{N}k_{i,j}p_{j}\right]p_{i},$
(4.13)
$\displaystyle\frac{da_{j}}{dt}=\left[g^{A}\left(\frac{1}{N}\sum_{i=0}^{N}c_{ij}p_{i}\right)-\frac{1}{N}\sum_{i=0}^{N}h_{i,j}a_{i}\right]a_{j},$

where the initial data has been defined as $p_{i}(0)=\bar{p}_{0}(x_{i})$,
$p_{i}(0)=\bar{p}_{0}(x_{i})$, $a_{j}(0)=\bar{a}_{0}(x_{j})$. We recognize the
ODE system (2.3) obtained when the numbers of plant and pollinator species are
finite with $n=m=N$.

We shall consider the case of growth functions (4.2) with $\alpha_{A}=3$,
$\alpha_{P}=25$, $\beta_{A}=\beta_{P}=\gamma_{P}=1$, $\gamma_{A}=0.3$,
$d_{P}=1$, $d_{A}=3$ and $\delta_{P}=3$, with $N=500$ species of plants and
pollinators. In order to fit the hypotheses of Proposition 4.5 and insure the
convergence towards an explicit equilibrium of the dynamics, the interaction
matrix $(c_{ij})_{ij}$ is nondecreasing in each index $i$ and $j$, and the
competition kernels will be chosen constant. We consider the simple case

$c_{ij}=\frac{(i+1)(j+1)}{2N^{2}}.$

We represent in Figure 4 (in logarithmic scale) the convergence towards the
equilibrium state reached when no evolution occurs anymore on the discrete
dynamics, for a random initial datum supported in $[0,1]$. We observe that the
conclusion of proposition 4.5 are in agreement, namely the system converges
towards an equilibrium which consists of a Dirac delta centered inside the
domain for the plants, and on the right border for the pollinators: the last
remaining plant specie is moderately specialist ($x\simeq 0.6$) whereas the
remaining pollinator is genuinely generalist ($y=1$).

Figure 4: Solutions to (4.13) at time $100$ (blue), $500$ (orange) and $1500$
(green). The left panel represents the distribution of plant species and the
right panel the distribution of pollinator species. Parameters are
$\alpha_{A}=3$, $\alpha_{P}=25$, $\beta_{A}=\beta_{P}=\gamma_{P}=1$,
$\gamma_{A}=0.3$, $d_{P}=1$, $d_{A}=3$ and $\delta_{P}=3$, with $N=500$

##  5 Conclusion

Two of the main goals of ecology are i) to explain how emerged the size,
composition and structure of communities, for example plant-pollinator
communities, and ii) predict the dynamics and stability of a given community.
Such questions involve different hierarchical scales, from the individuals
which effectively interact (the microscopic scale), to the species and the
whole community (the macroscopic scale).

Theoretical ecology mostly address these questions by studying a system of
ODEs where one equation refers to a given species. This approach has several
methodological and conceptual drawbacks. First, the system is discretized and
structured by species, which precludes the possibility of within-species
variability regarding the rate and intensity of interactions, in particular
between-species overlaps. Second, the interaction graphs underlying the system
of ODEs are most often arbitrarily given without specified mechanisms.
Finally, due to the high-dimensionality of the ODEs system, it is difficult to
obtain general properties of communities.

Here, we modelled an individual-based plant-pollinator network, where
interactions are structured by an individual trait. We found continuous limits
of the microscopic system and finally showed that it could be approximated by
PDEs. We finally studied the dynamics, stationary state and stability of this
continuous limits. Our approach allows to address several limits exposed
before: the relationships between the individuals, species and community
scales are explicit; interactions variability within and between-species is
taken into account; the interaction graph is based on individuals’ traits;
PDEs approximations allows an explicit and analytical study of the property of
a community. We showed in particular that a nested plant-pollinator network is
expected to collapse, a phenomenon already observed by previous works, but to
our knowledge for the first time formally demonstrated under general
conditions. Our approach can thus provide a new and original theoretical
framework for ecologists to address long-standing questions.

One of the actual limit of our model is that interactions take place through a
“mass-action model", in particular with the use of the resources (1.5)-(1.6)
which are already at a macroscopic scale. Returning to an event-based
modelling for establishing these functional responses (e.g. [5]) and showing
how interactions at the individual level would translate into an interaction
graph and into the dynamics of the whole community is an open question.

## Code

The code for simulations is available here:
`https://gitlab.com/thoma.rey/PlantPollinatorsNetwork/`

##  Appendix A Proofs

###  A.1 Bounds on the microscopic representation

Here, we give a sketch of proof of Proposition 1.7. The proof follows from
usual stochastic calculus with Poisson point processes, as developed in [25]
for example.

###### Proof of Proposition 1.7.

Step 1: let us first prove some moment estimates, including (1.12). First,
notice that the process can be stochastically bounded by a process where only
births occur at the maximal birth rate, denoted by
$M^{P}:=\sup_{r\in\mathbb{R}^{+}}b^{P}(r)$:

$\displaystyle\langle\mathbf{P}^{K,n}_{t},\mathds{1}\rangle$
$\displaystyle=\frac{1}{nK}\sum_{i=1}^{n}P^{K,i}_{t}$
$\displaystyle\leq\langle\tilde{\mathbf{P}}^{K,n}_{t},\mathds{1}\rangle=\langle\mathbf{P}^{K,n}_{0},\mathds{1}\rangle+\int_{0}^{t}\int_{E}\frac{1}{nK}\sum_{i=1}^{n}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{i=k}{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}_{\theta\leq
M^{P}\tilde{P}^{K,i}_{s-}}Q^{P}_{B}(ds,dk,d\theta).$

Since

$\mathbb{E}\Big{(}\sup_{s\leq
t\wedge\tau^{K}_{N}}\langle\mathbf{P}^{K,n}_{s},\mathds{1}\rangle\Big{)}\leq
1+\mathbb{E}\Big{(}\sup_{s\leq
t\wedge\tau^{K}_{N}}\langle\mathbf{P}^{K,n}_{s},\mathds{1}\rangle^{2}\Big{)}\leq
1+\mathbb{E}\Big{(}\sup_{s\leq
t\wedge\tau^{K}_{N}}\langle\tilde{\mathbf{P}}^{K,n}_{s},\mathds{1}\rangle^{2}\Big{)},$

we only have to find a bound for the second moment of
$\langle\tilde{\mathbf{P}}^{K,n}_{s},\mathds{1}\rangle$.
For a constant $N>0$, we introduce the stopping time

$\tau^{K}_{N}=\inf\Big{\\{}t\geq 0,\
\langle\tilde{\mathbf{P}}^{K,n}_{t},\mathds{1}\rangle\geq N\mbox{ or
}\langle\tilde{\mathbf{A}}^{K,m}_{t},\mathds{1}\rangle\geq N\Big{\\}}.$

Using Itô’s formula for jump processes (e.g. [34, Th.5.1 P.66]), we have:

$\displaystyle\langle$
$\displaystyle\tilde{\mathbf{P}}^{K,n}_{t},\mathds{1}\rangle^{2}$
$\displaystyle=$
$\displaystyle\langle\mathbf{P}^{K,n}_{0},\mathds{1}\rangle^{2}+\int_{0}^{t}\int_{E}\sum_{i=1}^{n}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{k=i}\Big{(}\big{(}\langle\tilde{\mathbf{P}}^{K,n}_{s-},\mathds{1}\rangle+\frac{1}{nK}\big{)}^{2}-\langle\tilde{\mathbf{P}}^{K,n}_{s-},\mathds{1}\rangle^{2}\Big{)}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\theta\leq
M^{P}\tilde{P}^{K,i}_{s-}}Q^{P}_{B}(ds,dk,d\theta)$ $\displaystyle=$
$\displaystyle\langle\mathbf{P}^{K,n}_{0},\mathds{1}\rangle^{2}+\int_{0}^{t}\int_{E}\sum_{i=1}^{n}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{k=i}\Big{(}\frac{2}{nK}\sup_{u\leq
s{-}}\langle\tilde{\mathbf{P}}^{K,n}_{u},\mathds{1}\rangle+\frac{1}{n^{2}K^{2}}\Big{)}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\theta\leq
M^{P}\tilde{P}^{K,i}_{s-}}Q^{P}_{B}(ds,dk,d\theta).$

Since the right hand side is an increasing process, we can again replace the
left hand side by $\sup_{s\leq
t\wedge\tau_{N}^{K}}\langle\tilde{\mathbf{P}}^{K,n}_{s},\mathds{1}\rangle$.
Then, taking the expectation gives:

$\displaystyle\mathbb{E}\Big{(}\sup_{s\leq
t\wedge\tau^{K}_{N}}\langle\tilde{\mathbf{P}}^{K,n}_{s},\mathds{1}\rangle^{2}\Big{)}\leq$
$\displaystyle\mathbb{E}\big{(}\langle\mathbf{P}^{K,n}_{0},\mathds{1}\rangle^{2}\big{)}+\mathbb{E}\Big{[}\int_{0}^{t\wedge\tau^{K}_{N}}\sum_{i=1}^{n}\Big{(}\frac{2}{nK}\sup_{u\leq
s}\langle\tilde{\mathbf{P}}^{K,n}_{u},\mathds{1}\rangle+\frac{1}{n^{2}K^{2}}\Big{)}M^{P}\tilde{P}^{K,i}_{s}\
ds\Big{]}$ $\displaystyle\leq$
$\displaystyle\mathbb{E}\big{(}\langle\mathbf{P}^{K,n}_{0},\mathds{1}\rangle^{2}\big{)}+\int_{0}^{t}M^{P}\Big{[}(2+\frac{1}{nK})\mathbb{E}\Big{(}1+\sup_{u\leq
s\wedge\tau^{K}_{N}}\langle\tilde{\mathbf{P}}^{K,n}_{u},\mathds{1}\rangle^{2}\Big{)}\Big{]}\
ds.$

Using Gronwall’s lemma, we deduce

$\displaystyle\mathbb{E}\Big{(}\sup_{s\leq
t\wedge\tau^{K}_{N}}\langle\tilde{\mathbf{P}}^{K,n}_{s},\mathds{1}\rangle^{2}\Big{)}\leq$
$\displaystyle\Big{[}\mathbb{E}\big{(}\langle\mathbf{P}^{K,n}_{0},\mathds{1}\rangle^{2}\big{)}+3M^{P}\Big{]}e^{3M^{P}t}<+\infty.$

Since the bound does not depend on $N$, we make $N$ grow to $\infty$ and
obtain

$\displaystyle\mathbb{E}\Big{(}\sup_{s\leq
t}\langle\tilde{\mathbf{P}}^{K,n}_{s},\mathds{1}\rangle^{2}\Big{)}\leq$
$\displaystyle\Big{[}\mathbb{E}\big{(}\langle\mathbf{P}^{K,n}_{0},\mathds{1}\rangle^{2}\big{)}+3M^{P}\Big{]}e^{3M^{P}t}<+\infty.$

Similarly,

$\mathbb{E}\Big{(}\sup_{s\leq
t}\langle\tilde{\mathbf{A}}^{K,m}_{s},\mathds{1}\rangle^{2}\Big{)}\leq\Big{[}\mathbb{E}\big{(}\langle\mathbf{A}^{K,m}_{0},\mathds{1}\rangle^{2}\big{)}+3M^{A}\Big{]}e^{3M^{A}t}<+\infty.$
(A.1)

Notice that the upper bound that we obtain do not depend on $K$. Once these
estimates have been obtained, existence and uniqueness of a strong solution to
the stochastic differential equation is given by [34, Th. IV.9.1].

Step 2: Let us now consider a measurable real function $f$ on $[0,1]$. Using
Itô’s formula for SDEs with jumps (see e.g. [34, p. 66-67]) and (1.10)
provides (1.13) in which:

$\displaystyle M^{K,i}_{t}=$
$\displaystyle\int_{0}^{t}\int_{E}{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}_{i=k}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\theta\leq b^{P}\big{(}\sum_{j\sim
i}R^{A,K,ij}_{s-}\big{)}P^{K,i}_{s-}}Q^{P}_{B}(ds,dk,d\theta)-\int_{0}^{t}b^{P}\big{(}\sum_{j\sim
i}R^{K,ij}_{s}\big{)}P^{K,i}_{s}ds$ $\displaystyle-$
$\displaystyle\int_{0}^{t}\int_{E}{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}_{i=k}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\theta\leq\big{[}d^{P}\big{(}\sum_{j\sim
i}R^{A,K,ij}_{s-}\big{)}-k\star\mathbf{P}^{K,n}_{s-}(x^{i})\big{]}P^{K,i}_{s-}}Q^{P}_{D}(ds,dk,d\theta)$
$\displaystyle+$ $\displaystyle\int_{0}^{t}\big{[}d^{P}\big{(}\sum_{j\sim
i}R^{A,K,ij}_{s}\big{)}-k\star\mathbf{P}^{K,n}_{s}(x^{i})\big{]}P^{K,i}_{s}ds$
(A.2)

is a square integrable martingale with bracket given by (1.14). A similar
computation can be done for $\mathbf{A}^{K,m}$ to obtain (1.15). ∎

###  A.2 Large number of species limit

###### Proof of Proposition 3.2.

Recall that for $n,m\geq 1$, we have:

$\displaystyle\int_{0}^{1}f(x)d\widetilde{\mathbf{P}}^{n,m}_{t}(x)=$
$\displaystyle\int_{0}^{1}f(x)d\widetilde{\mathbf{P}}^{n,m}_{0}(x)dx$ (A.3)
$\displaystyle+\int_{0}^{t}\frac{1}{n}\sum_{i=1}^{n}f(x^{i})\Big{[}g^{P}\Big{(}\frac{1}{m}\sum_{j=1}^{m}mc^{n,m}_{i,j}G_{ij}\widetilde{A}^{j}_{s}\Big{)}-\frac{1}{n}\sum_{\ell=1}^{n}k(x^{i},x^{\ell})\widetilde{P}^{\ell}_{s}\Big{]}\widetilde{P}^{i}_{s}ds,$
$\displaystyle\int_{0}^{1}f(y)d{\widetilde{\mathbf{A}}}^{n,m}_{t}(y)=$
$\displaystyle\int_{0}^{1}f(y)d\widetilde{\mathbf{A}}^{n,m}_{0}(y)dy$
$\displaystyle+\int_{0}^{t}\frac{1}{m}\sum_{j=1}^{m}f(y^{j})\Big{[}g^{A}\Big{(}\frac{1}{n}\sum_{i=1}^{n}nc^{n,m}_{i,j}G_{ij}\widetilde{P}^{i}_{s}\Big{)}-\frac{1}{m}\sum_{\ell=1}^{m}h(y^{j},y^{\ell})\widetilde{A}^{\ell}_{s}\Big{]}\widetilde{A}^{j}_{s}\
ds.$

where conditionally on $(x^{i},y^{j})_{1\leq i\leq n,1\leq j\leq m}$,
$G_{ij}\sim Bern(\phi(x^{i},y^{j}))$ and are independent random variables, and
$c^{n,m}_{i,j}$ are independent random variables satisfying Equations (1.2)
and (1.3).

Step 1: Processes are bounded. Since $g^{P}$ and $g^{A}$ are bounded from
above respectively by $M_{P}$ and $M_{A}$ and since the competition terms for
the pollinators are non-positive, choosing $f\equiv 1$ in the previous
equation (A.3), we find that a.s.,

$\displaystyle\langle\widetilde{\mathbf{P}}^{n,m}_{t},1\rangle\leq\langle\widetilde{\mathbf{P}}^{n,m}_{0},1\rangle
e^{M_{P}t}\quad\mbox{ and
}\quad\langle\widetilde{\mathbf{A}}^{n,m}_{t},1\rangle\leq\langle\widetilde{\mathbf{A}}^{n,m}_{0},1\rangle
e^{M_{A}t},$ (A.4)

for all $n$ and $m$. Moreover, from Assumption (3.2), we have that

$\lim_{n,m\to+\infty}\langle\widetilde{\mathbf{P}}^{n,m}_{0},1\rangle=\int_{[0,1]}\bar{p}_{0}(x)dx\quad\text{and}\quad\lim_{n,m\to+\infty}\langle\widetilde{\mathbf{A}}^{n,m}_{0},1\rangle=\int_{[0,1]}\bar{a}_{0}(x)dx\quad\text{a.s.}.$

In other words,
$\sup_{n,m}(\langle\widetilde{\mathbf{P}}^{n,m}_{0},1\rangle+\langle\widetilde{\mathbf{A}}^{n,m}_{0},1\rangle)$
is bounded by a finite constant a.s. and for $T>0$, we define

$C(T)=\sup_{n,m}\left(\langle\widetilde{\mathbf{P}}^{n,m}_{0},1\rangle
e^{M_{P}T}+\langle\widetilde{\mathbf{A}}^{n,m}_{0},1\rangle
e^{M_{A}T}\right)<+\infty\text{ a.s.}.$ (A.5)

Thus, the processes
$\big{(}\widetilde{\mathbf{P}}^{n,m}_{t}(dx),\widetilde{\mathbf{A}}_{t}^{n,m}(dy)\big{)}_{t\in[0,T]}$
take their values in $(\mathcal{M}_{\leq C(T)}([0,1]))^{2}$ which is a compact
set.

Step 2: Relative compactness in $\mathcal{C}\Big{(}[0,T],\mathcal{M}_{\leq
C(T)}([0,1])^{2}\Big{)}$

Following [11, Theorem 7.3] and proof of Theorem 5.3 in [25], it is sufficient
to prove that

  * (i)

$\sup_{n,m\in\mathbb{N}}\left(\langle\widetilde{\mathbf{P}}^{n,m}_{0},1\rangle+\langle\widetilde{\mathbf{A}}^{n,m}_{0},1\rangle\right)$
is bounded almost surely, and

  * (ii)

for any continuous function $f$ on $[0,1]$, and for any $\varepsilon>0$,

$\displaystyle\lim_{\delta\to
0}\limsup_{n,m\in\mathbb{N}}\mathbb{P}\left(\sup_{|t-s|\leq\delta}|\langle\widetilde{\mathbf{P}}^{n,m}_{s},f\rangle-\langle\widetilde{\mathbf{P}}^{n,m}_{t},f\rangle|\geq\varepsilon\right)=0$
$\displaystyle\lim_{\delta\to
0}\limsup_{n,m\in\mathbb{N}}\mathbb{P}\left(\sup_{|t-s|\leq\delta}|\langle\widetilde{\mathbf{A}}^{n,m}_{s},f\rangle-\langle\widetilde{\mathbf{A}}^{n,m}_{t},f\rangle|\geq\varepsilon\right)=0$

Point (i) follows from Equation (A.5). Concerning point (ii), notice first
that

$\mathbb{E}[\widetilde{R}^{A,n,m}]=\mathbb{E}\left[\frac{1}{m}\sum_{j=1}^{m}mc^{n,m}_{i,j}G_{ij}\widetilde{A}^{j}_{s}\right]\leq\frac{1}{m}\sum_{j=1}^{m}\mathbb{E}[mc^{n,m}_{i,j}]\
C(T)\leq\max_{n,m\in\mathbb{N}^{*}}\|mc^{n,m}\|_{\infty}\ C(T)<\infty,$

which is finite from Assumption 3.1. Thus, there exists $\tilde{C}>0$ such
that $\mathbb{P}(\widetilde{R}^{A,n,m}>\tilde{C})$ is sufficiently small. Then
$g^{P}$ is bounded on $[0,\tilde{C}]$ by some constant $\tilde{M}$ (since it
is a Lypschitz function). Thus, on $\\{\tilde{R}^{A,n,m}\leq\tilde{C}\\}$, for
any positive continuous function $f$ on $[0,1]$ and any $0\leq s\leq t\leq T$,

$\displaystyle|\langle\widetilde{\mathbf{P}}^{n,m}_{t},f\rangle-\langle\widetilde{\mathbf{P}}^{n,m}_{s},f\rangle|$
$\displaystyle\leq\int_{s}^{t}(\tilde{M}+\|k\|_{\infty}C(T))\|f\|_{\infty}C(T)du.$

Same computations can be done for the sequence of processes
$(\widetilde{\mathbf{A}}^{n,m})_{m\in\mathbb{N}^{*}}$. And these are
sufficient to conclude point (ii).

Step 3: Density of the limiting values. Let us then prove that every limiting
values of
$\big{(}\widetilde{\mathbf{P}}^{n,m},\widetilde{\mathbf{A}}^{n,m}\big{)}$,
denoted here by $(\bar{\mathbf{P}},\bar{\mathbf{A}})$, is a process whose time
marginals at $t>0$ admit densities $\bar{p}_{t}$ and $\bar{a}_{t}$ with
respect to the Lebesgue measure on $[0,1]$. To achieve this, we dominate the
measures $\widetilde{\mathbf{P}}^{n,m}_{t}(dx)$ and
$\widetilde{\mathbf{A}}_{t}^{n,m}(dy)$ by measures with densities. For all
positive continuous function $f$ on $[0,1]$, with arguments similar to these
of (A.4), we prove that

$\displaystyle\langle\bar{\mathbf{P}}_{t},f\rangle\leq\langle\bar{\mathbf{P}}_{0},f\rangle
e^{M_{P}t}=\int_{0}^{1}f(x)\big{(}p_{0}(x)e^{M_{P}t}\big{)}dx,\quad\text{a.s.}$
$\displaystyle\langle\bar{\mathbf{A}}_{t},f\rangle\leq\langle\bar{\mathbf{A}}_{0},f\rangle
e^{M_{A}t}=\int_{0}^{1}f(y)\big{(}a_{0}(y)e^{M_{A}t}\big{)}dy,\quad\text{a.s.}$

which is sufficient to deduce the existence of densities with respect to
Lebesgue measure.

Step 4: Equation satisfied by the limiting value. Let us fix $t>0$ and $f$ a
continuous positive function on $[0,1]$. We define the map $\psi$ on
$\mathcal{C}([0,1],\mathcal{M}_{\leq C(T)}([0,1])^{2})$, endowed with the
uniform convergence, defined by

$\psi(\mathbf{P},\mathbf{A})=\begin{pmatrix}\langle\mathbf{P}_{t},f\rangle-\langle\mathbf{P}_{0},f\rangle-\int_{0}^{t}\int_{0}^{1}f(x)\left[g^{P}\left((c\phi)\star\mathbf{A}_{s}(x)\right)-k\star\mathbf{P}_{s}(x)\right]d\mathbf{P}_{s}(x)ds\\\\[5.0pt]
\langle\mathbf{A}_{t},f\rangle-\langle\mathbf{A}_{0},f\rangle-\int_{0}^{t}\int_{0}^{1}f(x)\left[g^{A}\left((c\phi)\star\mathbf{P}_{s}(x)\right)-h\star\mathbf{A}_{s}(x)\right]d\mathbf{A}_{s}(x)ds\end{pmatrix}.$

The continuity of map $\psi$ is straightforward once noticing that the weak
convergence on $\mathcal{M}_{\leq C(T)}^{2}([0,1])$ is equivalent to
convergence of measures in the Kantorovich–Rubinstein distance (see for
example Theorem 6.9 of [48])

$\mathcal{W}_{1}(\nu,\mu)=\sup\left\\{\int_{0}^{1}f(x)d(\nu-\mu)(x)\Big{|}f:[0,1]\mapsto\mathbb{R}^{+}\text{
Lipschitz },Lip(f)\leq 1\right\\},$

since $[0,1]$ is a compact space. The function $\psi$ is also bounded since
all measures are bounded by $C(T)$.
Then, as $(\bar{\mathbf{P}},\bar{\mathbf{A}})$ is a weak limiting value of
$\big{(}\widetilde{\mathbf{P}}^{n,m},\widetilde{\mathbf{A}}^{n,m}\big{)}$, we
obtain by definition of the weak convergence that

$\mathbb{E}\left[\|\psi(\widetilde{\mathbf{P}}^{n,m},\widetilde{\mathbf{A}}^{n,m})\|_{1}\right]\underset{n,m\to+\infty}{\longrightarrow}\mathbb{E}\left[\|\psi(\bar{\mathbf{P}},\bar{\mathbf{A}})\|_{1}\right].$
(A.6)

On the other hand,

$\displaystyle\mathbb{E}$
$\displaystyle\left[\|\psi(\widetilde{\mathbf{P}}^{n,m},\widetilde{\mathbf{A}}^{n,m})\|_{1}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\left|\int_{0}^{t}\int_{0}^{1}f(x)\bigg{(}g^{P}\Big{(}\frac{1}{m}\sum_{j=1}^{m}G_{ij}mc^{n,m}_{i,j}\widetilde{A}^{j}_{s}\Big{)}-g^{P}\left((c\phi)\star\widetilde{\mathbf{A}}^{n,m}\right)\bigg{)}d\widetilde{\mathbf{P}}^{n,m}_{s}(x)ds\right|\right]$
(A.7)
$\displaystyle+\mathbb{E}\left[\left|\int_{0}^{t}\int_{0}^{1}f(x)\bigg{(}g^{A}\Big{(}\frac{1}{n}\sum_{j=1}^{n}G_{ij}nc^{n,m}_{i,j}\widetilde{P}^{j}_{s}\Big{)}-g^{A}\left((c\phi)\star\widetilde{\mathbf{P}}^{n,m}\right)\bigg{)}d\widetilde{\mathbf{A}}^{n,m}_{s}(x)ds\right|\right].$
(A.8)

Using Assumptions 1.5, the first term (A.7) of the r.h.s can be bounded from
above by

$\displaystyle\eqref{term1}$
$\displaystyle\leq\int_{0}^{t}\frac{C}{n}\sum_{i=1}^{n}\left[\mathbb{E}\left|\frac{1}{m}\sum_{j=1}^{m}\widetilde{A}^{j}_{s}\Big{(}G_{ij}\big{(}mc^{n,m}_{i,j}-c(x^{i},y^{j})\big{)}+c(x^{i},y^{j})\big{(}G_{ij}-\phi(x^{i},y^{j})\big{)}\Big{)}\right|\right]ds$
$\displaystyle\leq C\Bigg{[}TC(T)\|mc^{n,m}-c\|_{\infty}$
$\displaystyle\qquad+\int_{0}^{t}\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left|\frac{1}{m}\sum_{j=1}^{m}\widetilde{A}^{j}_{s}G_{ij}m\big{(}c^{n,m}_{i,j}-c^{n,m}(x^{i},y^{j})\big{)}+\widetilde{A}^{j}_{s}c(x^{i},y^{j})\big{(}G_{ij}-\phi(x^{i},y^{j})\big{)}\right|\Bigg{]},$

with $C:=\|f\|_{\infty}C(T)L^{P}$ and $L^{P}$ is the Lipschitz constant of
function $g^{P}$. Then $\|mc^{n,m}-c\|_{\infty}$ converges to $0$ when $n,m$
grow to $\infty$, and it will be sufficiently small for $n,m$ large enough.
Finally, from Cauchy-Schwarz inequality, Assumptions 3.1 and 1.1, and denoting
by $Var_{x,y}(A)=Var(A|(x^{i},y^{j})_{i\leq n,j\leq m})$, we find for $n,m$
sufficiently large,

$\displaystyle\eqref{term1}$
$\displaystyle\leq\varepsilon+\int_{0}^{t}\frac{C}{n}\sum_{i=1}^{n}\mathbb{E}\left[\text{Var}_{x,y}\left(\frac{1}{m}\sum_{j=1}^{m}\widetilde{A}^{j}_{s}G_{ij}mc^{n,m}_{i,j}\right)\right]^{1/2}+\mathbb{E}\left[\text{Var}_{x,y}\left(\frac{1}{m}\sum_{j=1}^{m}\widetilde{A}^{j}_{s}c(x^{i},y^{j})G_{ij}\right)\right]^{1/2}$
$\displaystyle\leq\varepsilon+\int_{0}^{t}\frac{C}{n}\sum_{i=1}^{n}\mathbb{E}\left[\frac{C(T)^{2}}{m^{2}}\sum_{j=1}^{m}V_{max}\right]^{1/2}+\mathbb{E}\left[\frac{C(T)^{2}\|c\|_{\infty}^{2}}{m^{2}}\sum_{j=1}^{m}\phi(x^{i},y^{j})(1-\phi(x^{i},y^{j}))\right]^{1/2}$
$\displaystyle\leq\varepsilon+CT\frac{C(T)(V_{max}^{1/2}+\|c\|_{\infty}\|\phi(1-\phi)\|_{\infty}^{1/2})}{m^{1/2}}$
$\displaystyle\leq 2\varepsilon\text{ for }n,m\text{ sufficiently large}.$

The second term (A.8) can be treated similarly. This finally proves that
$\mathbb{E}\left[\|\psi(\widetilde{\mathbf{P}}^{n,m},\widetilde{\mathbf{A}}^{n,m})\|_{1}\right]$
converges to $0$ when $n,m$ go to $\infty$. In addition with (A.6), we deduce
that $\psi(\bar{\mathbf{P}},\bar{\mathbf{A}})=(0,0)$ a.s. In other words, the
limiting value satisfies Equation (3.3) a.s.

Step 5: Uniqueness of the solution of (3.3). It remains to prove the
uniqueness of the solution to (3.3) to conclude. To this aim, let us take
$(\mathbf{P}^{1},\mathbf{A}^{1})$ and $(\mathbf{P}^{2},\mathbf{A}^{2})$ two
deterministic measures, weak solutions to (3.3) and with identical initial
conditions. Let us denote by $\rho$ the Radon metric between two measures
$\nu$ and $\mu$ on $\mathcal{M}_{F}([0,1])^{2}$:

$\rho(\nu,\mu)=\sup\left\\{\int_{0}^{1}f(x)d(\nu-\mu)(x)\Big{|}f:[0,1]\mapsto[-1,1]\text{
continuous}\right\\}.$

With straightforward computations, for any continuous function
$f:[0,1]\mapsto[-1,1]$, we can find some finite constants $C_{1},C_{2}$,
independent from $t$ and $f$, such that

$\displaystyle|\langle\mathbf{P}^{1}_{t}-\mathbf{P}^{2}_{t},f\rangle|$
$\displaystyle=\left|\int_{0}^{t}\int_{0}^{1}f(x)\left[g^{P}(c\phi\star\mathbf{A}^{1}_{s}(x))-k\star\mathbf{P}^{1}_{s}\right]d\mathbf{P}^{1}_{s}(x)ds\right.$
$\displaystyle\qquad\left.-\int_{0}^{t}\int_{0}^{1}f(x)\left[g^{P}(c\phi\star\mathbf{A}^{2}_{s}(x))-k\star\mathbf{P}^{2}_{s}\right]d\mathbf{P}^{2}_{s}(x)ds\right|$
$\displaystyle\leq
C_{1}\int_{0}^{t}\rho(\mathbf{P}^{1}_{s},\mathbf{P}^{2}_{s})ds+\int_{0}^{t}\int_{0}^{1}\left(L^{P}\left|c\phi\star(\mathbf{A}^{1}_{s}-\mathbf{A}^{2}_{s}(x))\right|+\left|k\star(\mathbf{P}^{1}_{s}-\mathbf{P}^{2}_{s})(x)\right|\right)d\mathbf{P}^{2}_{s}(x)ds$
$\displaystyle\leq
C_{2}\int_{0}^{t}(\rho(\mathbf{P}^{1}_{s},\mathbf{P}^{2}_{s})+\rho(\mathbf{A}^{1}_{s},\mathbf{A}^{2}_{s}))ds,$

since $x\mapsto
f(x)(g^{P}(c\phi\star\mathbf{A}^{1}_{s}(x))-k\star\mathbf{P}^{1}_{s}(x))$ are
continuous bounded functions on $[0,1]$ for all $s\in[0,T]$, and $k$ and
$c\phi$ are continuous bounded functions on $[0,1]^{2}$. Same computations can
be done with $\mathbf{A}^{1}$ and $\mathbf{A}^{2}$, then taking the supremum
over all $f$, we find a constant $C_{3}>0$ such that

$\rho(\mathbf{P}^{1}_{t},\mathbf{P}^{2}_{t})+\rho(\mathbf{A}^{1}_{t},\mathbf{A}^{2}_{t})\leq
C_{2}\int_{0}^{t}(\rho(\mathbf{P}^{1}_{s},\mathbf{P}^{2}_{s})+\rho(\mathbf{A}^{1}_{s},\mathbf{A}^{2}_{s}))ds.$

We conclude with Gronwall inequality that
$\rho(\mathbf{P}^{1}_{t},\mathbf{P}^{2}_{t})+\rho(\mathbf{A}^{1}_{t},\mathbf{A}^{2}_{t})=0$
for all $t\in[0,T]$. In other words, $(\mathbf{P}^{1},\mathbf{A}^{1})$ and
$(\mathbf{P}^{2},\mathbf{A}^{2})$ are identical, and the solution to (3.3) is
unique.

This ends the proof. ∎

## References

  * [1] E. Abbe.  Community detection and stochastic block models: recent development.  Journal of Machine Learning Research, 18(177):1–86, 2018.
  * [2] I. Akjouj, M. Barbier, M. Clénet, W. Hachem, M. Ma\̈mathbf{i}da, F. Massol, J. Najim, and V.C. Tran.  Complex systems in ecology: A guided tour with the lotka-volterra model.  review in progress, 2021.
  * [3] G. Ben Arous, Y.V. Fyodorov, and B.A. Khoruzhenko.  Counting equilibria of large complex systems by instability index.  Proceedings of the National Academy of Sciences, 118(34):e2023719118, August 2021.
  * [4] Blanca Arroyo-Correa, Ignasi Bartomeus, and Pedro Jordano.  Individual-based plant–pollinator networks are structured by phenotypic and microsite plant traits.  Journal of Ecology, 109(8):2832–2844, 2021.
  * [5] V. Bansaye, S. Billiard, and J.R. Chazottes.  Rejuvenating functional responses avec renewal theory.  Journal of the Royal Society Interface, 15(146):20180239, 2018.
  * [6] V. Bansaye and S. Méléard.  Stochastic models for structured populations. Scaling limits and long time behavior.  MBI Lecture Series 1.4. Springer, 2015.
  * [7] M. Barbier, J.F. Arnoldi, G. Bunin, and M. Loreau.  Generic assembly patterns in complex ecological communities.  Proceedings of the National Academy of Sciences, 115(9):2156–2161, 2018.
  * [8] J. Bascompte, P. Jordano, C. Melian, and J.M. Olesen.  The nested assembly of plant-animal mutualistic networks.  Proceedings of the National Academy of Sciences, 100(16):9383–9387, 2003.
  * [9] Nicola Bellomo, Ahmed Elaiw, Abdullah M Althiabi, and Mohammed Ali Alghamdi.  On the interplay between mathematics and biology: Hallmarks toward a new systems biology.  Physics of life reviews, 12:44–64, 2015.
  * [10] Cecilia Berardo, Stefan Geritz, Mats Gyllenberg, and Gaël Raoul.  Interactions between different predator–prey states: a method for the derivation of the functional and numerical response.  Journal of Mathematical Biology, 80:2431–2468, 2020.
  * [11] P. Billingsley.  Convergence of Probability Measures.  Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., New York, second edition edition, 1999.
  * [12] B. Bollobás.  Random graphs.  Cambridge University Press, 2 edition, 2001.
  * [13] C. Borgs, J. Chayes, L. Lovász, V. Sós, and K. Vesztergombi.  Limits of randomly grown graph sequences.  European Journal of Combinatorics, 32(7):985–999, 2011.
  * [14] G. Bunin.  Ecological communities with Lotka-Volterra dynamics.  Physical Review E, 95(4):042414, 2017.
  * [15] S.A. Chamberlain, R.V. Cartar, A.C. Worley, S.J. Semmler, G. Gielens, S. Elwell, M.E. Evans, J.C. Vamosi, and E. Elle.  Traits and phylogenetic history contribute to network structure across canadian plant-pollinator communities.  Oecologia, 176:545–556, 2014.
  * [16] N. Champagnat, R. Ferrière, and S. Méléard.  Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models via timescale separation.  Theoretical Population Biology, 69:297–321, 2006.
  * [17] J.E. Cohen and C.M. Newman.  When will a large complex system be stable?  Journal of Theoretical Biology, 113:153–156, 1985.
  * [18] Eva Delmas, Mathilde Besson, Marie-Hélène Brice, Laura A. Burkle, Giulio V. Dalla Riva, Marie-Josée Fortin, Dominique Gravel, Paulo R. Guimarães Jr., David H. Hembry, Erica A. Newman, Jens M. Olesen, Mathias M. Pires, Justin D. Yeakel, and Timothée Poisot.  Analysing ecological networks of species interactions.  Biological Reviews, 94(1):16–36, 2019.
  * [19] R. Van der Hofstad.  Random Graphs and Complex Networks, volume 1 of Cambridge Series in Statistical and Probabilistic Mathematics.  Cambridge University Press, Cambridge, 2017.
  * [20] Laurent Desvillettes, Pierre Emmanuel Jabin, Stéphane Mischler, and Gaël Raoul.  On selection dynamics for continuous structured populations.  Communications in Mathematical Sciences, 6(3):729–747, 2008.
  * [21] R. Durrett.  Random graph dynamics.  Cambridge University Press, New York, 2007.
  * [22] S.N. Ethier and T.G. Kurtz.  Markov Processes, Characterization and Convergence.  John Wiley & Sons, New York, 1986.
  * [23] R. Ferrière and V.C. Tran.  Stochastic and deterministic models for age-structured populations with genetically variable traits.  ESAIM: Proceedings, 27:289–310, 2009.
  * [24] C. Fontaine, P.R. Guimaraes Jr., S. Kéfi, and E. Thébault.  The ecological and evolutionary implications of merging different types of networks.  Ecology Letters, 14(11):1170–1181, 2011.
  * [25] N. Fournier and S. Méléard.  A microscopic probabilistic description of a locally regulated population and macroscopic approximations.  Ann. Appl. Probab., 14(4):1880–1919, 2004.
  * [26] Y.V. Fyodorov and B.A. Khoruzhenko.  Nonlinear analogue of the May-Wigner instability transition.  Proceedings of the National Academy of Sciences, 113(25):6827–6832, 2016.
  * [27] D.T. Gillespie.  A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.  Journal of Computational Physics, 22(4):403–434, 1976.
  * [28] D.T. Gillespie.  Exact stochastic simulation of coupled chemical reactions.  Journal of Physical Chemistry, 81(25):2340–2361, 1977.
  * [29] Paulo R. Guimarães.  The structure of ecological networks across levels of organization.  Annual Review of Ecology, Evolution, and Systematics, 51(1):433–460, 2020.
  * [30] J.N. Holland and J.L. Bronstein.  Mutualism.  In Sven Erik Jørgensen and Brian D. Fath, editors, Encyclopedia of Ecology, pages 2485–2491. Academic Press, Oxford, 2008.
  * [31] JN Holland and Judith L Bronstein.  Mutualism.  In Encyclopedia of Ecology, Five-Volume Set, pages 2485–2491. Elsevier Inc., 2008.
  * [32] J.N. Holland and D.L DeAngelis.  A consumer–resource approach to the density-dependent population dynamics of mutualism.  Ecology, 91:1286–1295, 2010.
  * [33] P. Holland, K. Laskey, and S. Leinhardt.  Stochastic blockmodels: some first steps.  Social Networks, 5:109–137, 1983.
  * [34] N. Ikeda and S. Watanabe.  Stochastic Differential Equations and Diffusion Processes, volume 24.  North-Holland Publishing Company, 1989.  Second Edition.
  * [35] Pierre-Emmanuel Jabin and Gaël Raoul.  On selection dynamics for competitive interactions.  Journal of Mathematical Biology, 63(3):493–517, 2011.
  * [36] A. Joffe and M. Métivier.  Weak convergence of sequences of semimartingales with applications to multitype branching processes.  Advances in Applied Probability, 18:20–65, 1986.
  * [37] E. Knop, L. Zoller, R. Ryser, C. Gerpe, M. Hörler, and C. Fontaine.  Artificial light at night as a new threat to pollination.  Nature, 548(7666):206–209, 2017.
  * [38] J. Jelle Lever, Egbert H. van Nes, Marten Scheffer, and Jordi Bascompte.  The sudden collapse of pollinator communities.  Ecology letters, 17(3):350–359, 2014.
  * [39] L. Lovász.  Large networks and graph limits.  (60), 2012.
  * [40] R.M. May.  Will a large complex system be stable?  Nature, 238:413–414, 1972.
  * [41] R.M. May.  Qualitative stability in model ecosystems.  Ecology, 54(3):638–641, 1973.
  * [42] Camille Pouchol and Emmanuel Trélat.  Global stability with selection in integro-differential lotka-volterra systems modelling trait-structured populations.  Journal of Biological Dynamics, 12(1):872–893, 2018.
  * [43] Jonathan Roughgarden.  Theory of Population Genetics and Evolutionary Ecology: An Introduction.  Macmillan, 1979.
  * [44] Y. Takeuchi.  Global dynamical properties of Lotka-Volterra systems.  World Scientific, 1996.
  * [45] S. Tang and S. Allesina.  Reactivity and stability of large ecosystems.  Frontiers in Ecology and Evolution, 2(21), 2014.
  * [46] E. Thébault and C. Fontaine.  Stability of ecological communities and the architecture of mutualistic and trophic networks.  Science, 329(5993):853–856, 2010.
  * [47] S. Villalobos, J.M. Sevenello-Montagner, and J.C. Vamosi.  Specialization in plant–pollinator networks: insights from local-scale interactions in glenbow ranch provincial park in alberta, canada.  BMC Ecology, 19:34, 2019.
  * [48] C. Villani.  Optimal transport: old and new.  Springer, 338, 2009.
  * [49] Stella Watts, Carsten F. Dormann, Ana M. Mart\́mathbf{i}n González, and Jeff Ollerton.  The influence of floral traits on specialization and modularity of plant–pollinator networks in a biodiversity hotspot in the Peruvian Andes.  Annals of Botany, 118(3):415–429, 07 2016.