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

HAL
open science

# Continuous limits of large plant-pollinator random networks and some applications

Sylvain Billiard, Hélène Leman, Thomas Rey, Viet-Chi Tran

► To cite this version:

Sylvain Billiard, Hélène Leman, Thomas Rey, Viet-Chi Tran. Continuous limits of large plant-pollinator random networks and some applications. MathematicS In Action, In press. hal-03525607v2

HAL Id: hal-03525607

https://cnrs.hal.science/hal-03525607v2

Submitted on 23 Jun 2022

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

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

HAL
open science

# Continuous limits of large plant-pollinator random networks and some applications

Sylvain Billiard, Hélène Leman, Thomas Rey, Viet-Chi Tran

► To cite this version:

Sylvain Billiard, Hélène Leman, Thomas Rey, Viet-Chi Tran. Continuous limits of large plant-pollinator random networks and some applications. MathematicS In Action, 2025. hal-03525607v3

HAL Id: hal-03525607

https://cnrs.hal.science/hal-03525607v3

Submitted on 26 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Continuous limits of large plant-pollinator random networks and
some applications*

Sylvain Billiard†, Hélène Leman‡, Thomas Rey§and Viet Chi Tran¶

June 23, 2022

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

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 inter-actions: 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.

**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

<footnote>*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).</footnote>

<footnote>†Univ. Lille, CNRS, UMR 8198 - Evo-Eco-Paleo, F-59000 Lille, France. sylvain.billiard@univ-lille.fr</footnote>

<footnote>‡Inria, ENSL, UMPA, CNRS UMR 5669, 69364 Lyon, France. helene.leman@inria.fr</footnote>

<footnote>§Univ. Lille, CNRS, Inria, UMR 8524 - Laboratoire P. Painlevé, F-59000 Lille, France. thomas.rey@univ-lille.fr</footnote>

<footnote>¶LAMA, Univ Gustave Eiffel, Univ Paris Est Creteil, CNRS, F-77454 Marne-la-Vallée, France. chi.tran@univ-eiffel.fr</footnote>
# Continuous limits of large plant-pollinator random networks and some applications*

Sylvain Billiard†, Hélène Leman‡, Thomas Rey§and Viet Chi Tran¶

August 26, 2025

## 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 the densities of pollinators. 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 interaction graph and the quantities 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.

**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.] may be studied. In this paper, we address this question by challenging because the number of species in a community is large, and so is the number of possible

<footnote>*This work is dedicated to S. Méléard. 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. This work is funded by the European Union (ERC, SINGER, 101054787). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This work has also been supported by the Chair "Modélisation Mathématique et Biodiversité" of Veolia Environnement-Ecole Polytechnique-Museum National d'Histoire Naturelle-Fondation X. T.R. and V.C.T. acknowledge the support of the CDP C2EMPI, as well as the French State under the France-2030 programme, the University of Lille, the Initiative of Excellence of the University of Lille, the European Metropolis of Lille for their funding and support of the R-CDP-24-004-C2EMPI project.</footnote>

<footnote>†Univ. Lille, CNRS, UMR 8198 - Evo-Eco-Paleo, F-59000 Lille, France. sylvain.billiard@univ-lille.fr</footnote>

<footnote>‡Inria, ENSL, UMPA, CNRS UMR 5669, 69364 Lyon, France. helene.leman@inria.fr</footnote>

<footnote>§Univ. Côte d'Azur, CNRS, LJAD, F-06108 Nice, France. thomas.rey@univ-lille.fr</footnote>

<footnote>¶Univ. Lille, CNRS, Inria, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille, France. viet-chi.tran@inria.fr</footnote>
observations thus needs much effort, without being ensured that the inferred networks are representa-
tive 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. [39, 42, 0, 3, 32, 27] and [2] for a re-
view). 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. [14, 47, 39, 42]. 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* × *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
edges. Describing ecological networks from observations thus requires much effort, without being en-
sured 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 com-
munities 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, 45, 46, 47, 38] and [2] for a review). This makes the identification of general prop-
erties of ecological communities difficult. Overall, despite being more than 50 years old, the question of
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 trait variations
between or within species such as size, symmetry, or time of activity, see e.g. [15, 50, 48, 4]. Starting
from a stochastic individual-based model, our goal is to provide simplifications of an ecological network
by 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 due to resource exchanges. On the one hand, plants benefit from pollinator visits 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 plants, which can negatively affect plant demographic rates, proportionally to pollinator densities [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 generalizes Erdös-Rényi graphs that correspond to the case 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 [47, 29].

In Section 1, we present the stochastic individual-based model, how the random interaction 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 those 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 [43], is not a novelty for modelling 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 [44] for the mathematical modeling of evolutionary ecology. However, the
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 species is characterized by a trait $x \in [0, 1]$ and each pollinator species by a trait $y \in [0, 1]$. For $i \in \{1, \dots, n\} = [1, n]$ and $j \in \{1, \dots, m\} = [1, m]$, 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_{ij}^{n,m} = \begin{cases} 1 & \text{if the pollinator species } j \text{ can interact with the plant species } i, \\ 0 & \text{otherwise.} \end{cases}$$

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 [1, n] \times [1, m]}$, the random variables $(G_{ij}^{n,m})_{(i,j) \in [1, n] \times [1, m]}$ are independent and distributed as Bernoulli random variables with parameters $(\phi(x^i, y^j))_{(i,j) \in [1, n] \times [1, m]}$.
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 flowers 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 positions in the plant-pollinator network, determined by their species traits, 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 [25, 23, 42] 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 species is characterized by a trait $x \in [0, 1]$ and each pollinator species by a trait $y \in [0, 1]$. For $i \in \{1, \dots, n\} = [1, n]$ and $j \in \{1, \dots, m\} = [1, m]$, we respectively denote by $x^i$ and $y^j$ the traits of the plant species $i$ and of the pollinator species $j$. These traits can represent for instance their degrees of generalism, i.e. their tendencies 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_{ij}^{n,m} = \begin{cases} 1 & \text{if the plant species } i \text{ can interact with the pollinator species } j, \\ 0 & \text{otherwise.} \end{cases}$$

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 the plant species $i$ and the 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 $\phi$ is Lipschitz continuous with respect to (w.r.t.) both variables and such that conditionally on the
<figure><img src="image_2.png" /><figcaption>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]</figcaption></figure>

**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. \qquad (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_{ij}^{n,m} = 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[1,n]\times[1,m]}$, $(c_{ij}^{n,m})_{(i,j)\in[1,n]\times[1,m]}$ 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 \to \mathbb{R}$: for all $i \in [1, n]$, $j \in [1, m]$:

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

with $\mathbb{E}_{x,y}[A] = \mathbb{E}[A|(x^i, y^j)_{(i,j)\in[1,n]\times[1,m]}]$,
traits $(x^i, y^j)_{i,j \in [1,n] \times [1,m]}$, the random variables $(G_{ij}^{n,m})_{(i,j) \in [1,n] \times [1,m]}$ are independent and distributed as Bernoulli random variables with parameters $(\phi(x^i, y^j))_{(i,j) \in [1,n] \times [1,m]}$.

<figure><img src="image_3.png" /><figcaption>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]</figcaption></figure>

**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 degree distributions 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. \qquad (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 j$ and $J \ni i$. 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). In this case, the function $\phi$ is not Lipschitz continuous, which is a technical assumption needed in Section 3 but it is possible to smooth it in this case.

If species $i$ and $j$ interact, i.e. if $G_{ij}^{n,m} = 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 [1,n] \times [1,m]}$, $(c_{ij}^{n,m})_{(i,j) \in [1,n] \times [1,m]}$ is assumed to be a sequence of independent non-negative r. v. such that
• 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 [1, n] \times [1, m]} \mathrm{Var}_{x,y}((m+n)c_{i,j}^{n,m}) < \infty. \quad (1.3)$$

**Example 1.4.** Depending on plants and pollinators' traits considered, different function $c_{i,j}^{n,m}$ be considered:

• The intensity and frequency of the interaction can be independent of the species traits, i.e. $c_{i,j}^{n,m} = \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_{i,j}^{n,m}$ a function measuring the overlap of activity for both species; Similarly, $c_{i,j}^{n,m}$ 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_{i,j}^{n,m}$ some morphological preference; For instance, assuming $(x^i, y^j) \in [0, 1]^2$, $c_{i,j}^{n,m} = \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_{i,j}^{n,m} = \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 \to \infty$.

Given the scaling factor K, we denote $P_t^{K,i}$ and $A_t^{K,j}$ 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}_t^{K,n,m}(dx) = \frac{1}{nK} \sum_{i=1}^{n} P_t^{K,i} \delta_{x^i}(dx), \qquad \mathbf{A}_t^{K,n,m}(dy) = \frac{1}{mK} \sum_{j=1}^{m} A_t^{K,j} \delta_{y^j}(dy). \quad (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_t^{K,j}/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_t^{A,K,i} := \sum_{j \sim i} \frac{c_{ij}^{n,m} A_t^{K,j}}{K}. \quad (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_t^{P,K,j} = \sum_{i \sim j} \frac{c_{ij}^{n,m} P_t^{K,i}}{K}. \quad (1.6)$$
* the expected values only depend on the traits of the plants and pollinators through a function $c^{n,m} : [0, 1]^2 \to \mathbb{R}$: for all $i \in [1, n]$, $j \in [1, m]$:

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

with $\mathbb{E}_{x,y}[A] = \mathbb{E}[A|(x^i, y^j)_{(i,j)\in[1,n]\times[1,m]}]$ for any bounded random variable A;

* there exists $K_c > 0$ such that almost surely,

$$ \sup_{n \in \mathbb{N}, m \in \mathbb{N}} \sup_{(i,j) \in [1, n] \times [1, m]} (m+n)c_{i,j}^{n,m} < K_c, $$

in particular, this implies that

$$ V_{max} := \sup_{n \in \mathbb{N}, m \in \mathbb{N}} \sup_{(i,j) \in [1, n] \times [1, m]} \mathrm{Var}_{x,y}((m+n)c_{i,j}^{n,m}) < \infty, \qquad (1.3) $$

where $\mathrm{Var}_{x,y}[A] = \mathrm{Var}[A|(x^i, y^j)_{(i,j)\in[1,n]\times[1,m]}]$.

**Example 1.4.** Depending on plants and pollinators' traits considered, different functions $c^{n,m}$ be considered:

* The intensity and frequency of the interaction can be independent of the species traits, i.e. $c_{i,j}^{n,m} = \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_{i,j}^{n,m}$ a measure of the overlap of activities for both species; Similarly, $c_{i,j}^{n,m}$ can represent geographical proximity 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' bodies and the r.v. $c_{i,j}^{n,m}$ some morphological preferences. For instance, assuming $(x^i, y^j) \in [0, 1]^2$, $c_{i,j}^{n,m} = \frac{x^i y^j}{n+m}$ would mean that pollinators and plants tend to interact more when the traits are both large, while assuming $c_{i,j}^{n,m} = \frac{x^i (1-y^j)}{n+m}$ would mean that they will interact more when the plant has a large trait value and the pollinator a low trait value.

## 1.3 Stochastic dynamics of the plant-pollinator community

The dynamics of 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 \to \infty$.

Given the scaling factor $K$, we denote $P_t^{K,i}$ and $A_t^{K,j}$ the sizes 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}_t^{K,n,m}(dx) = \frac{1}{nK} \sum_{i=1}^{n} P_t^{K,i} \delta_{x^i}(dx), \qquad \mathbf{A}_t^{K,n,m}(dy) = \frac{1}{mK} \sum_{j=1}^{m} A_t^{K,j} \delta_{y^j}(dy). \qquad (1.4) $$

The interactions between plants and pollinators are characterized by exchanges of resources [32]. The quantities of resources gained by plants or pollinators $R$ are modeled through a "mass-action
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), \qquad (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}_t^{K,n,m}(x) := \int_{[0,1]} k(x, x') d\mathbf{P}_t^{K,n,m}(dx') = \frac{1}{nK} \sum_{i=1}^{n} k(x, x^i) P_t^{K,i}, \qquad (1.8)$$

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

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

where $h(y, y')$ quantifies the competition of individuals with trait $y'$ 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}_t^{K,n,m}, \mathbf{A}_t^{K,n,m})_{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_D^A(ds, dk, d\theta)$) be Poisson point measures on $\mathbb{R}_+ \times E := \mathbb{R}_+ \times [1, n] \times \mathbb{R}_+$ (resp. $\mathbb{R}_+ \times F := \mathbb{R}_+ \times [1, m] \times \mathbb{R}_+$) with intensity measure $q(ds, dk, d\theta) = ds \mathbf{n}(dk) d\theta$ where $ds$ and $d\theta$ are Lebesgue measures on $\mathbb{R}_+$ and where $\mathbf{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

$$P_t^{K,i} = P_0^{K,i} + \int_0^t \int_E \mathbb{1}_{[i=k]} \mathbb{1}_{\theta \le b^P(R_{s-}^{A,K,i})} P_{s-}^{K,i} Q_B^P(ds, dk, d\theta) \\ - \int_0^t \int_E \mathbb{1}_{[i=k]} \mathbb{1}_{\theta \le [d^P(R_{s-}^{A,K,i}) + k \star \mathbf{P}_{s-}^{K,n}(x^i)]} P_{s-}^{K,i} Q_D^P(ds, dk, d\theta) \qquad (1.10)$$

$$A_t^{K,j} = A_0^{K,j} + \int_0^t \int_F \mathbb{1}_{[j=k]} \mathbb{1}_{\theta \le b^A(R_{s-}^{P,K,j})} A_{s-}^{K,j} Q_B^A(ds, dk, d\theta) \\ - \int_0^t \int_F \mathbb{1}_{[j=k]} \mathbb{1}_{\theta \le [d^A(R_{s-}^{P,K,j}) + h \star \mathbf{A}_{s-}^{K,m}(y^j)]} A_{s-}^{K,j} Q_D^A(ds, dk, d\theta). \qquad (1.11)$$
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_t^{K,j}/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_t^{A,K,i} := \sum_{j \sim i} \frac{c_{ij}^{n,m} A_t^{K,j}}{K}. \qquad (1.5)$$

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

$$R_t^{P,K,j} = \sum_{i \sim j} \frac{c_{ij}^{n,m} P_t^{K,i}}{K}. \qquad (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 quantities 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), \qquad (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}_t^{K,n,m}(x) := \int_{[0,1]} k(x, x') d\mathbf{P}_t^{K,n,m}(dx') = \frac{1}{nK} \sum_{i=1}^{n} k(x, x^i) P_t^{K,i}, \qquad (1.8)$$

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

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

where $h(y, y')$ quantifies the competition of individuals with trait $y'$ on individuals with trait $y$.

The following assumptions on the functions involved in the model will be needed, both for modelling 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 $\mathbb{R}^+$.

(ii) The competition kernels $k$ and $h$ are assumed to be Lipschitz continuous w.r.t. both variables 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}_t^{K,n,m}, \mathbf{A}_t^{K,n,m})_{t \in \mathbb{R}_+}$ defined in (1.4) by SDEs driven by Poisson point measures. Let us first present these SDEs and then explain their heuristic meaning.
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 $1_{i=k}$ or $1_{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_s^{A,K,i}$ and $R_s^{P,K,j}$ (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}_0^{K,n,m}, 1 \rangle^2 + \langle \mathbf{A}_0^{K,n,m}, 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 \ge 0$,

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

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

$$ \begin{aligned} \langle \mathbf{P}_t^{K,n,m}, f \rangle &= \frac{1}{n} \sum_{i=1}^{n} \frac{1}{K} P_t^{K,i} f(x^i) \\ &= \langle \mathbf{P}_0^{K,n,m}, f \rangle + \int_0^t \frac{1}{nK} \sum_{i=1}^{n} f(x^i) [g^P(R_{s-}^{A,K,i}) - k \star \mathbf{P}_s^{K,n,m}(x^i)] P_s^{K,i} ds \\ &\quad + \frac{1}{nK} \sum_{i=1}^{n} f(x^i) M_t^{K,i}, \end{aligned} \quad (1.13) $$

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

$$ \langle M^{K,i} \rangle_t = \int_0^t (b^P(R_s^{A,K,i}) + d^P(R_s^{A,K,i}) + k \star \mathbf{P}_s^{K,n,m}(y^i)) P_s^{K,i} ds. \quad (1.14) $$

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

$$ \begin{aligned} \langle \mathbf{A}_t^{K,n,m}, f \rangle &= \frac{1}{m} \sum_{j=1}^{m} \frac{1}{K} A_t^{K,j} f(y^j) \\ &= \langle \mathbf{A}_0^{K,n,m}, f \rangle + \int_0^t \frac{1}{mK} \sum_{j=1}^{m} f(y^j) [g^A(R_{s-}^{P,K,j}) - h \star \mathbf{A}_s^{K,n,m}(y^j)] A_s^{K,j} ds \\ &\quad + \frac{1}{mK} \sum_{j=1}^{m} f(y^j) M_t^{K,j}, \end{aligned} \quad (1.15) $$

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

$$ \langle M^{K,j} \rangle_t = \int_0^t (b^A(R_s^{P,K,j}) + d^A(R_s^{P,K,j}) + h \star \mathbf{A}_s^{K,n,m}(y^j)) A_s^{K,j} ds. $$
**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_D^A(ds, dk, d\theta)$) be Poisson point measures on $\mathbb{R}_+ \times E := \mathbb{R}_+ \times [1, n] \times \mathbb{R}_+$ (resp. $\mathbb{R}_+ \times F := \mathbb{R}_+ \times [1, m] \times \mathbb{R}_+$) with intensity measure $q(ds, dk, d\theta) = ds \mathbf{n}(dk) d\theta$ where $ds$ and $d\theta$ are Lebesgue measures on $\mathbb{R}_+$ and where $\mathbf{n}(dk)$ is the counting measure on $\mathbb{N}^* = \{1, 2, \dots\}$. Consider initial conditions $P_0^{K,i}$'s and $A_0^{K,j}$'s that are independent from the Poisson point measures and from the $G_{ij}^{n,m}$'s and $c_{ij}^{n,m}$'s.

For $i \in [1, n]$ and $j \in [1, m]$, one has

$$ P_t^{K,i} = P_0^{K,i} + \int_0^t \int_E \mathbb{1}_{i=k} \mathbb{1}_{\theta \le b^P(R_{s-}^{A,K,i})} P_{s-}^{K,i} Q_B^P(ds, dk, d\theta) \\ - \int_0^t \int_E \mathbb{1}_{i=k} \mathbb{1}_{\theta \le [d^P(R_{s-}^{A,K,i}) + k \star \mathbf{P}_{s-}^{K,n,m}(x^i)]} P_{s-}^{K,i} Q_D^P(ds, dk, d\theta) \quad (1.10) $$

$$ A_t^{K,j} = A_0^{K,j} + \int_0^t \int_F \mathbb{1}_{j=k} \mathbb{1}_{\theta \le b^A(R_{s-}^{P,K,j})} A_{s-}^{K,j} Q_B^A(ds, dk, d\theta) \\ - \int_0^t \int_F \mathbb{1}_{j=k} \mathbb{1}_{\theta \le [d^A(R_{s-}^{P,K,j}) + h \star \mathbf{A}_{s-}^{K,n,m}(y^j)]} A_{s-}^{K,j} Q_D^A(ds, dk, d\theta). \quad (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 $\mathbb{1}_{i=k}$ or $\mathbb{1}_{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_s^{A,K,i}$ and $R_s^{P,K,j}$ (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). For a finite measure $\mu$ and a positive measurable function $f$, we will denote by $\langle \mu, f \rangle = \int f d\mu$ the integral of $f$ w.r.t. $\mu$.

**Proposition 1.7.** Consider here $n$, $m$ and $K$ fixed. If Assumption 1.5 holds and if

$$ \mathbb{E} [\langle \mathbf{P}_0^{K,n,m}, \mathbb{1} \rangle^2 + \langle \mathbf{A}_0^{K,n,m}, \mathbb{1} \rangle^2] < +\infty, $$

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

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

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

$$ \begin{aligned} \langle \mathbf{P}_t^{K,n,m}, f \rangle &= \frac{1}{n} \sum_{i=1}^n \frac{1}{K} P_t^{K,i} f(x^i) \\ &= \langle \mathbf{P}_0^{K,n,m}, f \rangle + \int_0^t \frac{1}{nK} \sum_{i=1}^n f(x^i) [g^P(R_{s-}^{A,K,i}) - k \star \mathbf{P}_s^{K,n,m}(x^i)] P_s^{K,i} ds \\ &\quad + \frac{1}{nK} \sum_{i=1}^n f(x^i) M_t^{K,i}, \end{aligned} \quad (1.13) $$
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}_t^{K,n,m})_{t \in \mathbb{R}_+}$ and $(\mathbf{A}_t^{K,n,m})_{t \in \mathbb{R}_+}$ are well-defined in the set $\mathcal{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}_0^{K,n,m}, \mathbf{1} \rangle^3 + \langle \mathbf{A}_0^{K,n,m}, \mathbf{1} \rangle^3 \right] < +\infty, \quad (2.1) $$

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

$$ \lim_{K \to +\infty} \frac{P_0^{K,i}}{K} = \tilde{P}_0^i, \quad \lim_{K \to +\infty} \frac{A_0^{K,j}}{K} = \tilde{A}_0^j. $$

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

$$ \lim_{K \to \infty} \sup_{t \le T} \sup_{i,j} \left\{ \left| \frac{P_t^{K,i}}{K} - \tilde{P}_t^i \right|, \left| \frac{A_t^{K,j}}{K} - \tilde{A}_t^j \right| \right\} = 0, \quad (2.2) $$

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

$$ \forall 1 \le i \le n, \quad \frac{d\tilde{P}_t^i}{dt} = \left( g^P \left( \sum_{j \sim i} \tilde{R}_t^{A,ij} \right) - \frac{1}{n} \sum_{\ell=1}^n k(x^i, x^\ell) \tilde{P}_t^\ell \right) \tilde{P}_t^i \quad (2.3) $$

$$ \forall 1 \le j \le m, \quad \frac{d\tilde{A}_t^j}{dt} = \left( g^A \left( \sum_{i \sim j} \tilde{R}_t^{P,ij} \right) - \frac{1}{m} \sum_{\ell=1}^m h(y^j, y^\ell) \tilde{A}_t^\ell \right) \tilde{A}_t^j, $$

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

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.
where $(M^{K,i})_{i \in \{1, \dots, n\}}$ are square integrable martingales with predictable quadratic variation processes:

$$ \langle M^{K,i} \rangle_t = \int_0^t (b^P(R_s^{A,K,i}) + d^P(R_s^{A,K,i}) + k \star \mathbf{P}_s^{K,n,m}(x^i)) P_s^{K,i} ds. \quad (1.14) $$

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

$$ \begin{aligned} \langle \mathbf{A}_t^{K,n,m}, f \rangle &= \frac{1}{m} \sum_{j=1}^{m} \frac{1}{K} A_t^{K,j} f(y^j) \\ &= \langle \mathbf{A}_0^{K,n,m}, f \rangle + \int_0^t \frac{1}{mK} \sum_{j=1}^{m} f(y^j) [g^A(R_{s-}^{P,K,j}) - h \star \mathbf{A}_s^{K,n,m}(y^j)] A_s^{K,j} ds \\ &\quad + \frac{1}{mK} \sum_{j=1}^{m} f(y^j) M_t^{K,j}, \end{aligned} \quad (1.15) $$

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

$$ \langle M^{K,j} \rangle_t = \int_0^t (b^A(R_s^{P,K,j}) + d^A(R_s^{P,K,j}) + h \star \mathbf{A}_s^{K,n,m}(y^j)) A_s^{K,j} 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}_t^{K,n,m})_{t \in \mathbb{R}_+}$ and $(\mathbf{A}_t^{K,n,m})_{t \in \mathbb{R}_+}$ are well-defined in the set $\mathcal{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 parameter $K$ rescaling the population sizes of plants and pollinators tends to $+\infty$. In this section, since $n$ and $m$ are fixed, these letters will be omitted 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}_0^{K,n,m}, \mathbf{1} \rangle^3 + \langle \mathbf{A}_0^{K,n,m}, \mathbf{1} \rangle^3 \right] < +\infty, \quad (2.1) $$

and such that there exist $(\tilde{P}_0^1, \dots, \tilde{P}_0^n)$ and $(\tilde{A}_0^1, \dots, \tilde{A}_0^m)$, independent from the $G_{ij}^{n,m}$'s and $c_{ij}^{n,m}$'s, satisfying the following convergences almost surely (a.s.) :

$$ \lim_{K \to +\infty} \frac{P_0^{K,i}}{K} = \tilde{P}_0^i, \quad \lim_{K \to +\infty} \frac{A_0^{K,j}}{K} = \tilde{A}_0^j. $$

Then, for all $T \ge 0$, the following convergence holds a.s.:

$$ \lim_{K \to \infty} \sup_{t \le T} \sup_{i,j} \left\{ \left| \frac{P_t^{K,i}}{K} - \tilde{P}_t^i \right|, \left| \frac{A_t^{K,j}}{K} - \tilde{A}_t^j \right| \right\} = 0, \quad (2.2) $$
Proof. Since the numbers $n$ and $m$ of species are supposed constant, working with the vector processes $(P^{K,j}, A^{K,j})_{i \in [1,n], j \in [1,m]}$ 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} \left[ \sup_{t \in [0,T]} \langle \mathbf{P}_t^{K,n,m}, \mathbf{1} \rangle^3 + \langle \mathbf{A}_t^{K,n,m}, \mathbf{1} \rangle^3 \right] < +\infty. \quad (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 $\mathcal{D}([0, T], \mathcal{M}_F([0, 1])^2)$ to the process $(\tilde{\mathbf{P}}^n, \tilde{\mathbf{A}}^m)$ defined by:

$$ \tilde{\mathbf{P}}_t^n(dx) = \frac{1}{n} \sum_{i=1}^{n} \tilde{P}_t^i \delta_{x^i}(dx), \quad \tilde{\mathbf{A}}_t^m(dy) = \frac{1}{m} \sum_{j=1}^{m} \tilde{A}_t^j \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_1^{P,n,m}, \dots, \Phi_n^{P,n,m})$ and $\Phi^{A,n,m} = (\Phi_1^{A,n,m}, \dots, \Phi_m^{A,n,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\}$:

$$ \frac{d\tilde{P}_t^i}{dt} = \Phi_i^{P,n,m}(\tilde{P}_t, \tilde{A}_t) \quad \text{and} \quad \frac{d\tilde{A}_t^j}{dt} = \Phi_j^{A,n,m}(\tilde{P}_t, \tilde{A}_t). \quad (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 \to +\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_t^{K,P} = \sqrt{K} \left( \begin{array}{c} P_t^{K,1}/K - \tilde{P}_t^1 \\ \vdots \\ P_t^{K,n}/K - \tilde{P}_t^n \end{array} \right), \quad \text{and} \quad \eta_t^{K,A} = \sqrt{K} \left( \begin{array}{c} A_t^{K,1}/K - \tilde{A}_t^1 \\ \vdots \\ A_t^{K,m}/K - \tilde{A}_t^m \end{array} \right). \quad (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 $C^1$, and that the initial conditions $(\eta_0^{K,P}, \eta_0^{K,A})_{K \ge 0}$ converge in distribution towards a
where $(\tilde{P}_t^1, \dots, \tilde{P}_t^n)_{t \ge 0}$ and $(\tilde{A}_t^1, \dots, \tilde{A}_t^m)_{t \ge 0}$ are the unique solution of the system, which exists on $\mathbb{R}^{n+m}$:

$$ 
\begin{aligned}
\forall i \in [1, n], \quad \frac{d\tilde{P}_t^i}{dt} &= \left( g^P(\tilde{R}_t^{A,i}) - \frac{1}{n} \sum_{\ell=1}^n k(x^i, x^\ell) \tilde{P}_t^\ell \right) \tilde{P}_t^i \\
j \in [1, m], \quad \frac{d\tilde{A}_t^j}{dt} &= \left( g^A(\tilde{R}_t^{P,j}) - \frac{1}{m} \sum_{\ell=1}^m h(y^j, y^\ell) \tilde{A}_t^\ell \right) \tilde{A}_t^j,
\end{aligned}
\quad (2.3) $$

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

Equations (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 traits $x^i$ of plants and $y^j$ of pollinators, and these traits 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). These functions can take any form (see Section 4 for some examples). As a consequence, (2.3) can capture a large variety of ecological situations. In particular, many ODE models published in the ecological literature are special cases of (2.3) (e.g. [47, 8, 38]) and therefore can be seen as limits of stochastic individual-based models as defined in Section 1. Notice that coefficients $c_{ij}^{n,m}$ may be random coefficients as explained in the previous section. However, these coefficients remain constant over 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 [1,n], j \in [1,m]}$ 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} \left[ \sup_{t \in [0,T]} \langle \mathbf{P}_t^{K,n,m}, \mathbf{1} \rangle^3 + \langle \mathbf{A}_t^{K,n,m}, \mathbf{1} \rangle^3 \right] < +\infty. \quad (2.4) $$

The global existence and uniqueness of solution to (2.3) can be proved by classical results for ODEs: 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 $\mathcal{D}([0, T], \mathcal{M}_F([0, 1])^2)$ to the process $(\tilde{\mathbf{P}}^n, \tilde{\mathbf{A}}^m)$ defined by:

$$ \tilde{\mathbf{P}}_t^n(dx) = \frac{1}{n} \sum_{i=1}^n \tilde{P}_t^i \delta_{x^i}(dx), \quad \tilde{\mathbf{A}}_t^m(dy) = \frac{1}{m} \sum_{j=1}^m \tilde{A}_t^j \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_1^{P,n,m}, \dots, \Phi_n^{P,n,m})$ and $\Phi^{A,n,m} = (\Phi_1^{A,n,m}, \dots, \Phi_m^{A,n,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, n]$ and $j \in [1, m]$:

$$ \frac{d\tilde{P}_t^i}{dt} = \Phi_i^{P,n,m}(\tilde{P}_t, \tilde{A}_t) \quad \text{and} \quad \frac{d\tilde{A}_t^j}{dt} = \Phi_j^{A,n,m}(\tilde{P}_t, \tilde{A}_t). \quad (2.5) $$
deterministic vector $(\tilde{\eta}_0^P, \tilde{\eta}_0^A)$ when $K \to +\infty$. Under the same assumptions as Proposition 2.1, one has

$$ (\eta_t^{K,P}, \eta_t^{K,A})_{t \ge 0} \underset{K \to \infty}{\Rightarrow} (\tilde{\eta}_t^P, \tilde{\eta}_t^A)_{t \ge 0}, $$

where the converge holds in law in the Skorohod space $D(\mathbb{R}_+, \mathbb{R}^{n+m})$ and where the processes $(\tilde{\eta}_i^P, \tilde{\eta}_i^A) = (\tilde{\eta}_i^{P,i}, \tilde{\eta}_i^{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 \le i \le n, 1 \le j \le m}$ in $\mathbb{R}$:
for all $1 \le i \le n$ and $1 \le j \le m$,

$$ \begin{aligned} \tilde{\eta}_t^{P,i} = & \tilde{\eta}_0^{P,i} + \int_0^t \sqrt{ \left( b^P(\tilde{R}_s^{A,i}) + d^P(\tilde{R}_s^{A,i}) + \frac{1}{n} \sum_{\ell=1}^n k(x^i, x^\ell) \tilde{P}_s^\ell \right) \tilde{P}_s^i } dW_s^{P,i} \\ & + \int_0^t \left[ \sum_{\ell=1}^n \frac{\partial \Phi_i^{P,n,m}}{\partial p^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{P,\ell} + \sum_{\ell=1}^m \frac{\partial \Phi_i^{P,n,m}}{\partial a^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{A,\ell} \right] ds, \end{aligned} \tag{2.7} $$ $$ \begin{aligned} \tilde{\eta}_t^{A,j} = & \tilde{\eta}_0^{A,j} + \int_0^t \sqrt{ \left( b^A(\tilde{R}_s^{P,j}) + d^A(\tilde{R}_s^{P,j}) + \frac{1}{m} \sum_{\ell=1}^m h(y^j, y^\ell) \tilde{A}_s^\ell \right) \tilde{A}_s^j } dW_s^{A,j} \\ & + \int_0^t \left[ \sum_{\ell=1}^n \frac{\partial \Phi_j^{A,n,m}}{\partial p^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{P,\ell} + \sum_{\ell=1}^m \frac{\partial \Phi_j^{A,n,m}}{\partial a^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{A,\ell} \right] ds. \end{aligned} $$

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 \to +\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_t^{K,i}$ and $A_t^{K,j}$ (1.13)-(1.15) and (2.5). To understand (2.7), consider for example the ith component of $\eta_t^{K,P,i}$, $\eta_t^{K,P,i}$. The other terms can be treated similarly. We sketch here this alternative proof:

$$ \begin{aligned} \eta_t^{K,P,i} = & \eta_0^{K,P,i} + \sqrt{K} \int_0^t \frac{1}{n} \sum_{\ell=1}^n k(x^\ell, x^i) \left( \tilde{P}_s^\ell \tilde{P}_s^i - \frac{P_s^{K,\ell}}{K} \frac{P_s^{K,i}}{K} \right) ds + \sqrt{K} \frac{M_t^{K,i}}{K} \\ & + \sqrt{K} \int_0^t \left[ g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \frac{A_s^{K,j}}{K} \right) \frac{P_s^{K,i}}{K} - g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_s^j \right) \tilde{P}_s^i \right] ds \\ = & \eta_0^{K,P,i} + \frac{1}{n} \sum_{\ell=1}^n k(x^\ell, x^i) \int_0^t \left[ \eta_s^{K,P,\ell} \frac{P_s^{K,i}}{K} + \tilde{P}_s^\ell \eta_s^{K,P,i} \right] ds + \sqrt{K} \frac{M_t^{K,i}}{K} \tag{2.8} \\ & + \int_0^t \left[ g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \frac{A_s^{K,j}}{K} \right) \eta_s^{K,P,i} + \sqrt{K} \tilde{P}_s^i \left( g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \frac{A_s^{K,j}}{K} \right) - g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_s^j \right) \right) \right] ds, \end{aligned} $$

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 $C^1$, so is $g^P$ and a Taylor expansion can be used for the last
## 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, while the ODEs (2.3) are often used in the modelling papers. The law of large number (Proposition 2.1) shows that both points of view are related in the limit $K \to +\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_t^{K,P} = \sqrt{K} \begin{pmatrix} P_t^{K,1}/K - \tilde{P}_t^1 \\ \vdots \\ P_t^{K,n}/K - \tilde{P}_t^n \end{pmatrix}, \quad \text{and} \quad \eta_t^{K,A} = \sqrt{K} \begin{pmatrix} A_t^{K,1}/K - \tilde{A}_t^1 \\ \vdots \\ A_t^{K,m}/K - \tilde{A}_t^m \end{pmatrix}. \quad (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 $C^1$, and that the initial conditions $(\eta_0^{K,P}, \eta_0^{K,A})_{K \ge 0}$ converge in distribution towards a deterministic vector $(\tilde{\eta}_0^P, \tilde{\eta}_0^A)$ when $K \to +\infty$. Under the same assumptions as Proposition 2.1, one has

$$ (\eta_t^{K,P}, \eta_t^{K,A})_{t \ge 0} \underset{K \to \infty}{\Rightarrow} (\tilde{\eta}_t^P, \tilde{\eta}_t^A)_{t \ge 0}, $$

where the converge holds in distribution in the Skorohod space $D(\mathbb{R}_+, \mathbb{R}^{n+m})$ and where the processes $(\tilde{\eta}_t^P, \tilde{\eta}_t^A) = (\eta_t^{P,i}, \eta_t^{A,j})_{i \in [1,n], j \in [1,m]}$ are solutions of the following SDEs driven by $n+m$ independent standard real-valued Brownian motions $(W^{P,i}, W^{A,j})_{i \in [1,n], j \in [1,m]}$: for all $i \in [1,n]$ and $j \in [1,m]$,

$$ \begin{aligned} \tilde{\eta}_t^{P,i} = & \tilde{\eta}_0^{P,i} + \int_0^t \sqrt{ \left( b^P(\tilde{R}_s^{A,i}) + d^P(\tilde{R}_s^{A,i}) + \frac{1}{n} \sum_{\ell=1}^n k(x^i, x^\ell) \tilde{P}_s^\ell \right) \tilde{P}_s^i } dW_s^{P,i} \\ & + \int_0^t \left[ \sum_{\ell=1}^n \frac{\partial \Phi_i^{P,n,m}}{\partial p^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{P,\ell} + \sum_{\ell=1}^m \frac{\partial \Phi_i^{P,n,m}}{\partial a^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{A,\ell} \right] ds, \end{aligned} \quad (2.7) \\ \tilde{\eta}_t^{A,j} = & \tilde{\eta}_0^{A,j} + \int_0^t \sqrt{ \left( b^A(\tilde{R}_s^{P,j}) + d^A(\tilde{R}_s^{P,j}) + \frac{1}{m} \sum_{\ell=1}^m h(y^j, y^\ell) \tilde{A}_s^\ell \right) \tilde{A}_s^j } dW_s^{A,j} \\ & + \int_0^t \left[ \sum_{\ell=1}^n \frac{\partial \Phi_j^{A,n,m}}{\partial p^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{P,\ell} + \sum_{\ell=1}^m \frac{\partial \Phi_j^{A,n,m}}{\partial a^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{A,\ell} \right] ds. \end{aligned} $$

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 \to +\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 other papers, 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_t^{K,i}$ and $A_t^{K,j}$ (1.13)-(1.15) and (2.5). To
term:

$$
\begin{align}
\sqrt{K} \left( g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \frac{A_s^{K,j}}{K} \right) - g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_s^j \right) \right) &= (g^P)' \left( \sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_s^j \right) \times \sum_{j \sim i} c_{ij}^{n,m} \sqrt{K} \left( \frac{A_s^{K,j}}{K} - \tilde{A}_s^j \right) + \varepsilon_K \nonumber \\
&= (g^P)' \left( \sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_s^j \right) \times \sum_{j \sim i} c_{ij}^{n,m} \eta_s^{K,A,j} + \varepsilon_K, \tag{2.9}
\end{align}
$$

where $\varepsilon_K$ is a remainder term. From (2.8) and (2.9), we recognize that:

$$
\eta_t^{K,P,i} = \eta_0^{K,P,i} + \sqrt{K} \frac{M_t^{K,i}}{K} + \int_0^t \left[ \sum_{\ell=1}^n \frac{\partial \Phi_i^{P,n,m}}{\partial p^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{P,\ell} + \sum_{\ell=1}^m \frac{\partial \Phi_i^{P,n,m}}{\partial a^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{A,\ell} \right] 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 \to +\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 cumula-
tive 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 \ge 1}$ and $(\tilde{v}_j)_{j \ge 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 \ge 1}$ and
$(v^{j,m})_{j=1...m}$ be the ordered $m^{th}$ first values of $(\tilde{v}_j)_{j \ge 1}$, then for any $i \in \{1, .., n\}$ and $j \in \{1, .., m\}$,

$$
x^{i,n} = F_P^{-1}(u^{i,n}) \quad \text{and} \quad 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_{ij}^{n,m})_{(i,j) \in [1,n] \times [1,m]}$ and the weight of $c_{ij}^{n,m}$ of the interaction $i \sim j$ whenever $G_{ij}^{n,m} = 1$. The random variables $G_{ij}^{n,m}$ 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 [1, n] \times [1, m]$. For the harvesting coefficients $c_{ij}^{n,m}$, 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 \to \mathbb{R}$ to which the sequences of functions $nc^{n,m}(\cdot, \cdot)$ and $mc^{n,m}(\cdot, \cdot)$ 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(\cdot, \cdot)$ is a “harvesting function” underlying the
matrix $(c_{ij}^{n,m})_{i \in [1,n], j \in [1,m]}$ 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 \ge 1}$ such that

$$
m = \alpha_n n \quad \text{and} \quad \lim_{n \to \infty} \alpha_n = 1. \tag{3.1}
$$

In what follows, we write $n \to \infty$, $m \to \infty$ or $n, m \to \infty$ indistinctly.
understand (2.7), consider for example the ith component $\eta_t^{K,P,i}$ of $\eta_t^{K,P}$. The other terms can be treated similarly. We sketch here this alternative proof:

$$
\begin{align}
\eta_t^{K,P,i} &= \eta_0^{K,P,i} + \sqrt{K} \int_0^t \frac{1}{n} \sum_{\ell=1}^n k(x^\ell, x^i) \left( \tilde{P}_s^\ell \tilde{P}_s^i - \frac{P_s^{K,\ell}}{K} \frac{P_s^{K,i}}{K} \right) ds + \sqrt{K} \frac{M_t^{K,i}}{K} \nonumber \\
&\quad + \sqrt{K} \int_0^t \left[ g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \frac{A_s^{K,j}}{K} \right) \frac{P_s^{K,i}}{K} - g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_s^j \right) \tilde{P}_s^i \right] ds \nonumber \\
&= \eta_0^{K,P,i} + \frac{1}{n} \sum_{\ell=1}^n k(x^\ell, x^i) \int_0^t \left[ \eta_s^{K,P,\ell} \frac{P_s^{K,i}}{K} + \tilde{P}_s^\ell \eta_s^{K,P,i} \right] ds + \sqrt{K} \frac{M_t^{K,i}}{K} \tag{2.8} \\
&\quad + \int_0^t \left[ g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \frac{A_s^{K,j}}{K} \right) \eta_s^{K,P,i} + \sqrt{K} \tilde{P}_s^i \left( g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \frac{A_s^{K,j}}{K} \right) - g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_s^j \right) \right) \right] ds, \nonumber
\end{align}
$$

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

$$
\begin{align}
\sqrt{K} \left( g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \frac{A_s^{K,j}}{K} \right) - g^P \left( \sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_s^j \right) \right) &= (g^P)' \left( \sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_s^j \right) \times \sum_{j \sim i} c_{ij}^{n,m} \sqrt{K} \left( \frac{A_s^{K,j}}{K} - \tilde{A}_s^j \right) + \varepsilon_K \nonumber \\
&= (g^P)' \left( \sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_s^j \right) \times \sum_{j \sim i} c_{ij}^{n,m} \eta_s^{K,A,j} + \varepsilon_K, \tag{2.9}
\end{align}
$$

where $\varepsilon_K$ is a remainder term. From (2.8) and (2.9), we recognize that:

$$
\eta_t^{K,P,i} = \eta_0^{K,P,i} + \sqrt{K} \frac{M_t^{K,i}}{K} + \int_0^t \left[ \sum_{\ell=1}^n \frac{\partial \Phi_i^{P,n,m}}{\partial p^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{P,\ell} + \sum_{\ell=1}^m \frac{\partial \Phi_i^{P,n,m}}{\partial a^\ell} (\tilde{P}_s, \tilde{A}_s) \tilde{\eta}_s^{A,\ell} \right] 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 \to +\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 cumu-
lative distribution functions $F_P$ and $F_A$ respectively. To order species according to their respective
traits, we proceed as follows: let $(\tilde{u}_i)_{i \ge 1}$ and $(\tilde{v}_j)_{j \ge 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 \in [1,n]}$ be the ordered $n^{th}$ first values of $(\tilde{u}_i)_{i \ge 1}$
and $(v^{j,m})_{j \in [1,m]}$ be the ordered $m^{th}$ first values of $(\tilde{v}_j)_{j \ge 1}$, then for any $i \in [1, n]$ and $j \in [1, m]$,

$$
x^{i,n} = F_P^{-1}(u^{i,n}) \quad \text{and} \quad y^{j,m} = F_A^{-1}(v^{j,m}).
$$

The indices *n* and *m* will be dropped when no confusion is possible.
**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 \to +\infty} \frac{1}{n} \sum_{i=1}^{n} \tilde{P}_0^i \delta_{x^i} \stackrel{w}{=} \bar{p}_0(x) dx, \quad \lim_{m \to +\infty} \frac{1}{m} \sum_{j=1}^{m} \tilde{A}_0^j \delta_{y^j} \stackrel{w}{=} \bar{a}_0(y) dy \quad \text{a.s.} \quad (3.2) $$

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

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

converge in law in $C([0, T], \mathcal{M}_F([0, 1])^2)$ to a deterministic process $(\bar{\mathbf{P}}, \bar{\mathbf{A}})$ in $C([0, T], \mathcal{M}_F([0, 1])^2)$, such that:

(i) for all $t \ge 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 C([0, 1], \mathbb{R})$, of

$$ \int_0^1 f(x) d\bar{\mathbf{P}}_t(x) = \int_0^1 f(x) \bar{p}_0(x) dx \\ + \int_0^t \int_0^1 f(x) \left[ g^P \left( \int_0^1 c(x, y) \phi(x, y) \bar{a}_s(y) dy \right) - k \star \bar{p}_s(x) \right] \bar{p}_s(x) dx \, ds, \quad (3.3) $$

$$ \int_0^1 f(y) d\bar{\mathbf{A}}_t(y) = \int_0^1 f(y) \bar{a}_0(y) dy \\ + \int_0^t \int_0^1 f(y) \left[ g^A \left( \int_0^1 c(x, y) \phi(x, y) \bar{p}_s(x) dx \right) - h \star \bar{a}_s(y) \right] \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 \to [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$:

$$ \partial_t \bar{p}_t(x) = \left[ g^P \left( \int_0^1 c(x, y) \phi(x, y) \bar{a}_t(y) dy \right) - k \star \bar{p}_t(x) \right] \bar{p}_t(x), \quad (3.4) $$

$$ \partial_t \bar{a}_t(y) = \left[ g^A \left( \int_0^1 c(x, y) \phi(x, y) \bar{p}_t(x) dx \right) - h \star \bar{a}_t(y) \right] \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 \to +\infty$, and the kinetic PDEs (3.3) when additionally $n, m \to +\infty$.
Recall that the plant-pollinator network is defined by its adjacency matrix $(G_{ij}^{n,m})_{(i,j)\in[1,n]\times[1,m]}$ and the weights $c_{ij}^{n,m}$ of the interaction $i \sim j$ whenever $G_{ij}^{n,m} = 1$. The random variables $G_{ij}^{n,m}$ 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 r.v. with parameters $\phi(x^i, y^j)$ for the term $(i, j) \in [1, n] \times [1, m]$. For the harvesting coefficients $c_{ij}^{n,m}$, 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 Lipschitz continuous function $c: [0, 1]^2 \to \mathbb{R}$ to which the sequences of functions $nc^{n,m}(\cdot, \cdot)$ and $mc^{n,m}(\cdot, \cdot)$ 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(\cdot, \cdot)$ is a “harvesting function” underlying the matrix $(c_{ij}^{n,m})_{i \in [1,n], j \in [1,m]}$ of harvesting coefficients. These assumptions are satisfied by Examples 1.4.

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 \ge 1}$ such that

$$m = \alpha_n n \quad \text{and} \quad \lim_{n \to \infty} \alpha_n = 1. \qquad (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 \to +\infty} \frac{1}{n} \sum_{i=1}^{n} \tilde{P}_{0}^{i} \delta_{x^{i}} \stackrel{w}{=} \bar{p}_{0}(x) dx, \quad \lim_{m \to +\infty} \frac{1}{m} \sum_{j=1}^{m} \tilde{A}_{0}^{j} \delta_{y^{j}} \stackrel{w}{=} \bar{a}_{0}(y) dy \quad \text{a.s.} \qquad (3.2)$$

such that for all $n, m \in \mathbb{N}$, the sequences $(\tilde{P}_0^i)_{i \in [1,n]}$ and $(\tilde{A}_0^j)_{j \in [1,m]}$ only depend on $(x^i, y^j)_{(i,j) \in [1,n] \times [1,m]}$. We also assume that there is a constant $C > 0$ such that a.s.,

$$\sup_{n,m} \langle \tilde{\mathbf{P}}_0^{n,m}, 1 \rangle + \langle \tilde{\mathbf{A}}_0^{n,m}, 1 \rangle < C. \qquad (3.3)$$

Then, for any $T \ge 0$, and for $n, m \to +\infty$, the sequence of measure-valued processes

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

converge in law in $C([0, T], \mathcal{M}_F([0, 1])^2)$ to a deterministic process $(\bar{\mathbf{P}}, \bar{\mathbf{A}})$ in $C([0, T], \mathcal{M}_F([0, 1])^2)$, such that:

(i) for all $t \ge 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 C([0, 1], \mathbb{R})$, of

$$\begin{align}
\int_0^1 f(x) d\bar{\mathbf{P}}_t(x) &= \int_0^1 f(x) \bar{p}_0(x) dx \nonumber \\
&\quad + \int_0^t \int_0^1 f(x) \left[ g^P \left( \int_0^1 c(x, y) \phi(x, y) \bar{a}_s(y) dy \right) - k \star \bar{p}_s(x) \right] \bar{p}_s(x) dx \, ds, \tag{3.4} \\
\int_0^1 f(y) d\bar{\mathbf{A}}_t(y) &= \int_0^1 f(y) \bar{a}_0(y) dy \nonumber \\
&\quad + \int_0^t \int_0^1 f(y) \left[ g^A \left( \int_0^1 c(x, y) \phi(x, y) \bar{p}_s(x) dx \right) - h \star \bar{a}_s(y) \right] \bar{a}_s(y) dy \, ds, \nonumber
\end{align}$$

where $k \star \nu$ denotes the convolution product 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.
## 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 $(\tilde{P}_{\infty}^1, \dots, \tilde{P}_{\infty}^n, \tilde{A}_{\infty}^1, \dots, \tilde{A}_{\infty}^m)$ of (2.3) are solutions to the following system:

$$ 
\begin{array}{l}
\forall 1 \le i \le n, \quad \tilde{P}_{\infty}^{i} = 0 \quad \text{or} \quad g^{P}\left(\sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_{\infty}^{j}\right) = \frac{1}{n} \sum_{\ell=1}^{n} k(x^{i}, x^{\ell}) \tilde{P}_{\infty}^{\ell}, \\
\forall 1 \le j \le m, \quad \tilde{A}_{\infty}^{j} = 0 \quad \text{or} \quad g^{A}\left(\sum_{i \sim j} c_{ij}^{n,m} \tilde{P}_{\infty}^{i}\right) = \frac{1}{m} \sum_{\ell=1}^{m} h(y^{j}, y^{\ell}) \tilde{A}_{\infty}^{\ell}.
\end{array}
\qquad (4.1) $$

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

$$ \tilde{P}_{\infty}^{i} = 0 \text{ for all } i \in \{1, ..n\} \quad \text{and} \quad \tilde{A}_{\infty}^{j} = 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. $\tilde{P}_{\infty}^i > 0$ and $\tilde{A}_{\infty}^j > 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) \lor g^A(0) < 0$, the null equilibrium is locally stable. If $g^P(0) \lor 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 & \dots & 0 & 0 \\ 0 & g^P(0) & 0 & \dots & 0 & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & 0 & g^A(0) & 0 \\ 0 & 0 & \dots & 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:

$$ 
\begin{align}
g^P(R^A) &= 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), \tag{4.2} \\
g^A(R^P) &= b^A(R^P) - d^A(R^P) = \frac{\alpha_A R^P}{\beta_A + \gamma_A R^P} - d_A,
\end{align}
$$

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.
The proof of this Proposition is given in Appendix A.2. Notice that assuming (3.3) is not very restrictive in our biological context, since by (3.2), $\langle \tilde{\mathbf{P}}_0^{n,m}, 1 \rangle$ (resp. $\langle \tilde{\mathbf{A}}_0^{n,m}, 1 \rangle$) converges a.s. to $\|\bar{p}_0\|_1$ (resp. $\|\bar{a}_0\|_1$) which is finite. This just ensures that the constant C in (3.3) is not a r.v.

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

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

$$ 
\begin{aligned}
\partial_t \bar{p}_t(x) &= [g^P(\int_0^1 c(x, y)\phi(x, y)\bar{a}_t(y)dy) - k \star \bar{p}_t(x)] \bar{p}_t(x), \\
\partial_t \bar{a}_t(y) &= [g^A(\int_0^1 c(x, y)\phi(x, y)\bar{p}_t(x)dx) - h \star \bar{a}_t(y)] \bar{a}_t(y),
\end{aligned} 
\quad (3.5) $$

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 \to +\infty$, and the kinetic PDEs (3.4) when additionally $n, m \to +\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 $(\tilde{P}_\infty^1, \dots, \tilde{P}_\infty^n, \tilde{A}_\infty^1, \dots, \tilde{A}_\infty^m)$ of (2.3) are solutions to the following system:

$$ 
\begin{aligned}
\forall i \in [1, n], \quad \tilde{P}_\infty^i &= 0 \quad \text{or} \quad g^P\left(\sum_{j \sim i} c_{ij}^{n,m} \tilde{A}_\infty^j\right) = \frac{1}{n} \sum_{\ell=1}^n k(x^i, x^\ell) \tilde{P}_\infty^\ell, \\
\forall j \in [1, m], \quad \tilde{A}_\infty^j &= 0 \quad \text{or} \quad g^A\left(\sum_{i \sim j} c_{ij}^{n,m} \tilde{P}_\infty^i\right) = \frac{1}{m} \sum_{\ell=1}^m h(y^j, y^\ell) \tilde{A}_\infty^\ell.
\end{aligned} 
\quad (4.1) $$

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

$$ \tilde{P}_\infty^i = 0 \quad \text{for all } i \in \{1, \dots, n\} \quad \text{and} \quad \tilde{A}_\infty^j = 0 \quad \text{for all } j \in \{1, \dots, 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. $\tilde{P}_\infty^i > 0$ and $\tilde{A}_\infty^j > 0$, as soon as the interactions between plants and pollinators translate into positive growth rate for both plants and pollinators. Despite that this condition is little
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. The form given here to the growth rates (or

<figure><img src="image_3.png" /><figcaption>Figure 2: Graphical representation of g_A and g_P</figcaption></figure>

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 $\tilde{P}_t^1$ 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:

$$\begin{align}
\frac{dP_t}{dt} &= (g^P(cA_t) - kP_t) P_t \tag{4.3} \\
\frac{dA_t}{dt} &= (g^A(cP_t) - hA_t) A_t \nonumber
\end{align}$$

(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).

Precisely, we prove the following lemma that gives the number of equilibria and their stability.
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 obtain 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 & \dots & 0 & 0 \\ 0 & g^P(0) & 0 & \dots & 0 & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & 0 & g^A(0) & 0 \\ 0 & 0 & \dots & 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 out further computation.

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

$$ \begin{align} g^P(R^A) &= 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), \tag{4.2} \\ g^A(R^P) &= b^A(R^P) - d^A(R^P) = \frac{\alpha_A R^P}{\beta_A + \gamma_A R^P} - d_A, \notag \end{align} $$

where $R^A$ and $R^P$ are the total resources collected respectively by plants and pollinators (Eqs. (1.5) and (1.6)). The 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 for plants $R^A$. 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. 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 $\tilde{P}_t^1$ becomes $P_t$ and $c_{11}^{1,1}$ becomes $c$ for
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

**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 equilibria 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{array}{l} P_{\infty} = \frac{1}{k} g^P(cA_{\infty}) \\ A_{\infty} = \frac{1}{h} g^A \left( \frac{c}{k} g^P(cA_{\infty}) \right), \text{ and } cA_{\infty} \in (C_0^-, C_0^+) \end{array} \right. \qquad (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

$$ \begin{array}{l} g^P(cA_\infty) = kP_\infty \quad \text{or} \quad P_\infty = 0 \\ g^A(cP_\infty) = hA_\infty \quad \text{or} \quad A_\infty = 0. \end{array} \qquad (4.5) $$

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(\frac{c}{k}g^P(cx)) - x$ on $(\frac{C_0^-}{c}, \frac{C_0^+}{c})$. A rapid differentiation gives

$$ f''(x) = \frac{c^3}{hk}(g^P)''(cx)(g^A)'\left(\frac{c}{k}g^P(cx)\right) + \frac{c^3}{(hk)^2}(g^P)'(cx)^2(g^A)''\left(\frac{c}{k}g^P(cx)\right) < 0, $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Graphical representation of gA and gP</figcaption></figure>

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:

$$ 
\begin{aligned}
\frac{dP_t}{dt} &= (g^P(cA_t) - kP_t) P_t \\
\frac{dA_t}{dt} &= (g^A(cP_t) - hA_t) A_t
\end{aligned}
\qquad (4.3) $$

(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. (see Figure 3).

<figure><img src="image_5.png" /><figcaption>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. For the plots, the parameters are $\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$. 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.</figcaption></figure>

Precisely, we prove the following lemma that gives the number of equilibria and their stability.
since $g^A$ is increasing and $g^P$ and $g^A$ are concave functions. Thus, $f'$ is decreasing and

$$f'(x) = \frac{c^2}{hk}(g^P)'(cx)(g^A)'\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 equilibria,

* 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 equilibria, 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_{\infty}^{\pm} & c P_{\infty}^{\pm} (g^P)' (c A_{\infty}^{\pm}) \\ c A_{\infty}^{\pm} (g^A)' (c P_{\infty}^{\pm}) & -h A_{\infty}^{\pm} \end{pmatrix}.$$

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

$$\det(J^{\pm}) = P_{\infty}^{\pm} A_{\infty}^{\pm} (hk - c^2(g^P)'(cA_{\infty}^{\pm})(g^A)'(cP_{\infty}^{\pm})) = -P_{\infty}^{\pm} A_{\infty}^{\pm}hk f'(A_{\infty}^{\pm}),$$

and its sign depends only on the sign of $f'(A_\infty^\pm)$. 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 $M_F([0, 1])$ such that, for all positive, bounded and continuous function $f$ on $[0, 1]$:

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

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

where

$$\psi(x, y) := c(x, y)\phi(x, y). \quad (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;
**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.3) has 2 positive equilibria 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 they are the two positive solutions to

$$ \left\{ \begin{array}{l} P_{\infty} = \frac{1}{k} g^P(cA_{\infty}) \\ A_{\infty} = \frac{1}{h} g^A \left( \frac{c}{k} g^P(cA_{\infty}) \right), \text{ and } cA_{\infty} \in (C_0^-, C_0^+) \end{array} \right. \qquad (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

$$ \begin{array}{l} g^P(cA_\infty) = kP_\infty \quad \text{or} \quad P_\infty = 0 \\ g^A(cP_\infty) = hA_\infty \quad \text{or} \quad A_\infty = 0. \end{array} \qquad (4.5) $$

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(\frac{c}{k}g^P(cx)) - x$ on $(\frac{C_0^-}{c}, \frac{C_0^+}{c})$. Since $g^A(c/k \times \cdot)$ and $g^P(c \times \cdot)$ are concave functions with $g^A$ non-decreasing, their composition is also concave and so is $f$. Since $f(C_0^\pm/c) = -d_A/h - C_0^\pm/c < 0$, it admits

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

* 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 equilibria, 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} -kP_{\infty}^{\pm} & cP_{\infty}^{\pm}(g^{P})'(cA_{\infty}^{\pm}) \\ cA_{\infty}^{\pm}(g^{A})'(cP_{\infty}^{\pm}) & -hA_{\infty}^{\pm} \end{pmatrix}. $$

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

$$ \det(J^{\pm}) = P_{\infty}^{\pm} A_{\infty}^{\pm} (hk - c^2(g^P)'(cA_{\infty}^{\pm})(g^A)'(cP_{\infty}^{\pm})) = -P_{\infty}^{\pm} A_{\infty}^{\pm} h k f'(A_{\infty}^{\pm}), $$

and its sign depends only on the sign of $f'(A_\infty^\pm)$. According to the previous study of the 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. ■
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

$$ \exists \bar{a}_0, \bar{p}_1 \in \mathbb{R}_+^*, \bar{p}_2 \in \mathbb{R}_+, \bar{x}_1, \bar{x}_2, \bar{y}_0 \in [0, 1], \quad \begin{cases} \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{cases} \\ \text{with} \quad \begin{cases} g^P(\bar{a}_0 \psi(\bar{x}_1, \bar{y}_0)) = g^P(\bar{a}_0 \psi(\bar{x}_2, \bar{y}_0)) = k(\bar{p}_1 + \bar{p}_2) \\ g^A(\bar{p}_1 \psi(\bar{x}_1, \bar{y}_0) + \bar{p}_2 \psi(\bar{x}_2, \bar{y}_0)) = h \bar{a}_0. \end{cases} \tag{4.9} $$

All these stationary states are unstable, except the state

$$ \begin{cases} \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{cases} \tag{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, \tag{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, \cdot) \bar{p}_t(x) dx$, $\int_0^1 \bar{p}_t(x) dx$, $\int_0^1 \psi(\cdot, 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 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$).
### 4.2.2 Behaviour of the kinetic equations (3.4)

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

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

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

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

where

$$\psi(x, y) := c(x, y)\phi(x, y). \quad (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.4) does not admit any 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

$$\exists \bar{a}_0, \bar{p}_1 \in \mathbb{R}_+^*, \bar{p}_2 \in \mathbb{R}_+, \bar{x}_1, \bar{x}_2, \bar{y}_0 \in [0, 1], \quad \begin{cases} \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{cases} \\ \text{with} \quad \begin{cases} g^P(\bar{a}_0 \psi(\bar{x}_1, \bar{y}_0)) = g^P(\bar{a}_0 \psi(\bar{x}_2, \bar{y}_0)) = k(\bar{p}_1 + \bar{p}_2) \\ g^A(\bar{p}_1 \psi(\bar{x}_1, \bar{y}_0) + \bar{p}_2 \psi(\bar{x}_2, \bar{y}_0)) = h \bar{a}_0. \end{cases} \quad (4.9)$$

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. \quad (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, \quad (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, \cdot) \bar{p}_t(x)dx$, $\int_0^1 \bar{p}_t(x)dx$, $\int_0^1 \psi(\cdot, 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).

In the proof of Proposition 4.5, we will see that the result is not restricted to the specified forms of $g^P$ and $g^A$ given in (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 the functions specified in (4.2).

We will also see that in general, System (3.4) has no stationary state with densities, all stationary states will be composed of Dirac measures. Moreover, the maximal number of Dirac measures found in
Proof. Any stationary state in $L^1([0, 1]^2)$ satisfies that $(\bar{P}^\infty, \bar{A}^\infty)$

$$ \left\{ \begin{array}{l} g^P \left( \int_0^1 \psi(x, y) \bar{A}^\infty(dy) \right) = \int_0^1 k \bar{P}^\infty(dx'), \quad \text{for all } x \in \text{supp} \bar{P}^\infty, \\ \\ g^A \left( \int_0^1 \psi(x, y) \bar{P}^\infty(dx') \right) = \int_0^1 h \bar{A}^\infty(dy'), \quad \text{for all } y \in \text{supp} \bar{A}^\infty \end{array} \right. \qquad (4.12) $$

Since $(\bar{P}^\infty, \bar{A}^\infty) \in L^1([0, 1]^2)$, $\int_0^1 k \bar{P}^\infty(dx')$ and $\int_0^1 h \bar{A}^\infty(dy')$ 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 $S$ the set of stationary states given by (4.9) excluding equilibrium (4.10). Let us prove that all states in $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 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 \ge 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(\bar{p}_t(x)) = g^P \left( \int_0^1 \psi(x, y) \bar{a}_t(y) dy \right) - \int_0^1 k \, \text{hel} \bar{p}_0(t(x')) dx'. $$

Under our assumption, $S$ admits at least a stable state and since the initial state is close to this set, we deduce that $(\bar{p}_t, \bar{a}_t)$ will converge to some $(\hat{P}^\infty, \hat{A}^\infty) \in 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]$,

$$ \begin{aligned} \frac{d}{dt} \log(\bar{p}_t(x)) &\xrightarrow{t \to \infty} g^P \left( \int_0^1 \psi(x, y) \hat{a}^\infty(y) dy \right) - \int_0^1 k \hat{p}^\infty(x') dx' \\ &= g^P \left( \int_0^1 \psi(x, y) \hat{a}^\infty(y) dy \right) - g^P \left( \int_0^1 \psi(\hat{x}_1, y) \hat{a}^\infty(y) dy \right). \end{aligned} $$

The latter quantity is positive as soon as $x \in (\hat{x}_1, x_0)$, which contradicts the fact that $\bar{p}_t(x)$ converges to 0 for all $x \notin \{\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 \le i, j \le N$. Then the continuous model (3.4) can be approximated using the rectangular rule by the following system of coupled ODEs:

$$ \begin{aligned} \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, \\ \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, \end{aligned} \qquad (4.13) $$

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
such a stationary state corresponds to the maximal number of points that can be found in the inverse image of a positive real number for the functions $g^P$ and $g^A$. The stability of these stationary states can then be deduced using the same kind of arguments as the ones in the following proof.

Proposition 4.5 shows that the only possible stable equilibrium is composed of only one plant species and one pollinator species. Simulations below will illustrate this Proposition. 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 interactions. 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 of whether it is possible, in this type of ecological model, to maintain a stable coexistence of several 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 $(\bar{P}^{\infty}, \bar{A}^{\infty})$ in $L^1([0, 1]^2)$ satisfies that

$$ \left\{ \begin{array}{l} g^P\left(\int_0^1 \psi(x, y) \bar{A}^\infty(dy)\right) = \int_0^1 k \bar{P}^\infty(dx'), \quad \text{for all } x \in \text{supp} \bar{P}^\infty, \\ \\ g^A\left(\int_0^1 \psi(x, y) \bar{P}^\infty(dx')\right) = \int_0^1 h \bar{A}^\infty(dy'), \quad \text{for all } y \in \text{supp} \bar{A}^\infty \end{array} \right. \qquad (4.12) $$

Since $(\bar{P}^{\infty}, \bar{A}^{\infty}) \in L^1([0, 1]^2)$, $\int_0^1 k \bar{P}^{\infty}(dx')$ and $\int_0^1 h \bar{A}^{\infty}(dy')$ 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 one 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 us 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 \ge 0$ and we can study $\log(\bar{p}_t(x))$ for all $t \in \mathbb{R}_+$ and all $x \in [0, 1]$:

$$ \frac{d}{dt} \log(\bar{p}_t(x)) = g^P \left( \int_0^1 \psi(x, y) \bar{a}_t(y) dy \right) - \int_0^1 k, \bar{p}_0 t(x') dx'. $$

Under our assumption, $\mathcal{S}$ admits at least a stable state and since the initial state is close to this set, we deduce that $(\bar{p}_t, \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]$,

$$ \begin{aligned} \frac{d}{dt} \log(\bar{p}_t(x)) & \underset{t \to \infty}{\longrightarrow} g^P \left( \int_0^1 \psi(x, y) \hat{a}^\infty(y) dy \right) - \int_0^1 k \hat{p}^\infty(x') dx' \\ &= g^P \left( \int_0^1 \psi(x, y) \hat{a}^\infty(y) dy \right) - g^P \left( \int_0^1 \psi(\hat{x}_1, y) \hat{a}^\infty(y) dy \right). \end{aligned} $$

The latter quantity is positive as soon as $x \in (\hat{x}_1, x_0)$, which contradicts the fact that $\bar{p}_t(x)$ converges to 0 for all $x \notin \{\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. ■
<figure><img src="image_5.png" /><figcaption>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$</figcaption></figure>

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$).

# 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
<figure><img src="image_6.png" /><figcaption>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$ species of plants and pollinators.</figcaption></figure>

Let us illustrate numerically Proposition 4.5 (see Figure 4). To do so, we have discretized the continuous model (3.5). Let $N$ be a positive integer and let $(p_i(t))_i$, resp. $(a_j(t))_j$ be approximations of $(\bar{p}_t(x_i))_i$, resp. $(\bar{a}(y_j))_j$, at $x_i = i/N$ and $y_j = j/N$, $0 \le i, j \le N$. The solution of (3.5) can be approximated using the rectangular rule by the following system of coupled ODEs:

$$ \begin{aligned} \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, \\ \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, \end{aligned} \quad (4.13) $$

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). In order to fit the hypotheses of Proposition 4.5 and ensure the convergence towards an explicit equilibrium of the dynamics, the interaction matrix $(c_{ij})_{ij}$ is non-decreasing 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 checked numerically the conclusion of Proposition 4.5: 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 species is moderately specialist ($x \approx 0.6$) whereas the remaining pollinator is genuinely generalist ($y = 1$).

# 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
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/

# 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)$:

$$ \langle \mathbf{P}_t^{K,n}, \mathbf{1} \rangle = \frac{1}{nK} \sum_{i=1}^{n} P_t^{K,i} \\ \leq \langle \tilde{\mathbf{P}}_t^{K,n}, \mathbf{1} \rangle = \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle + \int_0^t \int_E \frac{1}{nK} \sum_{i=1}^{n} \mathbf{1}_{i=k} \mathbf{1}_{\theta \leq M^P \tilde{P}_{s-}^{K,i}} Q_B^P(ds, dk, d\theta). $$

Since

$$ \mathbb{E} \left( \sup_{s \leq t \wedge \tau_N^K} \langle \mathbf{P}_s^{K,n}, \mathbf{1} \rangle \right) \leq 1 + \mathbb{E} \left( \sup_{s \leq t \wedge \tau_N^K} \langle \mathbf{P}_s^{K,n}, \mathbf{1} \rangle^2 \right) \leq 1 + \mathbb{E} \left( \sup_{s \leq t \wedge \tau_N^K} \langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle^2 \right), $$

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

$$ \tau_N^K = \inf \left\{ t \geq 0, \langle \tilde{\mathbf{P}}_t^{K,n}, \mathbf{1} \rangle \geq N \text{ or } \langle \tilde{\mathbf{A}}_t^{K,m}, \mathbf{1} \rangle \geq N \right\}. $$

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

$$ \begin{aligned} \langle \tilde{\mathbf{P}}_t^{K,n}, \mathbf{1} \rangle^2 &= \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 + \int_0^t \int_E \sum_{i=1}^n \mathbf{1}_{k=i} \left( \left( \langle \tilde{\mathbf{P}}_{s-}^{K,n}, \mathbf{1} \rangle + \frac{1}{nK} \right)^2 - \langle \tilde{\mathbf{P}}_{s-}^{K,n}, \mathbf{1} \rangle^2 \right) \mathbf{1}_{\theta \le M^P \tilde{P}_{s-}^{K,i}} Q_B^P(ds, dk, d\theta) \\ &= \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 + \int_0^t \int_E \sum_{i=1}^n \mathbf{1}_{k=i} \left( \frac{2}{nK} \sup_{u \le s-} \langle \tilde{\mathbf{P}}_u^{K,n}, \mathbf{1} \rangle + \frac{1}{n^2 K^2} \right) \mathbf{1}_{\theta \le M^P \tilde{P}_{s-}^{K,i}} Q_B^P(ds, dk, d\theta). \end{aligned} $$
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 addresses these questions by studying systems of ODEs where each 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 these 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 allow 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/

# 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.* The proof is divided into several steps.

**Step 1:** let us first prove some moment estimates, including (1.12). First, notice that the process can be stochastically bounded by a process $(\hat{P}_t^{K,n}, \hat{A}_t^{K,m})$ where only births occur at the maximal birth
Since the right hand side is an increasing process, we can again replace the left hand side by $\sup_{s \le t \wedge \tau_N^K} \langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle$. Then, taking the expectation gives:

$$ 
\begin{aligned}
\mathbb{E} \left( \sup_{s \le t \wedge \tau_N^K} \langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle^2 \right) &\le \mathbb{E} \left( \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 \right) + \mathbb{E} \left[ \int_0^{t \wedge \tau_N^K} \sum_{i=1}^n \left( \frac{2}{nK} \sup_{u \le s} \langle \tilde{\mathbf{P}}_u^{K,n}, \mathbf{1} \rangle + \frac{1}{n^2 K^2} \right) M^P \tilde{P}_s^{K,i} ds \right] \\
&\le \mathbb{E} \left( \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 \right) + \int_0^t M^P \left[ \left( 2 + \frac{1}{nK} \right) \mathbb{E} \left( 1 + \sup_{u \le s \wedge \tau_N^K} \langle \tilde{\mathbf{P}}_u^{K,n}, \mathbf{1} \rangle^2 \right) \right] ds.
\end{aligned}
$$

Using Gronwall's lemma, we deduce

$$ \mathbb{E} \left( \sup_{s \le t \wedge \tau_N^K} \langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle^2 \right) \le \left[ \mathbb{E} \left( \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 \right) + 3M^P \right] e^{3M^Pt} < +\infty. $$

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

$$ \mathbb{E} \left( \sup_{s \le t} \langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle^2 \right) \le \left[ \mathbb{E} \left( \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 \right) + 3M^P \right] e^{3M^Pt} < +\infty. $$

Similarly,

$$ \mathbb{E} \left( \sup_{s \le t} \langle \tilde{\mathbf{A}}_s^{K,m}, \mathbf{1} \rangle^2 \right) \le \left[ \mathbb{E} \left( \langle \mathbf{A}_0^{K,m}, \mathbf{1} \rangle^2 \right) + 3M^A \right] e^{3M^At} < +\infty. \quad (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:

$$ 
\begin{aligned}
M_t^{K,i} = & \int_0^t \int_E \mathbf{1}_{i=k} \mathbf{1}_{\theta \le b^P} (\sum_{j \sim i} R_{s-}^{A,K,ij}) P_{s-}^{K,i} Q_B^P(ds, dk, d\theta) - \int_0^t b^P (\sum_{j \sim i} R_s^{K,ij}) P_s^{K,i} ds \\
& - \int_0^t \int_E \mathbf{1}_{i=k} \mathbf{1}_{\theta \le [d^P (\sum_{j \sim i} R_{s-}^{A,K,ij}) - k \star \mathbf{P}_{s-}^{K,n}(x^i)]} P_{s-}^{K,i} Q_D^P(ds, dk, d\theta) \\
& + \int_0^t [d^P (\sum_{j \sim i} R_s^{A,K,ij}) - k \star \mathbf{P}_s^{K,n}(x^i)] P_s^{K,i} ds
\end{aligned}
\quad (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 \ge 1$, we have:

$$ 
\begin{aligned}
\int_0^1 f(x) d\tilde{\mathbf{P}}_t^{n,m}(x) = & \int_0^1 f(x) d\tilde{\mathbf{P}}_0^{n,m}(x) dx \\
& + \int_0^t \frac{1}{n} \sum_{i=1}^n f(x^i) \left[ g^P \left( \frac{1}{m} \sum_{j=1}^m mc_{i,j}^{n,m} G_{ij} \tilde{A}_s^j \right) - \frac{1}{n} \sum_{\ell=1}^n k(x^i, x^\ell) \tilde{P}_s^\ell \right] \tilde{P}_s^i ds,
\end{aligned}
\quad (A.3) $$

$$ 
\begin{aligned}
\int_0^1 f(y) d\tilde{\mathbf{A}}_t^{n,m}(y) = & \int_0^1 f(y) d\tilde{\mathbf{A}}_0^{n,m}(y) dy \\
& + \int_0^t \frac{1}{m} \sum_{j=1}^m f(y^j) \left[ g^A \left( \frac{1}{n} \sum_{i=1}^n nc_{i,j}^{n,m} G_{ij} \tilde{P}_s^i \right) - \frac{1}{m} \sum_{\ell=1}^m h(y^j, y^\ell) \tilde{A}_s^\ell \right] \tilde{A}_s^j ds.
\end{aligned}
$$
rate, denoted by $M^P := \sup_{r \in \mathbb{R}^+} b^P(r)$:

$$ \langle \mathbf{P}_t^{K,n}, \mathbf{1} \rangle = \frac{1}{nK} \sum_{i=1}^{n} P_t^{K,i} \\ \leq \langle \tilde{\mathbf{P}}_t^{K,n}, \mathbf{1} \rangle = \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle + \int_0^t \int_E \frac{1}{nK} \sum_{i=1}^{n} \mathbf{1}_{i=k} \mathbf{1}_{\theta \leq M^P \tilde{P}_{s-}^{K,i}} Q_B^P(ds, dk, d\theta). $$

Since

$$ \mathbb{E} \left( \sup_{s \le t \wedge \tau_N^K} \langle \mathbf{P}_s^{K,n}, \mathbf{1} \rangle \right) \le 1 + \mathbb{E} \left( \sup_{s \le t \wedge \tau_N^K} \langle \mathbf{P}_s^{K,n}, \mathbf{1} \rangle^2 \right) \le 1 + \mathbb{E} \left( \sup_{s \le t \wedge \tau_N^K} \langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle^2 \right), $$

we only have to find a bound for the second moment of $\langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle$.

For a constant $N > 0$, we introduce the stopping time

$$ \tau_N^K = \inf \left\{ t \ge 0, \langle \tilde{\mathbf{P}}_t^{K,n}, \mathbf{1} \rangle \ge N \text{ or } \langle \tilde{\mathbf{A}}_t^{K,m}, \mathbf{1} \rangle \ge N \right\}. $$

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

$$ \begin{aligned} \langle \tilde{\mathbf{P}}_t^{K,n}, \mathbf{1} \rangle^2 &= \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 + \int_0^t \int_E \sum_{i=1}^n \mathbf{1}_{k=i} \left( \left( \langle \tilde{\mathbf{P}}_{s-}^{K,n}, \mathbf{1} \rangle + \frac{1}{nK} \right)^2 - \langle \tilde{\mathbf{P}}_{s-}^{K,n}, \mathbf{1} \rangle^2 \right) \mathbf{1}_{\theta \le M^P \tilde{P}_{s-}^{K,i}} Q_B^P(ds, dk, d\theta) \\ &= \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 + \int_0^t \int_E \sum_{i=1}^n \mathbf{1}_{k=i} \left( \frac{2}{nK} \sup_{u \le s-} \langle \tilde{\mathbf{P}}_u^{K,n}, \mathbf{1} \rangle + \frac{1}{n^2 K^2} \right) \mathbf{1}_{\theta \le M^P \tilde{P}_{s-}^{K,i}} Q_B^P(ds, dk, d\theta). \end{aligned} $$

Since the right hand side is an increasing process, we can replace the left hand side by $\sup_{s \le t \wedge \tau_N^K} \langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle^2$.
Then, taking the expectation gives:

$$ \begin{aligned} \mathbb{E} \left( \sup_{s \le t \wedge \tau_N^K} \langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle^2 \right) &\le \mathbb{E} \left( \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 \right) + \mathbb{E} \left[ \int_0^{t \wedge \tau_N^K} \sum_{i=1}^n \left( \frac{2}{nK} \sup_{u \le s} \langle \tilde{\mathbf{P}}_u^{K,n}, \mathbf{1} \rangle + \frac{1}{n^2 K^2} \right) M^P \tilde{P}_s^{K,i} \, ds \right] \\ &\le \mathbb{E} \left( \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 \right) + \int_0^t M^P \left[ \left( 2 + \frac{1}{nK} \right) \mathbb{E} \left( 1 + \sup_{u \le s \wedge \tau_N^K} \langle \tilde{\mathbf{P}}_u^{K,n}, \mathbf{1} \rangle^2 \right) \right] \, ds. \end{aligned} $$

Using Gronwall's lemma, we deduce

$$ \mathbb{E} \left( \sup_{s \le t \wedge \tau_N^K} \langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle^2 \right) \le \left[ \mathbb{E} \left( \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 \right) + 3M^P \right] e^{3M^Pt} < +\infty. $$

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

$$ \mathbb{E} \left( \sup_{s \le t} \langle \tilde{\mathbf{P}}_s^{K,n}, \mathbf{1} \rangle^2 \right) \le \left[ \mathbb{E} \left( \langle \mathbf{P}_0^{K,n}, \mathbf{1} \rangle^2 \right) + 3M^P \right] e^{3M^Pt} < +\infty. $$

Similarly,

$$ \mathbb{E} \left( \sup_{s \le t} \langle \tilde{\mathbf{A}}_s^{K,m}, \mathbf{1} \rangle^2 \right) \le \left[ \mathbb{E} \left( \langle \mathbf{A}_0^{K,m}, \mathbf{1} \rangle^2 \right) + 3M^A \right] e^{3M^At} < +\infty. \quad (A.1) $$

Notice that the upper bounds 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].
where conditionally on $(x^i, y^j)_{1\le i\le n, 1\le j\le m}$, $G_{ij} \sim \text{Bern}(\phi(x^i, y^j))$ and are independent random variables, and $c_{i,j}^{n,m}$ 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.,

$$ \langle \tilde{\mathbf{P}}_t^{n,m}, 1 \rangle \le \langle \tilde{\mathbf{P}}_0^{n,m}, 1 \rangle e^{M_P t} \quad \text{and} \quad \langle \tilde{\mathbf{A}}_t^{n,m}, 1 \rangle \le \langle \tilde{\mathbf{A}}_0^{n,m}, 1 \rangle e^{M_A t}, \quad (A.4) $$

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

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

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

$$ C(T) = \sup_{n,m} \left( \langle \tilde{\mathbf{P}}_0^{n,m}, 1 \rangle e^{M_P T} + \langle \tilde{\mathbf{A}}_0^{n,m}, 1 \rangle e^{M_A T} \right) < +\infty \quad \text{a.s..} \quad (A.5) $$

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

**Step 2: Relative compactness in $C([0, T], M_{\le C(T)}([0,1])^2)$**

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}} (\langle \tilde{\mathbf{P}}_0^{n,m}, 1 \rangle + \langle \tilde{\mathbf{A}}_0^{n,m}, 1 \rangle) \text{ is bounded almost surely, and} $$

$$ (ii) \text{ for any continuous function } f \text{ on } [0, 1], \text{ and for any } \varepsilon > 0, $$

$$ \lim_{\delta \to 0} \limsup_{n,m \in \mathbb{N}} \mathbb{P} \left( \sup_{|t-s| \le \delta} |\langle \tilde{\mathbf{P}}_s^{n,m}, f \rangle - \langle \tilde{\mathbf{P}}_t^{n,m}, f \rangle| \ge \varepsilon \right) = 0 $$

$$ \lim_{\delta \to 0} \limsup_{n,m \in \mathbb{N}} \mathbb{P} \left( \sup_{|t-s| \le \delta} |\langle \tilde{\mathbf{A}}_s^{n,m}, f \rangle - \langle \tilde{\mathbf{A}}_t^{n,m}, f \rangle| \ge \varepsilon \right) = 0 $$

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

$$ \mathbb{E}[\tilde{R}^{A,n,m}] = \mathbb{E} \left[ \frac{1}{m} \sum_{j=1}^{m} mc_{i,j}^{n,m} G_{ij} \tilde{A}_s^j \right] \le \frac{1}{m} \sum_{j=1}^{m} \mathbb{E}[mc_{i,j}^{n,m}] C(T) \le \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}(\tilde{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 Lypschiitz function). Thus, on $\{\tilde{R}^{A,n,m} \le \tilde{C}\}$, for any positive continuous function $f$ on $[0, 1]$ and any $0 \le s \le t \le T$,

$$ |\langle \tilde{\mathbf{P}}_t^{n,m}, f \rangle - \langle \tilde{\mathbf{P}}_s^{n,m}, f \rangle| \le \int_s^t (\tilde{M} + \|k\|_{\infty} C(T)) \|f\|_{\infty} C(T) du. $$

Same computations can be done for the sequence of processes $(\tilde{\mathbf{A}}_n^{n,m})_{n \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 $(\tilde{\mathbf{P}}_n^{n,m}, \tilde{\mathbf{A}}_n^{n,m})$, denoted here by $(\tilde{\mathbf{P}}, \tilde{\mathbf{A}})$, is a process whose time marginals at $t > 0$ admit densities $\tilde{p}_t$ and $\tilde{a}_t$ with
**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:

$$ M_t^{K,i} = \int_0^t \int_E \mathbf{1}_{i=k} \mathbf{1}_{\theta \le b^P(\sum_{j \sim i} R_{s-}^{A,K,ij})} P_{s-}^{K,i} Q_B^P(ds, dk, d\theta) - \int_0^t b^P(\sum_{j \sim i} R_s^{K,ij}) P_s^{K,i} ds \\ - \int_0^t \int_E \mathbf{1}_{i=k} \mathbf{1}_{\theta \le [d^P(\sum_{j \sim i} R_{s-}^{A,K,ij}) - k \star \mathbf{P}_{s-}^{K,n}(x^i)]} P_{s-}^{K,i} Q_D^P(ds, dk, d\theta) \\ + \int_0^t [d^P(\sum_{j \sim i} R_s^{A,K,ij}) - k \star \mathbf{P}_s^{K,n}(x^i)] P_s^{K,i} ds \quad (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 \ge 1$, we have:

$$ \int_0^1 f(x) d\tilde{\mathbf{P}}_t^{n,m}(x) = \int_0^1 f(x) d\tilde{\mathbf{P}}_0^{n,m}(x) \\ + \int_0^1 \frac{1}{n} \sum_{i=1}^n f(x^i) \left[ g^P \left( \frac{1}{m} \sum_{j=1}^m mc_{i,j}^{n,m} G_{ij} \tilde{A}_s^j \right) - \frac{1}{n} \sum_{\ell=1}^n k(x^i, x^\ell) \tilde{P}_s^\ell \right] \tilde{P}_s^i ds, \quad (A.3) $$

$$ \int_0^1 f(y) d\tilde{\mathbf{A}}_t^{n,m}(y) = \int_0^1 f(y) d\tilde{\mathbf{A}}_0^{n,m}(y) \\ + \int_0^1 \frac{1}{m} \sum_{j=1}^m f(y^j) \left[ g^A \left( \frac{1}{n} \sum_{i=1}^n nc_{i,j}^{n,m} G_{ij} \tilde{P}_s^i \right) - \frac{1}{m} \sum_{\ell=1}^m h(y^j, y^\ell) \tilde{A}_s^\ell \right] \tilde{A}_s^j ds. $$

where conditionally on $(x^i, y^j)_{i \in [1,n], j \in [1,m]}$, $G_{ij} \sim \text{Bern}(\phi(x^i, y^j))$ and are independent r.v., and $c_{i,j}^{n,m}$ are independent r.v. 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.,

$$ \langle \tilde{\mathbf{P}}_t^{n,m}, 1 \rangle \le \langle \tilde{\mathbf{P}}_0^{n,m}, 1 \rangle e^{M_P t} \quad \text{and} \quad \langle \tilde{\mathbf{A}}_t^{n,m}, 1 \rangle \le \langle \tilde{\mathbf{A}}_0^{n,m}, 1 \rangle e^{M_A t}, \quad (A.4) $$

for all $n$ and $m$. Notice that these upper-bounds hold whatever the values of $G_{ij}$ and $c_{ij}^{n,m}$. By Assumption (3.3), we can define a constant $C(T) < +\infty$ such that almost surely

$$ \sup_{n,m} \left( \langle \tilde{\mathbf{P}}_0^{n,m}, 1 \rangle e^{M_P T} + \langle \tilde{\mathbf{A}}_0^{n,m}, 1 \rangle e^{M_A T} \right) < C(T). \quad (A.5) $$

Thus, the processes $(\tilde{\mathbf{P}}_t^{n,m}(dx), \tilde{\mathbf{A}}_t^{n,m}(dy))_{t \in [0,T]}$ take a.s. their values in $(M_{\le C(T)}([0,1]))^2$ which is a compact set.

**Step 2: Relative compactness in $C([0, T], M_{\le C(T)}([0, 1]))$**

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

$$ (i) \sup_{n,m \in \mathbb{N}} \left( \langle \tilde{\mathbf{P}}_0^{n,m}, 1 \rangle + \langle \tilde{\mathbf{A}}_0^{n,m}, 1 \rangle \right) \text{ is bounded almost surely, and} $$
respect to the Lebesgue measure on $[0, 1]$. To achieve this, we dominate the measures $\tilde{\mathbf{P}}_t^{n,m}(dx)$ and $\tilde{\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

$$ \langle \bar{\mathbf{P}}_t, f \rangle \le \langle \bar{\mathbf{P}}_0, f \rangle e^{M_P t} = \int_0^1 f(x) (p_0(x) e^{M_P t}) dx, \quad \text{a.s.} $$

$$ \langle \bar{\mathbf{A}}_t, f \rangle \le \langle \bar{\mathbf{A}}_0, f \rangle e^{M_A t} = \int_0^1 f(y) (a_0(y) e^{M_A t}) 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 $C([0, 1], M_{\le 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) [g^P ((c\phi) \star \mathbf{A}_s(x)) - k \star \mathbf{P}_s(x)] d\mathbf{P}_s(x) ds \\ \langle \mathbf{A}_t, f \rangle - \langle \mathbf{A}_0, f \rangle - \int_0^t \int_0^1 f(x) [g^A ((c\phi) \star \mathbf{P}_s(x)) - h \star \mathbf{A}_s(x)] d\mathbf{A}_s(x) ds \end{pmatrix}. $$

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

$$ W_1(\nu, \mu) = \sup \left\{ \int_0^1 f(x) d(\nu - \mu)(x) \middle| f : [0, 1] \mapsto \mathbb{R}^+ \text{ Lipschitz, } Lip(f) \le 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 $(\mathbf{P}, \bar{\mathbf{A}})$ is a weak limiting value of $(\tilde{\mathbf{P}}^{n,m}, \tilde{\mathbf{A}}^{n,m})$, we obtain by definition of the weak convergence that

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

On the other hand,

$$ \mathbb{E} \left[ \| \psi(\tilde{\mathbf{P}}^{n,m}, \tilde{\mathbf{A}}^{n,m}) \|_1 \right] = \mathbb{E} \left[ \left\| \int_0^t \int_0^1 f(x) \left( g^P \left( \frac{1}{m} \sum_{j=1}^m G_{ij} m c_{i,j}^{n,m} \tilde{A}_s^j \right) - g^P \left( (c\phi) \star \tilde{\mathbf{A}}^{n,m} \right) \right) d\tilde{\mathbf{P}}_s^{n,m}(x) ds \right\| \right] \quad (A.7) $$

$$ + \mathbb{E} \left[ \left\| \int_0^t \int_0^1 f(x) \left( g^A \left( \frac{1}{n} \sum_{j=1}^n G_{ij} n c_{i,j}^{n,m} \tilde{P}_s^j \right) - g^A \left( (c\phi) \star \tilde{\mathbf{P}}^{n,m} \right) \right) d\tilde{\mathbf{A}}_s^{n,m}(x) ds \right\| \right]. \quad (A.8) $$

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

$$ \begin{aligned} (A.7) \le & \int_0^t \frac{C}{n} \sum_{i=1}^n \left\| \mathbb{E} \left| \frac{1}{m} \sum_{j=1}^m \tilde{A}_s^j \left( G_{ij} (mc_{i,j}^{n,m} - c(x^i, y^j)) + c(x^i, y^j) (G_{ij} - \phi(x^i, y^j)) \right) \right| \right\| ds \\ \le & C \left[ TC(T) \|mc^{n,m} - c\|_{\infty} \right. \\ & \left. + \int_0^t \frac{1}{n} \sum_{i=1}^n \mathbb{E} \left| \frac{1}{m} \sum_{j=1}^m \tilde{A}_s^j G_{ij} m (c_{i,j}^{n,m} - c^{n,m}(x^i, y^j)) + \tilde{A}_s^j c(x^i, y^j) (G_{ij} - \phi(x^i, y^j)) \right| \right], \end{aligned} $$
(ii) for any continuous function $f$ on $[0, 1]$, and for any $\varepsilon > 0$,

$$ \lim_{\delta \to 0} \limsup_{n,m \in \mathbb{N}} \mathbb{P} \left( \sup_{|t-s| \le \delta} |\langle \tilde{\mathbf{P}}_s^{n,m}, f \rangle - \langle \tilde{\mathbf{P}}_t^{n,m}, f \rangle| \ge \varepsilon \right) = 0 $$

$$ \lim_{\delta \to 0} \limsup_{n,m \in \mathbb{N}} \mathbb{P} \left( \sup_{|t-s| \le \delta} |\langle \tilde{\mathbf{A}}_s^{n,m}, f \rangle - \langle \tilde{\mathbf{A}}_t^{n,m}, f \rangle| \ge \varepsilon \right) = 0. $$

Point (i) corresponds to (3.3). Let us turn to point (ii). Recall the definition of $\tilde{R}_t^{A,i}$ in Proposition 2.1 and notice that $\tilde{R}_t^{A,i} = \frac{1}{m} \sum_{j=1}^m mc_{i,j}^{n,m} G_{ij} \tilde{A}_t^j$. Then,

$$ \mathbb{E}[\tilde{R}_t^{A,i}] \le \frac{1}{m} \sum_{j=1}^{m} \mathbb{E}[mc_{i,j}^{n,m}] C(T) \le \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}(\tilde{R}^{A,i} > \tilde{C})$ is sufficiently small. Then $g^P$ is bounded on $[0, \tilde{C}]$ by some constant $\tilde{M}$ (since it is a Lipschitz function). Thus, on $\{\tilde{R}^{A,i} \le \tilde{C}\}$, for any positive continuous function $f$ on $[0, 1]$ and any $0 \le s \le t \le T$,

$$ |\langle \tilde{\mathbf{P}}_t^{n,m}, f \rangle - \langle \tilde{\mathbf{P}}_s^{n,m}, f \rangle| \le \int_s^t (\tilde{M} + \|k\|_{\infty} C(T)) \|f\|_{\infty} C(T) du. $$

Same computations can be done for the sequence of processes $(\tilde{\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 value of $(\tilde{\mathbf{P}}^{n,m}, \tilde{\mathbf{A}}^{n,m})$, denoted here by $(\tilde{\mathbf{P}}, \tilde{\mathbf{A}})$, is a process whose time marginals at $t > 0$ admit densities $\tilde{p}_t$ and $\tilde{a}_t$ with respect to the Lebesgue measure on $[0, 1]$. To achieve this, we dominate the measures $\tilde{\mathbf{P}}_t^{n,m}(dx)$ and $\tilde{\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

$$ \langle \tilde{\mathbf{P}}_t, f \rangle \le \langle \tilde{\mathbf{P}}_0, f \rangle e^{M_P t} = \int_0^1 f(x) (p_0(x) e^{M_P t}) dx, \quad \text{a.s.} $$

$$ \langle \tilde{\mathbf{A}}_t, f \rangle \le \langle \tilde{\mathbf{A}}_0, f \rangle e^{M_A t} = \int_0^1 f(y) (a_0(y) e^{M_A t}) 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 for all $t \in [0, 1]$ the map $\psi_t$ on $C([0, 1], \mathcal{M}_{\le C(T)}([0, 1])^2)$, endowed with the uniform convergence, defined by

$$ \psi_t(\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) [g^P((c\phi) \star \mathbf{A}_s(x)) - k \star \mathbf{P}_s(x)] d\mathbf{P}_s(x) ds \\ \langle \mathbf{A}_t, f \rangle - \langle \mathbf{A}_0, f \rangle - \int_0^t \int_0^1 f(x) [g^A((c\phi) \star \mathbf{P}_s(x)) - h \star \mathbf{A}_s(x)] d\mathbf{A}_s(x) ds \end{pmatrix}. $$

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

$$ \mathcal{W}_1(\nu, \mu) = \sup \left\{ \int_0^1 f(x) d(\nu - \mu)(x) \middle| f : [0, 1] \mapsto [0, 1] \text{ Lipschitz, } Lip(f) \le 1 \right\}, $$
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 \le n, j \le m})$, we find for $n, m$ sufficiently large,

$$
\begin{align*}
(A.7) & \le \varepsilon + \int_0^t \frac{C}{n} \sum_{i=1}^n \mathbb{E} \left[ \mathrm{Var}_{x,y} \left( \frac{1}{m} \sum_{j=1}^m \tilde{A}_s^j G_{ij} mc_{i,j}^{n,m} \right) \right]^{1/2} + \mathbb{E} \left[ \mathrm{Var}_{x,y} \left( \frac{1}{m} \sum_{j=1}^m \tilde{A}_s^j c(x^i, y^j) G_{ij} \right) \right]^{1/2} \\
& \le \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} \\
& \le \varepsilon + C T \frac{C(T)(V_{max}^{1/2} + \|c\|_{\infty}\|\phi(1-\phi)\|_{\infty}^{1/2})}{m^{1/2}} \\
& \le 2\varepsilon \text{ for } n, m \text{ sufficiently large.}
\end{align*}
$$

The second term (A.8) can be treated similarly. This finally proves that $\mathbb{E}[\|\psi(\tilde{\mathbf{P}}^{n,m}, \tilde{\mathbf{A}}^{n,m})\|_1]$ 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 $M_F([0,1])^2$:

$$ \rho(\nu, \mu) = \sup \left\{ \int_0^1 f(x)d(\nu - \mu)(x) \middle| 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

$$
\begin{align*}
|\langle \mathbf{P}_t^1 - \mathbf{P}_t^2, f \rangle| &= \left| \int_0^t \int_0^1 f(x) [g^P(c\phi \star \mathbf{A}_s^1(x)) - k \star \mathbf{P}_s^1] d\mathbf{P}_s^1(x) ds \right. \\
&\qquad \left. - \int_0^t \int_0^1 f(x) [g^P(c\phi \star \mathbf{A}_s^2(x)) - k \star \mathbf{P}_s^2] d\mathbf{P}_s^2(x) ds \right| \\
&\le C_1 \int_0^t \rho(\mathbf{P}_s^1, \mathbf{P}_s^2) ds + \int_0^t \int_0^1 (L^P |c\phi \star (\mathbf{A}_s^1 - \mathbf{A}_s^2(x))| + |k \star (\mathbf{P}_s^1 - \mathbf{P}_s^2)(x)|) d\mathbf{P}_s^2(x) ds \\
&\le C_2 \int_0^t (\rho(\mathbf{P}_s^1, \mathbf{P}_s^2) + \rho(\mathbf{A}_s^1, \mathbf{A}_s^2)) ds,
\end{align*}
$$

since $x \mapsto f(x)(g^P(c\phi \star \mathbf{A}_s^1(x)) - k \star \mathbf{P}_s^1(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}_t^1, \mathbf{P}_t^2) + \rho(\mathbf{A}_t^1, \mathbf{A}_t^2) \le C_2 \int_0^t (\rho(\mathbf{P}_s^1, \mathbf{P}_s^2) + \rho(\mathbf{A}_s^1, \mathbf{A}_s^2)) ds. $$

We conclude with Gronwall inequality that $\rho(\mathbf{P}_t^1, \mathbf{P}_t^2) + \rho(\mathbf{A}_t^1, \mathbf{A}_t^2) = 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.

<figure><img src="image_6.png" /></figure>
since [0, 1] is a compact space. The function $\psi_t$ is also bounded since all measures are bounded by $C(T)$.

Consider a weak limiting value $(\bar{\mathbf{P}}, \bar{\mathbf{A}})$ of $(\tilde{\mathbf{P}}^{n,m}, \tilde{\mathbf{A}}^{n,m})$, and denote by $(\tilde{\mathbf{P}}^{n,m}, \tilde{\mathbf{A}}^{n,m})$ again (with an abuse of notation) the subsequence converging to $(\bar{\mathbf{P}}, \bar{\mathbf{A}})$. Using the Skorokhod representation theorem, we can assume that these processes are defined on the same probability space and that the convergence is almost sure:

$$ \lim_{n,m \to +\infty} \sup_{t \le T} (W_1(\tilde{\mathbf{P}}_t^{n,m}, \bar{\mathbf{P}}_t) + W_1(\tilde{\mathbf{A}}_t^{n,m}, \bar{\mathbf{A}}_t)) = 0, \quad \text{a.s.} \qquad (A.6) $$

Then, by the properties of $\psi_t$,

$$ \mathbb{E} \left[ \| \psi_t(\tilde{\mathbf{P}}^{n,m}, \tilde{\mathbf{A}}^{n,m}) \|_1 \right] \underset{n,m \to +\infty}{\longrightarrow} \mathbb{E} \left[ \| \psi_t(\bar{\mathbf{P}}, \bar{\mathbf{A}}) \|_1 \right]. \qquad (A.7) $$

Our purpose is now to prove that for all $t \in [0, T]$,

$$ \lim_{n,m \to +\infty} \mathbb{E} \left[ \| \psi_t(\tilde{\mathbf{P}}^{n,m}, \tilde{\mathbf{A}}^{n,m}) \|_1 \right] = 0. \qquad (A.8) $$

This will ensure, in addition with (A.7), that $\psi(\bar{\mathbf{P}}, \bar{\mathbf{A}}) = (0, 0)$ a.s. In other words, the limiting value satisfies Equation (3.4) a.s.

By (A.3), (A.4) and (A.5) and denoting by $L^P$ (resp. $L^A$) the Lipschitz constant of $g^P$ (resp. $g^A$):

$$ \begin{aligned} & \mathbb{E} \left[ \| \psi_t(\tilde{\mathbf{P}}^{n,m}, \tilde{\mathbf{A}}^{n,m}) \|_1 \right] \\ &= \mathbb{E} \left[ \left| \int_0^t \frac{1}{n} \sum_{i=1}^n f(x^i) \tilde{P}_s^i \left( g^P \left( \frac{1}{m} \sum_{j=1}^m G_{ij} m c_{i,j}^{n,m} \tilde{A}_s^j \right) - g^P \left( (c\phi) \star \tilde{\mathbf{A}}_s^{n,m}(x^i) \right) \right) ds \right| \right] \\ & \quad + \mathbb{E} \left[ \left| \int_0^t \frac{1}{m} \sum_{j=1}^m f(y^j) \tilde{A}_s^j \left( g^A \left( \frac{1}{n} \sum_{i=1}^n G_{ij} n c_{i,j}^{n,m} \tilde{P}_s^i \right) - g^A \left( (c\phi) \star \tilde{\mathbf{P}}_s^{n,m}(y^j) \right) \right) ds \right| \right], \\ & \le \int_0^t \frac{C(T)L^P \|f\|_\infty}{n} \sum_{i=1}^n \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^m G_{ij} m c_{i,j}^{n,m} \tilde{A}_s^j - (c\phi) \star \tilde{\mathbf{A}}_s^{n,m}(x^i) \right| \right] \\ & \quad + \int_0^t \frac{C(T)L^A \|f\|_\infty}{m} \sum_{j=1}^m \mathbb{E} \left[ \left| \frac{1}{n} \sum_{i=1}^n G_{ij} n c_{i,j}^{n,m} \tilde{P}_s^i - (c\phi) \star \tilde{\mathbf{P}}_s^{n,m}(y^j) \right| \right] ds. \quad (A.9) \end{aligned} $$

The two terms in (A.9) having similar expressions, we focus on the first one with the aim of showing that uniformly for all $t \in [0, T]$,

$$ \lim_{n,m \to +\infty} \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} \tilde{A}_t^j - (c\phi) \star \tilde{\mathbf{A}}_t^{n,m}(x^i) \right| \right] = 0. \qquad (A.10) $$

The proof of (A.10) is itself divided into smaller steps. To this aim, let us set $\varepsilon > 0$ and show that for $n$ and $m$ sufficiently large the term in the left hand side of (A.10) is smaller than $\varepsilon$. The difficulty is that $(\tilde{A}_t^j)$ are not independent from $(G_{ij})$ and $(c_{ij}^{n,m})$, except at time $t = 0$.
## 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ï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.
**Step 4.1:** (A.10) is satisfied for $t=0$. At time $t=0$, we have:

$$ 
\begin{aligned}
& \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} \tilde{A}_{0}^{j} - (c\phi) \star \tilde{\mathbf{A}}_{0}^{n,m}(x^{i}) \right| \right] \\
& \leq \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^{m} \left( G_{ij} m c_{i,j}^{n,m} - m c^{n,m}(x^{i}, y^{j}) \phi(x^{i}, y^{j}) \right) \tilde{A}_{0}^{j} \right| \right] \\
& \quad + \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^{m} \phi(x^{i}, y^{j}) \tilde{A}_{0}^{j} \left( m c^{n,m}(x^{i}, y^{j}) - c(x^{i}, y^{j}) \right) \right| \right]. \quad (A.11)
\end{aligned}
$$

By (A.5), the second term in the right hand side of (A.11) is upper bounded by $C(T)\|\phi\|_{\infty}\|mc^{n,m}-c\|_{\infty}$, which converges to zero by Assumption 3.1.

Let us now consider the first term in the right hand side of (A.11). Since at $t=0$ the random variables $(\tilde{A}_0^j)$, $(G_{ij})$ and $(c_{ij}^{n,m})$ are independent, we obtain by Assumption 1.1, the assumptions of Proposition 3.2 and (1.3) and by using Cauchy-Schwarz inequality that

$$ 
\begin{aligned}
& \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^{m} \left( G_{ij} m c_{i,j}^{n,m} - m c^{n,m}(x^i, y^j) \phi(x^i, y^j) \right) \tilde{A}_0^j \right| \right] \\
& \leq \frac{1}{m} \mathbb{E} \left[ \mathrm{Var}^{1/2} \left( \sum_{j=1}^{m} \tilde{A}_0^j G_{ij} m c_{i,j}^{n,m} \mid x^i, y^j, \tilde{A}_0^j, i \in [1, n], j \in [1, m] \right) \right].
\end{aligned}
$$

Then, recalling (1.3) and the fact that $\mathrm{Var}(XY) = \mathrm{Var}(X)\mathbb{E}[Y^2] + \mathrm{Var}(Y)\mathbb{E}[X]^2$ for any independent r.v. $X$ and $Y$:

$$ 
\begin{aligned}
& \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^{m} \left( G_{ij} m c_{i,j}^{n,m} - m c^{n,m}(x^i, y^j) \phi(x^i, y^j) \right) \tilde{A}_0^j \right| \right] \\
& \leq \frac{C(T)}{m} \mathbb{E} \left[ \left( \sum_{j=1}^{m} \mathrm{Var}_{x,y} \left( G_{ij} m c_{i,j}^{n,m} \right) \right)^{1/2} \right] \\
& \leq \frac{C(T)}{\sqrt{m}} \left( V_{\mathrm{max}}^{1/2} + \|mc^{n,m}\|_{\infty} \|\phi(1-\phi)\|_{\infty}^{1/2} \right).
\end{aligned}
$$

Thus, the expectation in the l.h.s of (A.11) can be made smaller than $\epsilon/3$ for $n,m$ larger than some $N_0, M_0$.

Along the proof of Step 4.1, the independence between the random network $(G_{ij}$'s and $c_{ij}^{n,m}$'s) and the abundancies $((\tilde{A}_0^j)$ and $(\tilde{P}_0^i)$'s) was a crucial ingredient, which does not hold any more for $t > 0$. Hence the difficulty. We will first show that (A.8) holds on a small time intervals $[0, \delta)$, for some $\delta > 0$. For this, we will construct auxiliary processes $(\tilde{A}, \tilde{P})$ that do not depend on the $G_{ij}$'s and $c_{ij}^{n,m}$'s. Then, for all $i \in [1, n]$,

$$ 
\begin{aligned}
& \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} \tilde{A}_t^j - (c\phi) \star \tilde{\mathbf{A}}_t^{n,m}(x^i) \right| \right] \\
& \leq \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} (\tilde{A}_t^j - \hat{A}_t^j) \right| \right. \\
& \quad \left. + \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{ij}^{n,m} \hat{A}_t^j - c\phi \star \hat{\mathbf{A}}_t^{n,m}(x_i) \right| + \left| c\phi \star (\hat{\mathbf{A}}_t^{n,m} - \tilde{\mathbf{A}}_t^{n,m})(x_i) \right| \right]. \quad (A.12)
\end{aligned}
$$

The second term of (A.12) will be treated using the exact same computation as above and thus it will be upper-bounded by $\epsilon/3$ as soon as $n,m$ are greater than some threshold $N_0, M_0$ and if the auxiliary
[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í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.
processes satisfy (A.5). Our aim is thus to find auxiliary processes such that the latter is true and such that the first and third terms of (A.12) are small.

Then, we will extend the result by recursion to any $t \in [0, T]$.

**Step 4.2: Construction of auxiliary processes $\hat{P}^i, \hat{A}^j$.** For $n, m$ given, and $t \ge 0$, let us define,

$$ \hat{\mathbf{P}}_t^{n,m}(dx) = \sum_{i=1}^{n} \hat{P}_t^i \delta_{x^i}, \quad \hat{\mathbf{A}}_t^{n,m}(dx) = \sum_{j=1}^{m} \hat{A}_t^j \delta_{y^j}, \qquad (A.13) $$

where for $i \in [1, n]$ and $j \in [1, m]$, the initial state of the auxiliary processes is constructed using a discretization of the limiting initial condition $(a_0(y)dy, p_0(x)dx)$:

$$ \hat{P}_0^i = n \int_{x^i}^{x^{i+1}} p_0(x) dx \quad \text{and} \quad \hat{A}_0^j = m \int_{x^j}^{x^{j+1}} a_0(y) dy \quad \text{with } x^{n+1} = y^{m+1} = 1. \qquad (A.14) $$

Then for all $t \ge 0$, we construct the $\hat{P}^i$'s and $\hat{A}^j$'s as the unique solutions of the following system

$$ \begin{aligned} \frac{d\hat{P}_t^i}{dt} &= \hat{P}_t^i \left( g^P(c\phi \star \hat{\mathbf{A}}_0^{n,m}(x^i)) - \frac{1}{n} \sum_{\ell=1}^n k(x^i, x^\ell) \hat{P}_t^\ell \right) \\ \frac{d\hat{A}_t^j}{dt} &= \hat{A}_t^j \left( g^A(c\phi \star \hat{\mathbf{P}}_0^{n,m}(y^j)) - \frac{1}{m} \sum_{\ell=1}^m h(y^j, y^\ell) \hat{A}_t^\ell \right). \end{aligned} $$

Notice that the processes $(\hat{P}_t^i, t \ge 0)$ and $(\hat{A}_t^j, t \ge 0)$ depend on $c, \phi$ and on the sequences $(x^i)_{i \le n}$ and $(y^j)_{j \le m}$, but not on the $G_{ij}$'s and the $c_{ij}^{n,m}$'s since their initial condition does not. Moreover, proceeding as in Step 1 (to obtain (A.5)), we can prove that for all $t \in [0, T]$,

$$ \sup_{n,m} (\langle \hat{\mathbf{P}}_t^{n,m}, 1 \rangle + \langle \hat{\mathbf{A}}_t^{n,m}, 1 \rangle) \le C(T), \quad \text{a.s.,} \qquad (A.15) $$

hence the second term of (A.12) is smaller than $\varepsilon/3$, as soon as $n, m \ge N_0, M_0$.

**Step 4.3: the third term of (A.12) is small.** Notice that, under the assumptions on $c$ and $\phi$, there exists $C_0 > 0$ such that for all measures $\nu$ and $\mu$, and for all $x \in [0, 1]$

$$ |c\phi \star (\nu - \mu)(x)| \le C_0 W_1(\nu, \mu). $$

Thus, to show this step, we will prove that $W_1(\hat{\mathbf{A}}_t^{n,m}, \tilde{\mathbf{A}}_t^{n,m})$ is small. To this aim, let $f$ be a Lipschitz function such that $\|f\|_\infty \le 1$ and such that its Lipschitz constant $\text{Lip}(f)$ is also bounded by 1, $\text{Lip}(f) \le 1$. According to the equations satisfied by all the processes, for any $t \ge 0$,

$$ \begin{aligned} & \left| \langle \hat{\mathbf{A}}_t^{n,m} - \tilde{\mathbf{A}}_t^{n,m}, f \rangle - \langle \hat{\mathbf{A}}_0^{n,m} - \tilde{\mathbf{A}}_0^{n,m}, f \rangle \right| \\ & \le \int_0^t \left| \frac{1}{m} \sum_{j=1}^m f(y^j) \left[ \tilde{A}_s^j \left( g^A(c\phi \star \hat{\mathbf{P}}_0^{n,m}(y^j)) - g^A\left(\frac{1}{n}\sum_{i=1}^n G_{ij} m c_{ij}^{n,m} \tilde{P}_s^i\right) \right) \right. \right. \\ & \qquad \left. \left. + g^A(c\phi \star \hat{\mathbf{P}}_0^{n,m}(y^j)) (\tilde{A}_s^j - \tilde{A}_s^j) - \frac{1}{m} \sum_{\ell=1}^m h(y^j, y^\ell) (\hat{A}_s^\ell \hat{A}_s^j - \tilde{A}_s^\ell \tilde{A}_s^j) \right] \right| ds \\ & \le \int_0^t \left| C(T) L^A \left( \max_{1 \le j \le m} \left| \langle \hat{\mathbf{P}}_0^{n,m} - \tilde{\mathbf{P}}_0^{n,m}, c\phi(\cdot, y^j) \rangle \right| + K_c \max_{1 \le i \le n} |\tilde{P}_s^i - \tilde{P}_0^i| \right) + \left| \langle \hat{\mathbf{A}}_s^{n,m} - \tilde{\mathbf{A}}_s^{n,m}, \bar{f} \rangle \right| \right. \\ & \qquad \left. + \frac{C(T)}{m} \sum_{j=1}^m \left| \langle \hat{\mathbf{A}}_s^{n,m} - \tilde{\mathbf{A}}_s^{n,m}, h(y^j, \cdot) \rangle \right| + \frac{C(T)}{m} \sum_{\ell=1}^m \left| \langle \tilde{\mathbf{A}}_s^{n,m} - \hat{\mathbf{A}}_s^{n,m}, f h(\cdot, y^\ell) \rangle \right| \right| ds, \end{aligned} $$
where $K_c$ corresponds to the bound defined in Assumption 1.3 and $\bar{f}(y) = f(y)g^A(c\phi \star \hat{\mathbf{P}}_0^{n,m}(y))$ is a Lipschitz continuous and bounded function such that Lip($\bar{f}$) is bounded independently from $n, m$, according to the assumptions on $c$ and $\phi$. Using also that $h$ is a Lipschitz function, we finally deduce that there exists $C_1 > 0$ such that

$$ \left|\langle \hat{\mathbf{A}}_t^{n,m} - \tilde{\mathbf{A}}_t^{n,m}, f \rangle\right| \le \mathcal{W}_1(\hat{\mathbf{A}}_0^{n,m}, \tilde{\mathbf{A}}_0^{n,m}) + \int_0^t C_1 \left( \mathcal{W}_1(\hat{\mathbf{P}}_0^{n,m}, \tilde{\mathbf{P}}_0^{n,m}) + \max_{1 \le i \le n, s \le t} |\tilde{P}_s^i - \tilde{P}_0^i| + \mathcal{W}_1(\hat{\mathbf{A}}_s^{n,m}, \tilde{\mathbf{A}}_s^{n,m}) \right) ds. $$

Taking the supremum in the left hand side w.r.t. $f$, and using Gronwall's lemma, we deduce that for all $t \ge 0$

$$ \mathcal{W}_1(\hat{\mathbf{A}}_t^{n,m}, \tilde{\mathbf{A}}_t^{n,m}) \le \left( \mathcal{W}_1(\hat{\mathbf{A}}_0^{n,m}, \tilde{\mathbf{A}}_0^{n,m}) + C_1 t \mathcal{W}_1(\hat{\mathbf{P}}_0^{n,m}, \tilde{\mathbf{P}}_0^{n,m}) + C_1 t \max_{1 \le i \le n, s \le t} |\tilde{P}_s^i - \tilde{P}_0^i| \right) e^{C_1 t}. \quad (A.16) $$

On the one hand, for the third term in the r.h.s., we have for any $s \in [0, T]$,

$$ \begin{aligned} |\tilde{P}_s^i - \tilde{P}_0^i| &= \left| \int_0^s \left( g^P \left( \frac{1}{m} \sum_{j=1}^m G_{ijm} c_{ij}^{n,m} \tilde{A}_u^j \right) - \frac{1}{n} \sum_{\ell=1}^n k(x^i, x^\ell) \tilde{P}_u^\ell \right) \tilde{P}_u^i du \right| \\ &\le s C(T) (|g^P(0)| + L^P C(T) K_c + \|k\|_\infty C(T)) = sC_2, \end{aligned} $$

where we set $C_2 := C(T)(|g^P(0)| + L^P C(T) K_c + \|k\|_\infty C(T))$, which is a constant independent from $n, m$ and $s$. On the other hand, $\tilde{\mathbf{P}}_0^{n,m}$ and $\tilde{\mathbf{A}}_0^{n,m}$ converge weakly towards $p_0(x)dx$ and $a_0(y)dy$ in probability. Moreover the sequence of r.v. $\mathcal{W}_1(\hat{\mathbf{A}}_0^{n,m}, a_0(y)dy)$ and $\mathcal{W}_1(\hat{\mathbf{P}}_0^{n,m}, p_0(x)dx)$ are bounded in $L^2$, and using the dominated convergence theorem, we deduce the convergences of the expectations to 0. In addition with the definition of $\hat{\mathbf{P}}_0^{n,m}$ and $\hat{\mathbf{A}}_0^{n,m}$ (see (A.14)), we can find $N_1, M_1 \ge N_0, M_0$ such that for all $n, m \ge N_1, M_1$,

$$ \mathbb{E}[\mathcal{W}_1(\hat{\mathbf{A}}_0^{n,m}, \tilde{\mathbf{A}}_0^{n,m})] \le \frac{\varepsilon}{12C_0} \quad \text{and} \quad \mathbb{E}[\mathcal{W}_1(\hat{\mathbf{P}}_0^{n,m}, \tilde{\mathbf{P}}_0^{n,m})] \le \frac{\varepsilon}{12C_0}. \quad (A.17) $$

In addition with (A.16), we deduce that there exists $\delta_1$ such that for all $t \le \delta_1$, all $x \in [0, 1]$ and all $n, m \ge N_1, M_1$,

$$ \begin{aligned} \mathbb{E} \left| c\phi \star (\hat{\mathbf{A}}_t^{n,m} - \tilde{\mathbf{A}}_t^{n,m})(x) \right| &\le C_0 \mathbb{E} \left[ \mathcal{W}_1(\hat{\mathbf{A}}_t^{n,m}, \tilde{\mathbf{A}}_t^{n,m}) \right] \\ &\le \left( \frac{\varepsilon}{12} + C_1 \delta_1 \frac{\varepsilon}{12} + C_0 C_1 C_2 \delta_1^2 \right) e^{C_1 \delta_1} \le \frac{\varepsilon}{3}. \end{aligned} \quad (A.18) $$

Using the same computations, we can find $\delta_0 \le \delta_1$ such that for all $\delta \le \delta_0$ and all $n, m \ge N_1, M_1$, the same result holds for the processes $\hat{\mathbf{P}}_0^{n,m}$ and $\tilde{\mathbf{P}}_0^{n,m}$.

**Step 4.4: the first term of (A.12) is small.** We will use similar computations as in Step 4.3 to reach this step. To this aim, we define for all $t \ge 0$

$$ \Sigma^{A,n,m}(t) = \sup \left\{ \frac{1}{m} \sum_{j=1}^{m} G_{ijm} c_{ij}^{n,m} f(y^j) (\tilde{A}_t^j - \hat{A}_t^j) \mid f : [0, 1] \to [0, 1] \text{ Lipschitz, } \text{Lip}(f) \le 1 \right\}, $$
and a similar definition for $\Sigma^{P,n,m}(t)$. Let $f$ be a Lipschitz function, bounded by 1. Then

$$ \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} f(y^j) (\tilde{A}_t^j - \hat{A}_t^j) \right| - \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} f(y^j) (\tilde{A}_0^j - \hat{A}_0^j) \right| \\ \le \int_0^t \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} f(y^j) g^A (c\phi \star \tilde{\mathbf{P}}_0^{n,m}(y^j)) (\tilde{A}_s^j - \hat{A}_s^j) \right. \\ + \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} f(y^j) \tilde{A}_s^j \left( g^A \left( \frac{1}{n} \sum_{i=1}^{n} G_{ij} n c_{ij}^{n,m} \tilde{P}_s^i \right) - g^A (c\phi \star \tilde{\mathbf{P}}_0^{n,m}(y^j)) \right) \\ + \frac{1}{m} \sum_{\ell=1}^{m} \hat{A}_s^\ell \left( \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} f(y^j) h(y^j, y^\ell) (\tilde{A}_s^j - \tilde{A}_s^\ell) \right) \\ \left. + \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} f(y^j) \tilde{A}_s^j \frac{1}{m} \sum_{\ell=1}^{m} h(y^j, y^\ell) (\tilde{A}_s^\ell - \tilde{A}_s^\ell) \right| ds. $$

Following the same argument as in the previous step, using Gronwall's lemma, we deduce that

$$ \Sigma^{A,n,m}(t) \leq (\Sigma^{A,n,m}(0) + C_1 t \Sigma^{P,n,m}(0) + C_1 t C_2) e^{C_1 t}, \quad (A.19) $$

where $C_1$ and $C_2$ have been defined in the previous step.

It remains to prove that $\mathbb{E}[\Sigma^{A,n,m}(0)]$ and $\mathbb{E}[\Sigma^{P,n,m}(0)]$ converge to 0 with $n$ and $m$. To this aim, we follow the exact same computations as in Step 4.1 to deduce that for any Lipschitz function $f$ with $\|f\|_{\infty} \le 1$,

$$ \mathbb{E} \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} f(y^j) A_0^j - (c\phi f) \star \mathbf{A}_0^{n,m}(x^i) \right| \le \frac{C(T)}{\sqrt{m}} \left( V_{\max}^{1/2} + \|mc^{n,m}\|_{\infty} \|\phi(1-\phi)\|_{\infty}^{1/2} \right) \\ + C(T) \|\phi\|_{\infty} \|mc^{n,m} - c\|_{\infty}, $$

with $A_0^j = \tilde{A}_0^j$ or $\hat{A}_0^j$. In addition with the weak convergence of $\tilde{\mathbf{A}}_0^{n,m}$ and $\hat{\mathbf{A}}_0^{n,m}$ in probability towards the same limits, and with similar computations for $P$, we deduce that there exists $N_2, M_2 \ge N_1, M_1$ such that

$$ \mathbb{E}[\Sigma^{A,n,m}(0)] \le \frac{\varepsilon}{12C_0} \quad \text{and} \quad \mathbb{E}[\Sigma^{P,n,m}(0)] \le \frac{\varepsilon}{12C_0}. $$

Using also (A.19), and remembering the definition of $\delta_0$ at the end of Step 4.2, we deduce that for all $n, m \ge N_2, M_2$ and all $t \le \delta_0$,

$$ \mathbb{E} \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} (\tilde{A}_t^j - \hat{A}_t^j) \right| \le C_0 \mathbb{E} [\Sigma^{A,n,m}(t)] \le \frac{\varepsilon}{3}. \quad (A.20) $$

Moreover, the same result holds for the processes $\tilde{\mathbf{P}}^{n,m}$ and $\hat{\mathbf{P}}^{n,m}$.

### Step 4.5: conclusion for $t \le \delta_0$.

Combining Equations (A.12), (A.20), (A.18) and the end of Step 4.2, we conclude that the limit in (A.10) is true uniformly for $t \in [0, \delta_0]$. More precisely, for all $n, m \ge N_2, M_2$ and all $t \le \delta_0$,

$$ \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} \tilde{A}_t^j - (c\phi) \star \tilde{\mathbf{A}}_t^{n,m}(x^i) \right| \right] \le \varepsilon. $$
Finally, (A.8) is valid for all $t \le \delta_0$ and any limiting value satisfies the equation (3.4) for $t \in [0, \delta_0]$. According to the uniqueness of the solution to this equation (see Step 5 below), we deduce in particular that

$$ \lim_{n,m \to +\infty} (W_1(\tilde{\mathbf{P}}_{\delta_0}^{n,m}, \bar{\mathbf{P}}_{\delta_0}) + W_1(\tilde{\mathbf{A}}_{\delta_0}^{n,m}, \bar{\mathbf{A}}_{\delta_0})) = 0, \quad \text{in probability,} \quad (A.21) $$

with $(\bar{\mathbf{P}}_t, \bar{\mathbf{A}}_t)_{t \in [0, \delta_0]}$ solution to (3.4). According to Step 3, $\bar{\mathbf{P}}_t$ and $\bar{\mathbf{A}}_t$ admit densities $\bar{p}_t$ and $\bar{a}_t$ with respect to the Lebesgue measure.

**Step 4.6: recursion.** We proceed now by recursion, by defining new auxiliary processes, for $n, m$ given, $t \ge \delta_0$,

$$ \hat{\mathbf{P}}_t^{n,m}(dx) = \sum_{i=1}^{n} \hat{P}^i \delta_{x^i}, \quad \hat{\mathbf{A}}_t^{n,m}(dx) = \sum_{j=1}^{m} \hat{A}^j \delta_{y^j}, $$

with for $i \in [1, n]$ and $j \in [1, m]$, $\hat{P}_{\delta_0}^i$ and $\hat{A}_{\delta_0}^j$ are defined in a similar way as the beginning of Step 4.2 using $\bar{p}_{\delta_0}(x)dx$ and $\bar{a}_{\delta_0}(y)dy$ instead of $p_0(x)dx$ and $a_0(y)dy$:

$$ \hat{P}_{\delta_0}^i = n \int_{x^i}^{x^{i+1}} p_{\delta_0}(x)dx \quad \text{and} \quad \hat{A}_{\delta_0}^j = m \int_{x^j}^{x^{j+1}} a_{\delta_0}(y)dy \quad \text{with } x^{n+1} = y^{m+1} = 1. \quad (A.22) $$

And for all $t \ge \delta_0$, we set

$$ \begin{aligned} \frac{d\hat{P}_t^i}{dt} &= \hat{P}_t^i \left( g^P(c\phi \star \hat{\mathbf{A}}_{\delta_0}(x^i)) - \frac{1}{n} \sum_{\ell=1}^n k(x^i, x^\ell) \hat{P}_t^\ell \right) \\ \frac{d\hat{A}_t^j}{dt} &= \hat{A}_t^j \left( g^A(c\phi \star \hat{\mathbf{P}}_{\delta_0}(y^j)) - \frac{1}{m} \sum_{\ell=1}^m h(y^j, y^\ell) \hat{A}_t^\ell \right). \end{aligned} $$

Notice again that the processes $(\hat{P}_t^i)$ and $(\hat{A}_t^j)$ do depend on $c$, $\phi$ and on the sequences $(x^i)_{i \le n}$ and $(y^j)_{j \le m}$, but not on the $G_{ij}$'s and the $c_{ij}^{n,m}$'s. Thus, we can proceed as in Step 4.2, 4.3, 4.4 and 4.5 and find $N_4, M_4 \ge N_2, M_2$ such that for all $n, m \ge N_4, M_4$ and all $t \in [\delta_0, 2\delta_0]$,

$$ \mathbb{E} \left[ \left| \frac{1}{m} \sum_{j=1}^{m} G_{ij} m c_{i,j}^{n,m} \tilde{A}_t^j - (c\phi) \star \tilde{\mathbf{A}}_t^{n,m}(x^i) \right| \right] \le \varepsilon. \quad (A.23) $$

In other words, we can conclude similarly as in Step 4.5 for all $t \in [0, 2\delta_0]$. With a direct recursive argument, we finally be able to find $N_{2\kappa}, M_{2\kappa}$, with $\kappa = [T/\delta_0]$, such that for all $n, m \ge N_{2\kappa}, M_{2\kappa}$ and all $t \in [0, T]$, (A.23) holds.

This finally proves (A.10) and thus (A.8), which finishes the proof of Step 4.

**Step 5: Uniqueness of the solution of (3.4).** It remains to prove the uniqueness of the solution to (3.4) 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.4) and with identical initial conditions. Let us denote by $\rho$ the Radon metric between two measures $\nu$ and $\mu$ on $M_F([0, 1])^2$:

$$ \rho(\nu, \mu) = \sup \left\{ \int_0^1 f(x) d(\nu - \mu)(x) \middle| 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