Mean field approximation of a heterogeneous population of plants in competition

TL;DR

This paper develops a mean-field approximation for modeling large, heterogeneous plant populations in competition, simplifying complex interactions into a manageable mathematical framework with potential applications in ecological inference.

Contribution

It introduces a mean-field approach to analyze large heterogeneous populations with pairwise interactions, enabling simplified modeling and inference in ecological systems.

Findings

Mean-field approximation effectively models large plant populations.

Simulation using semi-Lagrangian scheme and Gaussian process regression.

Asymptotic independence of individuals in large populations.

Abstract

The processes of interplant competition within a field are still poorly understood. However, they explain a large part of the heterogeneity in a field and may have longer-term consequences, especially in mixed stands. Modeling can help to better understand these phenomena but requires simulating the interactions between different individuals. In the case of large populations, assessing the parameters of a heterogeneous population model from experimental data is intractable computationally. This paper investigates the mean-field approximation of large dynamical systems with random initial conditions and individual parameters, and with interaction being represented by pairwise potentials between individuals. Under this approximation, each individual is in interaction with an infinitely-crowded population, summarized by a probability measure, the mean-field limit distribution, being itself the weak solution of a non-linear hyperbolic partial differential equation. In particular, the phenomenon of chaos propagation implies that the individuals are independent asymptotically when the size of the population tends towards infinity. This result provides perspectives for a possible simplification of the inference problem. The simulation of the mean-field distribution, consisting in a semi-Lagrangian scheme with an interpolation step using Gaussian process regression, is illustrated for a heterogeneous population model representing plants in competition for light.

Tables

Table 1. Table 1: Configuration of the parameters for the simulation of the system 9

$L$	1 m	$S_M$	1 m
$S_m$	0.8 m	${\gamma }_M$	1 day-1
${\gamma }_m$	0.1 day-1	${\sigma }_S$	${10}^{-2}$ m
${\sigma }_{\gamma }$	${10}^{-2}$ day-1	$s^0$	0.3 m
$s_m$	${5.10}^{-2}$ m	$R_M$	$\log \left(S_M/s_m\right)$
${\sigma }_x$	$L$	${\sigma }_r$	$\log \left(0.1/s_m\right)$
$\Delta t$		0.1 day

Equations

$$\forall i\in\llbracket 1;N\rrbracket,~{}\left\{\begin{array}[]{l}(X_{i}^{0},\theta_{i})\sim\mu_{0}\mathrm{~{}a~{}probability~{}measure}\\ X_{i}(0)=X_{i}^{0}\\ \displaystyle\frac{\mathrm{d}X_{i}(t)}{\mathrm{d}t}=F\left(X_{i}(t),\theta_{i};(X_{j}(t),\theta_{j})_{1\leq j\leq N}\right)\end{array}\right.$$

$$\forall i\in [[1;N]],F\left(X_i(t),{\theta }_i;(X_j(t),{\theta }_j{)}_{1\le j\le N}\right)=\frac{1}{N-1}\sum _{1\le j\le N,j\ne i}g(X_i(t),{\theta }_i,X_j(t),{\theta }_j)$$

$$\forall z=(X,\theta )\in \mathcal{Z},|z|=\sum _{i=1}^{d_{\mathcal{X}}+d_{\Theta }}\frac{|z_i|}{|z_i^{\ast }|}$$

$$(X_i^0,{\theta }_i{)}_{1\le i\le N}\sim {\mu }_0^{\otimes N}$$

$$\displaystyle\forall i\in\llbracket 1;N\rrbracket,~{}\left\{\begin{array}[]{l}X_{i}(0)=X_{i}^{0}\\ \forall t\in\mathbb{R}_{+},~{}\displaystyle\frac{\mathrm{d}X_{i}(t)}{\mathrm{d}t}=\frac{1}{N-1}\sum_{j\neq i}g(X_{i}(t),\theta_{i},X_{j}(t),\theta_{j})\end{array}\right.$$

$$|g(X_1,\theta ,X_1^{\prime },{\theta }^{\prime })-g(X_2,\theta ,X_2^{\prime },{\theta }^{\prime })|\le K_2(1+|X_1^{\prime }|+|X_2^{\prime }|)(|X_1-X_2|+|X_1^{\prime }-X_2^{\prime }|)$$

$$\forall (X_1,{\theta }_1,X_2,{\theta }_2)\in {\mathcal{Z}}^2,$$

$$\left|\frac{\partial g}{\partial X}(X_1,{\theta }_1,X_1,{\theta }_2)\right|=\underset{X\in \mathcal{X},|X|=1}{\sup }\left|\frac{\partial g}{\partial X}(X_1,{\theta }_1,X_1,{\theta }_2).X\right|\le K_3(1+|X_2|)$$

$$|g(X,{\theta }_1,X^{\prime },{\theta }_1^{\prime })-g(X,{\theta }_2,X^{\prime },{\theta }_2^{\prime })|\le K_4(1+|X|+|X^{\prime }|)(|{\theta }_1-{\theta }_2|+|{\theta }_1^{\prime }-{\theta }_2^{\prime }|)$$

$$\left\{\begin{array}[]{l}s(0)=s^{0}\\ \displaystyle\forall t\in\mathbb{R}_{+},~{}\frac{\mathrm{d}s(t)}{\mathrm{d}t}=\gamma s(t)\log\left(\frac{S}{s(t)}\right)\end{array}\right.\Rightarrow s(t)=S\exp\left(-e^{-\gamma t}\log\left(\frac{S}{s^{0}}\right)\right)$$

$$\forall i\in [[1;N]],$$

$$\displaystyle\left\{\begin{array}[]{l}s_{i}(0)=s_{i}^{0}\\ \displaystyle\forall t\in\mathbb{R}_{+},~{}\frac{\mathrm{d}s_{i}(t)}{\mathrm{d}t}=\gamma_{i}s_{i}(t)\left(\log\left(\frac{S_{i}}{s_{m}}\right)\left(1-\frac{1}{N-1}\sum_{j\neq i}C(s_{i}(t),s_{j}(t),|\vec{x}_{i}-\vec{x}_{j}|)\right)\right.\\ \displaystyle\left.-\log\left(\frac{s_{i}(t)}{s_{m}}\right)\right)\end{array}\right.$$

$$\forall i\in [[1;N]],\forall t\in {\mathbb{R}}_{+},{\gamma }_i{s}_i(t)\log \left(\frac{s_m}{s_i(t)}\right)\le \frac{\mathrm{d}{s}_i(t)}{\mathrm{d}t}\le {\gamma }_i{s}_i(t)\log \left(\frac{S_i}{s_i(t)}\right)$$

$$x\sim \mathcal{U}([0;L])y\sim \mathcal{U}([0;L])\text{ and }x,y\text{ are independent}$$

$$S|x\sim \mathcal{U}([S_1(x);S_2(x)])$$

$$\text{with }{S}_1(x)=S_m+\frac{x}{L}(S_M-{\sigma }_S-S_m),S_2(x)=S_m+{\sigma }_S+\frac{x}{L}(S_M-{\sigma }_S-S_m)$$

$$\text{with }{S}_m\text{ and }{\sigma }_S\text{ such that }{S}_m>0,{\sigma }_S>0,S_m+{\sigma }_S<S_M$$

$$\gamma |y\sim \mathcal{U}([{\gamma }_1(y);{\gamma }_2(y)])$$

$$\text{with }{\gamma }_1(y)={\gamma }_m+\frac{y}{L}({\gamma }_M-{\sigma }_{\gamma }-{\gamma }_m),{\gamma }_2(y)={\gamma }_m+{\sigma }_{\gamma }+\frac{y}{L}({\gamma }_M-{\sigma }_{\gamma }-{\gamma }_m)$$

$$s^0\sim {\delta }_{s^0}\text{(the initial size of the plants is a constant over the population)}$$

$$\forall (x,y,S,\gamma )\in {\mathbb{R}}^4,p_0^{\theta }(x,y,S,\gamma )=\frac{\mathbb{I}\{0\le x,y\le L\}\mathbb{I}\{S_1(x)\le S\le S_2(x)\}\mathbb{I}\{{\gamma }_1(y)\le \gamma \le {\gamma }_2(y)\}}{L^2{\sigma }_S{\sigma }_{\gamma }}$$

$$\forall z=(r,x,y,S,\gamma )\in \mathcal{Z},|z|=|r|+\frac{|x|+|y|}{L}+\frac{|S|}{s_m}+\frac{|\gamma |}{{\gamma }_m}$$

$$K_1={\gamma }_M\max (1,R_M)$$

$$K_3={\gamma }_M\max \left(1,\frac{1}{2{\sigma }_r}\right)$$

$$t\in {\mathbb{R}}_{+}\mapsto \mu [t,Z_N^0]=\frac{1}{N}\sum _{i=1}^N{\delta }_{z_i(t,Z_N^0)}$$

$${\mathcal{C}}_0^1({\mathbb{R}}_{+}\times \mathcal{Z}\to \mathbb{R})=\{\phi :{\mathbb{R}}_{+}\times \mathcal{Z}\to \mathbb{R}\text{ continuously differentiable}|$$

$$\underset{|z|\to +\infty }{\lim }|\phi (t,z)|+\left|\frac{\partial \phi }{\partial z}(t,z)\right|+\left|\frac{\partial \phi }{\partial t}(t,z)\right|=0\}$$

$${\mathbb{E}}_{\mu [t,Z_N^0]}(\phi (t,z))={\int }_{\mathcal{Z}}\phi (t,z)\mu [t,Z_N^0](\mathrm{d}z)=\frac{1}{N}\sum _{i=1}^N\phi (t,z_i(t,Z_N^0))$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\mathcal{Z}}\phi (t,z)\mu [t,Z_N^0](\mathrm{d}z)=\frac{1}{N}\sum _{i=1}^N\frac{\partial \phi }{\partial t}(t,z_i(t,Z_N^0))+\frac{\partial \phi }{\partial X}(t,z_i(t,Z_N^0){)}^{\text{T}}\frac{\mathrm{d}{X}_i}{\mathrm{d}t}(t,Z_N^0)$$

$$={\int }_{\mathcal{Z}}\frac{\partial \phi }{\partial t}(t,z)\mu [t,Z_N^0](\mathrm{d}z)+\frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}\frac{\partial \phi }{\partial X}(t,z_i(t,Z_N^0){)}^{\text{T}}g(z_i(t,Z_N^0),z_j(t,Z_N^0))$$

$$\frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}\frac{\partial \phi }{\partial X}(t,z_i(t,Z_N^0){)}^{\text{T}}g(z_i(t,Z_N^0),z_j(t,Z_N^0))=$$

$${\int }_{\mathcal{Z}}\frac{\partial \phi }{\partial X}(t,z{)}^{\text{T}}\left(\frac{N}{N-1}{\int }_{\mathcal{Z}}g(z,z^{\prime })\mu [t,Z_N^0](\mathrm{d}{z}^{\prime })-\frac{1}{N-1}g(z,z)\right)\mu [t,Z_N^0](\mathrm{d}z)$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\mathcal{Z}}\phi (t,z)\mu [t,Z_N^0](\mathrm{d}z)=$$

$${\int }_{\mathcal{Z}}\left(\frac{\partial \phi }{\partial t}(t,z)+\frac{\partial \phi }{\partial X}(t,z{)}^{\text{T}}\left(\frac{N}{N-1}{\int }_{\mathcal{Z}}g(z,z^{\prime })\mu [t,Z_N^0](\mathrm{d}{z}^{\prime })-\frac{1}{N-1}g(z,z)\right)\right)\mu [t,Z_N^0](\mathrm{d}z)$$

$$\forall z\in \mathcal{Z},{\mathcal{G}}_N(\mu [t,Z_N^0],z)=\frac{N}{N-1}{\int }_{\mathcal{Z}}g(z,z^{\prime })\mu [t,Z_N^0](\mathrm{d}{z}^{\prime })-\frac{1}{N-1}g(z,z)$$

$$\frac{\partial f}{\partial t}(t,z)+{\mathrm{div}}_X\left(f(t,z)\left(\frac{N}{N-1}{\int }_{\mathcal{Z}}g(z,z^{\prime })f(t,z^{\prime }){\lambda }^{\otimes d_z}(\mathrm{d}{z}^{\prime })-\frac{g(z,z)}{N-1}\right)\right)=0$$

$$\text{or }\frac{\partial f}{\partial t}(t,z)+{\mathrm{div}}_X\left(f(t,z){\mathcal{G}}_N\left(f(t,\cdot ){\lambda }^{\otimes d_z},z\right)\right)=0$$

$$\Pi ({\mu }_1,{\mu }_2)=\left\{\pi \in {\mathcal{P}}_1({\mathcal{Z}}^2)|{\mu }_1={\int }_{\mathcal{Z}}\pi (.,\mathrm{d}{z}_2),{\mu }_2={\int }_{\mathcal{Z}}\pi (\mathrm{d}{z}_1,.)\right\}$$

Full text

1 Introduction

The interest for modelling heterogeneous populations of plants is on the rise, especially due to the development of the practice of mixed cropping (Malézieux et al. [2009]). Mixing different varieties or different species (Tang et al. [2018]) may have various advantages, such as nitrogen transfer from one species to another, resistance of the population to disease and pests (Gurr et al. [2003]), or enhanced production quality (Gooding et al. [2007] for wheat). However, up to our knowledge, very few models are made to understand the emerging properties of such mixture and to design optimal crops (cf. Gaudio et al. [2019] for a review). A convenient framework for modelling heterogeneous populations is hierarchical modelling, also known as mixed-effects modelling (Schneider et al. [2006], Lv et al. [2008], Baey et al. [2016]). A classical formulation of hierarchical model of a population dynamics can be represented as a dynamical system, whose initial conditions and parameters are independent and identically distributed random variables.

[TABLE]

$N$ is the number of individuals in the population. Each individual is indexed by integer $i\in\llbracket 1;N\rrbracket$ and is described by a state variable $X_{i}$ and individual parameter $\theta_{i}$. The state variable $X_{i}$ represents time-varying features of the plant, e.g. the size of its aerial part, the total leaf area, etc. The individual parameter $\theta_{i}$ represents intrinsic characteristics of individual $i$, that are assumed to be constant throughout the considered time period, and that have influence on the population dynamics. $F$ is a function modelling the influence of the whole population, consisting in the collection $(X_{j},\theta_{j})_{1\leq j\leq N}$, on the individual development of each plant. A specific form of $F$ is going to be studied in the present article (see equation (2)). The heterogeneity of the population is represented by the probability measure $\mu_{0}$, that distributes the initial state variable and individual parameter to each individual at the population level. If the marginal distribution of variable $\theta$, $\mu_{0}^{\theta}$, is not reduced to a Dirac distribution, or equivalently if $\theta$ is not constant over the population, then the population is said to be heterogeneous, as it gathers individuals with different characteristics. The case of homogeneous population has been investigated for example in Cournède et al. [2007], Sievänen et al. [2008] focusing on the competition between plants.

The problem of statistical inference on such population model consists in identifying distribution $\mu_{0}$ and function $F$ from collected observation data. In the case where the plants do not interact with each other, various forms of Expectation - Maximization (EM) algorithm, introduced by Dempster et al. [1977], can be applied to estimate the parameters (Baey et al. [2016]),Baey et al. [2018], or direct Bayesian inference (Viaud [2018], chapter 4). Most common forms of EM algorithm and of direct Bayesian inference require a random exploration of the unknown parameter space using Metropolis-Hasting (MH) algorithm, or Metropolis-Hasting within Gibbs (MHWG) algorithmn, which are not suited for the exploration of high-dimensional space in terms of convergence time (Katafygiotis and Zuev [2008]). Nevertheless, EM algorithm or MHWG algorithm remain efficient tools for parameter estimation in a population model without interaction.

The relative effectiveness of these algorithms is challenged when taking into account interactions within the population model. The correlations between individuals hinder the distribution of the computation and the search space where MH algorithm is applied is of too high dimension, proportional to the number $N$ of individuals (cf. the computational issue encountered in Schneider et al. [2006]). The aim of this research is to suggest other methods more suited to this problem.

A possible research direction is given by variational Bayesian approximation (cf. Bishop [2006], chapter 10, for an introduction). This method consists in projecting the joint distribution of the random variables $(X_{i}(t),\theta_{i})_{1\leq i\leq N}$, which is non-factorized due to individuals interaction, onto a tensor product of parametric distributions. For specific expression of the function $F$, such as the interaction function used in the Cucker-Smale model (Cucker and Smale [2007], Carrillo et al. [2010]), any subset of the the population has a joint distribution asymptotically factorized as $N\rightarrow+\infty$, a phenomenon referred as chaos propagation in the literature (Bolley et al. [2011]). Qualitatively, when the population is infinite, the states at time $t$ and parameters of the individuals behave as if they were independent random variables distributed according to a single probability measure $\mu[t]$, that is called the mean field limit (MFL) distribution.

The question on how to integrate the MFL distribution $\mu[t]$ into the process of statistical inference is beyond the scope of this article. The first step is to check for which kind of heterogeneous population models we can obtain theoretical existence and uniqueness of MFL distribution, along with asymptotic factorization property. We have considered plant population models for which the interaction function $F$ can be decomposed as a sum of elementary interaction functions over the whole population.

[TABLE]

Such formulation is used in Schneider et al. [2006], Lv et al. [2008], Nakagawa et al. [2015]. These models focus on the competition for light within the plant population. Schneider et al. [2006] suggests various models coupling plant development and competition. Amongst the models being smooth enough, we have chosen the one with the most statistical relevance. This model is described in section 2. Equation (2) is quite close in its formulation to particle systems studied in kinetic equation theory (Carrillo et al. [2010]). The normalization by the size of the population is important for the study of the asymptotic behavior of the population as $N$ tends towards infinity. The derivative of an individual state remains of the same order of magnitude when the size of the population changes. In our case, this normalization is part of the model expression, but for other systems (such as the ones studied in statistical physics), this normalization can be interpreted as a change of time scale (Golse [2013a]). Other normalization can be considered in some flocking model, like Vicsek model (Vicsek et al. [1995], Degond [2018]) where the velocity is normalized by the sum of all the velocities in the population. We shall specify in the next section the assumptions to be made on $g$ and on $\mu_{0}$ to derive the MFL distribution. We give also an example of plant competition model from Schneider et al. [2006] to illustrate the theoretical development in section 3 to prove existence and uniqueness of MFL distribution, and finally chaos propagation. As our initial aim is to be able to use the MFL distribution for statistical inference, we present a preliminary work in section 4 to approximate this distribution.

2 Example and assumptions

2.1 Working assumptions and notations

In this subsection, we specify our assumptions on the systems (1). Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space. Let $\mathcal{X}$ be an Euclidean space of dimension $d_{\mathcal{X}}$ and $\Theta$ be a compact subset of an Euclidean space, such that the dimension of $\Theta$ is $d_{\Theta}$. The phase space is denoted by $\mathcal{Z}=\mathcal{X}\times\Theta$ and the set of probability measures defined over $\mathcal{Z}$ is denoted by $\mathcal{P}(\mathcal{Z})$. The set of probability measures $\mathcal{P}(\mathcal{Z})$ is associated to the space of random variables, i.e. the space of functions $f:\Omega\rightarrow\mathcal{Z}$ measurable for the measure $\mathbb{P}$. $\mathcal{Z}$ is often endowed with the Lebesgue measure, denoted by $\lambda^{\otimes d_{z}}$ ($d_{z}=\mathrm{dim}(\mathcal{Z})$). Unless otherwise stated, the metric used on $\mathcal{Z}$ is defined by

[TABLE]

where $z^{*}\in\mathcal{Z}$ is a reference vector with components all non-zero, and $|.|:a\in\mathbb{R}\mapsto|a|$ is the absolute norm over $\mathbb{R}$. The norm $|z|$ is therefore a dimensionless quantity. Similarly, we use the notation $\displaystyle\forall X\in\mathcal{X},~{}|X|=\sum_{i=1}^{d_{\mathcal{X}}}\frac{|X_{i}|}{|X_{i}^{*}|}$ and $\displaystyle\forall\theta\in\Theta,~{}|\theta|=\sum_{i=1}^{d_{\Theta}}\frac{|\theta_{i}|}{|\theta_{i}^{*}|}$. We consider the population model of initial distribution $\mu_{0}\in\mathcal{P}(\mathcal{Z})$ a probability measure over $\mathcal{Z}$ and of interaction function $g:\mathcal{Z}\rightarrow\mathcal{X}$.

[TABLE]

We shall use two notations for the interaction function $g$ : either $g:({X_{1},\theta_{1},X_{2},\theta_{2}})\in\mathcal{Z}^{2}\mapsto g({X_{1},\theta_{1},X_{2},\theta_{2}})\in\mathcal{X}$, either $g:({z_{1},z_{2}})\in\mathcal{Z}^{2}\mapsto g({z_{1},z_{2}})\in\mathcal{X}$. Here are some assumptions on the smoothness of function $g$.

(A1) Assumption 1 :

There exists $K_{1}>0$ such that for all $X_{1},X_{2}\in\mathcal{X}$ and all $\theta_{1},\theta_{2}\in\Theta$ $|g({X_{1},\theta_{1},X_{2},\theta_{2}})|\leq K_{1}(1+|X_{1}|+|X_{2}|)$. This assumption makes possible the existence of global solution over $\mathbb{R}_{+}$.

(A2) Assumption 2 :

There exists $K_{2}>0$ such that for all $X_{1},X_{1}^{\prime},X_{2},X_{2}^{\prime}\in\mathcal{X}$ and $\theta,\theta^{\prime}\in\Theta$ we have

[TABLE]

(A3) Assumption 3

: The transition function $g$ has a partial derivative with respect to the variable $X$, $(X,\theta,X^{\prime},\theta^{\prime})\in\mathcal{Z}^{2}\mapsto\displaystyle\frac{\partial g}{\partial X}(X,\theta,X^{\prime},\theta^{\prime})\in\mathcal{M}_{d_{\mathcal{X}}}(\mathbb{R})$, which is continuous and which is such that there exists $K_{3}>0$

[TABLE]

(A4) Assumption 4

: There exists a constant $K_{4}>0$ such that for all $X,X^{\prime}\in\mathcal{X}$ and $\theta_{1},\theta_{1}^{\prime},\theta_{2},\theta^{\prime}_{2}\in\Theta$111As $\Theta$ is not a vector space, we can have $\theta_{1}-\theta_{2}$ not belonging to $\Theta$. The notation $|\theta_{1}-\theta_{2}|$ has therefore to be understood as the norm of the vector $\theta_{1}-\theta_{2}$ in the Euclidean space containing the compact subset $\Theta$.

[TABLE]

The next subsection gives an example of differential system, where the interaction function $g$ satisfies all four assumptions listed above.

2.2 Example of Schneider model

The article of Schneider et al. [2006] studies a population of plants (Arabidopsis thaliana) in competition for light resources. A dozen of models, more or less empirical, are suggested in this paper to represent the growth of the aerial part of plants subject to the shade of its surroundings, and all these models are compared statistically against experimental data. A similar approach is carried out in Nakagawa et al. [2015] at the scale of a whole forest, observed for several decades. The population was then assumed to be homogeneous, certainly because of the computational issues previously mentioned.

In this model, the soil and water resources are assumed to be in abundance, so that the competition concerned only the light resource. Therefore, only the aerial part222In the case of A. thaliana, $s$ can be the diameter of the rosette (see figure 1) of the plant is represented by the model. A plant is described by the size of its aerial part $s$, its position $\vec{x}=(x,y)$ in the plane, and by two intrinsic factors $\gamma$ and $S$, determining properties of the individual growth. Over the time, only plants’ sizes change. The assumptions of the model are the following :

If the plant grows in isolation, or if the influence of competitors can be neglected, the dynamics of its growth is given by a Gompertz function (Paine et al. [2012]).

[TABLE]

The size of the plant converges towards an equilibrium size $S$ with rate $\gamma$. In a more accurate modelling, this equilibrium size should be a function of the environmental conditions, but they are not taken into account here (the light environment is assumed to be controlled). The initial size of the plant $s^{0}>0$ can be thought as the size of the sprout just after emergence. 2. 2.

If the plant grows in presence of competitors, in a population consisting of $N$ individuals, the equilibrium size $S_{i}$ of the individual $i\in\llbracket 1;N\rrbracket$ is perturbed by a factor representing the negative impact of the competition.

[TABLE]

where $\displaystyle C(s_{i},s_{j},|\vec{x}_{i}-\vec{x}_{j}|)=\frac{\log(s_{j}/s_{m})}{2R_{M}\displaystyle\left(1+\frac{|\vec{x}_{i}-\vec{x}_{j}|^{2}}{\sigma_{x}^{2}}\right)}\left(1+\tanh\left(\frac{1}{\sigma_{r}}\log\left(\frac{s_{j}}{s_{i}}\right)\right)\right)$ with $s_{m},\sigma_{x},\sigma_{r}$ being known positive constants and $R_{M}$ such that $\forall i\in\llbracket 1;N\rrbracket~{}\displaystyle\log\left(\frac{S_{i}}{s_{m}}\right)\leq R_{M}$.

In presence of competition, the available light environment of plant $i$, represented by the term $\displaystyle\log\left(\frac{S_{i}}{s_{m}}\right)$, is reduced by a competition factor $\displaystyle 1-\frac{1}{N-1}\sum_{1\leq j\leq N,j\neq i}C(s_{i}(t),s_{j}(t),|\vec{x}_{i}-\vec{x}_{j}|)$, which is dimensionless and takes values in $[0;1]$. The competition exerted on plant $i$ is all the more important than other plants are

tall in absolute terms, with the factor $\displaystyle\frac{\log(s_{j}/s_{m})}{R_{M}}$ 2. 2.

taller than plant $i$, with the factor $\displaystyle\frac{1}{2}\left(1+\tanh\left(\frac{1}{\sigma_{r}}\log\left(\frac{s_{j}}{s_{i}}\right)\right)\right)$ 3. 3.

close to plant $i$, with the factor $\displaystyle\frac{1}{\displaystyle\left(1+\frac{|\vec{x}_{i}-\vec{x}_{j}|^{2}}{\sigma_{x}^{2}}\right)}$

There exist more realistic and complex models to represent competition for light in plants population. Beyer et al. [2015] describes tree crowns development by a transport equation on foliage density. In this model, light ressource is allocated to the different individuals proportionally to their foliage volume. More mechanistic models can be found in the literature, namely the ones making use of Functional Structural Plant Models (FSPM), where the light environment is directly computed by ray tracing through a 3D reconstruction of the canopy (Cieslak et al. [2008]). Such models of competition are still too complex for the method we describe in this article.

Before going any further, we need to prove that system (9) is well-posed for any initial condition.

Proposition 1**.**

Let us consider the initial conditions $(s_{i}^{0})_{1\leq i\leq N}\in(\mathbb{R}_{+}^{*})^{N}$ and the collection of parameters $(x_{i},y_{i},S_{i},\gamma_{i})\in(\mathbb{R}^{2}\times(\mathbb{R}_{+}^{*})^{2})^{N}$. Then the system (9) has an unique solution $s_{1:N}:t\mapsto(s_{i}(t))_{1\leq i\leq N}$ defined over $\mathbb{R}_{+}$ taking positive values, i.e. verifying $\forall t\in\mathbb{R}_{+},~{}\forall i\in\llbracket 1;N\rrbracket,~{}s_{i}(t)>0$.

The existence of the solution of this system is a classical application of Cauchy-Lipschitz theorem to the system satisfied by the vector $r_{1:N}(t)=\left(\log\left(\frac{s_{i}(t)}{s_{m}}\right)\right)_{1\leq i\leq N}\in\mathbb{R}^{N}$. Details of the proof can be found in appendix 6.1. We need also to check for which conditions the global solution given by proposition 1 is consistent with the biological assumptions of the model. A solution $s_{1:N}:t\in\mathbb{R}_{+}\mapsto s_{1:N}(t)\in\mathbb{R}_{+}^{N}$ is consistent with the assumptions of the model if it meets the following constraints :

• •

The size of each individual must remain below its equilibrium size and above the minimal size $s_{m}$, i.e. for all $i\in\llbracket 1;N\rrbracket,~{}s_{m}<s_{i}(t)\leq S_{i}$.

• •

The competition factor must remain in $[0;1]$, i.e. for all $i,j\in\llbracket 1;N\rrbracket$, $C(s_{i}(t),s_{j}(t),|x_{i}-x_{j}|)\in[0;1]$.

These conditions can be met if we set some conditions on the support of the initial distribution $\mu_{0}$. The next proposition gives sufficient conditions on the support of $\mu_{0}$ for these constraints to be verified.

Proposition 2**.**

Let $\displaystyle\mathcal{D}=\{(s,x,y,S,\gamma)\in[s_{m};+\infty[\times\mathbb{R}\times\mathbb{R}\times[s_{m};+\infty[\times\mathbb{R}_{+}|s_{m}<S\leq s_{m}e^{R_{M}},~{}s_{m}<s\leq S\}$. Let $\mu_{0}$ be a probability over $\mathbb{R}^{5}$, i.e. $\mu_{0}\in\mathcal{P}(\mathbb{R}^{5})$, such that the support of $\mu_{0}$ is included in the interior of domain $\mathcal{D}$. Let $Z_{N}^{0}=(s_{i}^{0},x_{i},y_{i},S_{i},\gamma_{i})_{1\leq i\leq N}$ be a random variable of distribution $\mu_{0}^{\otimes N}$ and $t\in\mathbb{R}_{+}\mapsto s_{1:N}(t,Z_{N}^{0})$ the solution of system (9) with initial configuration $Z_{N}^{0}$. Then we have almost surely that for all time $t\in\mathbb{R}_{+}$, $(s_{i}(t,Z_{N}^{0}),x_{i},y_{i},S_{i},\gamma_{i})_{1\leq i\leq N}\in\mathring{\mathcal{D}}^{N}$, the interior of domain $\mathcal{D}^{N}$.

The proof of this proposition can be found in appendix 6.2. It is based on the fact that within domain $\mathcal{D}$, the evolution of each plant size is bounded between two growth rates, ensuring the size to remain within a biologically consistent interval.

[TABLE]

The trajectories associated respectively to the upper and the lower bound remain within the domain $\mathcal{D}$ when the initial condition is generated by a $\mu_{0}$ satisfying the assumptions of proposition 2.

For the sake of clarity, we give also an example of initial distribution $\mu_{0}$, that is the source of heterogeneity and randomness in the system represented by equation (9). Let $(s^{0},x,y,S,\gamma)$ be a random variable of distribution $\mu_{0}$. We have chosen the distribution $\mu_{0}$ such that the positions of plants are mutually independent, but with a spatial pattern on parameters $\gamma$ and $S$. In what follows, $\mathcal{U}([a;b])$ is the notation for the uniform distribution over the segment $[a;b]$.

[TABLE]

This initial distribution $\mu_{0}$ implies that the plants are evenly distributed over the square $[0;L]^{2}$, that the plants with large values of $x$ are likely to be tall, and the ones with high values of $y$ are likely to grow fast. In this example, the initial distribution $\mu_{0}$ is not absolutely continuous with respect to the Lebesgue measure. However, the marginal distribution of the intrisic parameters $\theta=(x,y,S,\gamma)$ is absolutely continous with respect to $\lambda^{\otimes 4}$, the Lebesgue measure over $\mathbb{R}^{4}$. Let $p_{0}^{\theta}:\mathbb{R}^{4}\rightarrow\mathbb{R}_{+}$ be the density of $\theta$.

[TABLE]

We can therefore simulate the population model by first drawing samples from the distribution $\mu_{0}$, and finally by solving the differential system (9) using standard numerical methods. In our case, we have used a simple Euler explicit method with a time step of $\Delta t=0.1$ day. The following table gives the configuration used for the simulation.

A visualization of the impact of competition on plant growth is presented in figure 2. Depending on its position and on its intrisic parameters $S$ and $\gamma$, the response of a plant to competition with the rest of the population can varie significantly. In the middle of the domain $[0;L]^{2}$, a plant is more subject to competition than a plant at the boundary, since it is surrounded by more competitors. The evolution of the size over the time depends also on the number $N$ of individuals. We can notice that the responses are quite different from $N=11$ to $N=101$, but there are very little changes from $N=101$ to $N=501$. This convergence constitutes a first visualization of the MFL distribution : as $N$ increases, the finite sample of competitors behaves more and more as a deterministic continuum. The next section gives a formal proof of this statement.

We can consider the change of variable $\displaystyle r=\log\left(\frac{s}{s_{m}}\right)$, so that the state variable $r$ lies in the vector space $\mathcal{X}=\mathbb{R}$, and $d_{\mathcal{X}}=1$. This change of variable is also applied to the initial distribution $\mu_{0}$, so that the marginal initial distribution of the state is for now on related to $r$ variable $\mu_{0}^{r}=\delta_{r^{0}}=\delta_{\log(s^{0}/s_{m})}$. The parameter space $\Theta$ can be chosen as $\Theta=[0;L]^{2}\times[s_{m};S_{M}]\times[\gamma_{m};\gamma_{M}]$. The reference vector to define a norm over $\mathcal{Z}$ can be chosen as $z^{*}=(1,L,L,s_{m},\gamma_{m})$.

[TABLE]

The interaction function $g$ has the expression of function $g_{r}$ defined in equation (38). Over $\mathcal{Z}$, all four assumptions are satisfied by function $g$. Possible choices of constants $K_{1},K_{2},K_{3},K_{4}$ are given below.

[TABLE]

3 Derivation of the mean-field limit

This section follows similar steps as in Golse [2013a] to establish the MFL distribution associated to system (4). We start by proving that system (4) implies a transport equation verified by the empirical measure of the population. From this transport equation, we derive the expression of the MFL transport equation monitoring the dynamics of the MFL distribution. Finally, the connection between the two transport equations is given by Dobrushin stability, which implies also chaos propagation.

3.1 Properties of the population empirical measure

The system (4) has an unique global solution if the interaction function $g$ satisfies assumptions (A1) and (A2). Let $Z_{N}^{0}=(z_{i}^{0})_{1\leq i\leq N}=((X_{i}^{0},\theta_{i}))_{1\leq i\leq N}$ be an initial configuration of the system (4). We introduce $t\in\mathbb{R}_{+}\mapsto Z_{N}(t,Z_{N}^{0})=(z_{i}(t,Z_{N}^{0}))_{1\leq i\leq N}=((X_{i}(t,Z_{N}^{0}),\theta_{i}))_{1\leq i\leq N}\in\mathcal{Z}^{N}$ the global solution of the system. The empirical measure of the population is defined as the map

[TABLE]

In the above equation, we use the notation $\forall z\in\mathcal{Z},\delta_{z}$ is the Dirac distribution centered at $z$, i.e. the distribution of the random variable which is almost surely constant equal to $z$. The empirical measure of the population is a dynamical probability distribution. Sampling this distribution at a fixed time $t$ corresponds to choose an individual uniformly over the population (with probability $\displaystyle\frac{1}{N}$). Interestingly, the empirical measure describes exhaustively the dynamics of the whole population, while remaining in a space $\mathcal{P}(\mathcal{Z})$, which does not depend of the population size $N$. However, there is a loss of information from vector $Z_{N}(t,Z_{N}^{0})$, where all individuals are labelled by indices in $\llbracket 1;N\rrbracket$, to the measure $\mu[t,Z_{N}^{0}]$ where all individuals are not distinguishable. In other words, a visualization of vector $Z_{N}^{0}$ in the phase space $\mathcal{Z}$ would be a cloud of points with all different colors, whereas a visualization of $\mu[t,Z_{N}^{0}]$ would be the same cloud of points with a single color. This indistinction of the individuals is a first step towards the mean-field limit, where individuals are punctual parts of a continuum.

Let us characterize the dynamics of $\mu[t,Z_{N}^{0}]$ using the system (4). We observe the dynamics of a probability measure through its action on test functions, which in our case is the functional space $\mathcal{C}^{1}_{0}(\mathbb{R}_{+}\times\mathcal{Z}\rightarrow\mathbb{R})$.

[TABLE]

In particular, the test functions considered in this article are bounded over their domain, and have bounded derivatives. We call action of $\mu[t,Z_{N}^{0}]$ on a test function $\varphi\in\mathcal{C}^{1}_{0}(\mathbb{R}_{+}\times\mathcal{Z}\rightarrow\mathbb{R})$ the dual pairing of $\mu[t,Z_{N}^{0}]$ and $\varphi$ or the expectaction of random variable $\varphi(t,z)$ where $z$ is a random variable of distribution $\mu[t,Z_{n}^{0}]$.

[TABLE]

The time evolution of $t\in\mathbb{R}_{+}\mapsto\mu[t,Z_{N}^{0}]$ can be studied by considering the differential equation satisfied by $t\in\mathbb{R}_{+}\mapsto\displaystyle\int_{\mathcal{Z}}\varphi(t,z)\mu[t,Z_{N}^{0}](\mathrm{d}z)$ for any test function $\varphi$. So let us express the derivative $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathcal{Z}}\varphi(t,z)\mu[t,Z_{N}^{0}](\mathrm{d}z)$ as an action of $\mu[t,Z_{N}^{0}]$ on some function depending on $\varphi$ and on the interaction function $g$.

[TABLE]

[TABLE]

In the above equation, we can interpret the term

[TABLE]

as the velocity field associated to the system (4), i.e. the one assigning to each individual its velocity according to the current state of the whole population333Similar developments can be found in Carrillo et al. [2010]. $\mathcal{G}_{N}$ is said to be a non-local velocity field, because it depends on the probability measure describing the state of the population, $\mu[t,Z_{N}^{0}]$ in this case. The velocity field can be associated to a conservative transport equation, having a formulation quite similar, in its principle at least, to Vlasov equations, where the velocity field depends on the unknown density (see Golse [2003], section 1.1.1).

[TABLE]

where $\mathrm{div}_{X}$ is the divergence operator with respect to state variable $X$, i.e. for any continuously differentiable map $F:\mathcal{X}\rightarrow\mathcal{X}$, $\mathrm{div}_{X}F(X)=\displaystyle\sum_{i=1}^{d_{X}}\frac{\partial F_{i}(X)}{\partial X_{i}}$, and $f(t,\cdot)\lambda^{\otimes d_{Z}}$ is the probability measure of density $f(t,\cdot):\mathcal{Z}\rightarrow\mathbb{R}_{+}$. The equation (17) is a weak formulation of equation (19), which is formally defined in definition 2. As the weak formulation deals with trajectories taking values in the space of probability measures, we need to introduce the Wasserstein distance to quantify the regularity of these trajectories.

Definition 1**.**

Let $\mu_{1},\mu_{2}\in\mathcal{P}_{1}(\mathcal{Z})$. Let $\Pi(\mu_{1},\mu_{2})$ the set of couplings of $\mu_{1}$ and $\mu_{2}$, i.e. the set of probability distributions having its first and second marginals equal to $\mu_{1}$ and $\mu_{2}$ respectively.

[TABLE]

The Wasserstein distance of first order between $\mu_{1}$ and $\mu_{2}$ is defined by

[TABLE]

or, equivalently, the Wasserstein distance of first order has a dual representation (Kantorovich and Rubinstein [1958])

[TABLE]

with $\mathcal{C}_{L}(\mathcal{Z})$ being the space of Lipschitz-continuous functions over $\mathcal{Z}$, taking values in $\mathbb{R}$, and $\mathrm{Lip}(\varphi)$ being the Lipschitz constant of $\varphi\in\mathcal{C}_{L}(\mathcal{Z})$.

Definition 2**.**

Let $\mathcal{G}:(\mu,z)\in\mathcal{P}_{1}(\mathcal{Z})\times\mathcal{Z}\rightarrow\mathcal{X}$ a non-local velocity field and $\mu_{0}\in\mathcal{P}_{1}(\mathcal{Z})$ a probability measure having first order moment, i.e. $\displaystyle\int_{\mathcal{Z}}|z|\mu_{0}(\mathrm{d}z)<+\infty$. We say that the trajectory $t\in\mathbb{R}_{+}\mapsto\mu[t]\in\mathcal{P}_{1}(\mathcal{Z})$ is a measure solution of the transport equation of velocity field $\mathcal{G}$ and of initial condition $\mu_{0}$ if

the trajectory $t\in\mathbb{R}_{+}\mapsto\mu[t]\in\mathcal{P}_{1}(\mathcal{Z})$ is continuous for the metric $W_{1}$. 2. 2.

for all test function $\varphi$, for all time $t\in\mathbb{R}_{+}$

[TABLE]

Proposition 3**.**

Let $g:\mathcal{Z}^{2}\rightarrow\mathcal{X}$ satisfying assumptions (A1) and (A2) and $Z_{N}^{0}=(z_{i}^{0})_{1\leq i\leq N}=((X_{i}^{0},\theta_{i}))_{1\leq i\leq N}\in\mathcal{Z}^{N}$. Then the empirical measure $t\in\mathbb{R}_{+}\mapsto\mu[t,Z_{N}^{0}]$, defined in equation (14) is a measure solution to the transport equation of velocity field $\mathcal{G}_{N}$, defined in equation (18), and of initial condition $\mu[0,Z_{N}^{0}]=\displaystyle\frac{1}{N}\sum_{i=1}^{N}\delta_{z_{i}^{0}}$.

This proposition summarizes the equation (17). It is also necessary to check the continuity of the trajectory $t\in\mathbb{R}_{+}\mapsto\mu[t,Z_{N}^{0}]\in\mathcal{P}_{1}(\mathcal{Z})$ is continuous for the metric $W_{1}$. This continuity is directly given by the continuity of the solution of the system (4) (see appendix section 7).

The transport equation (17) satisfied by the empirical measure leads to the transport equation describing the dynamics of the population with an infinite number of individuals by taking the limit $N\rightarrow+\infty$. The resulting equation, obtained informally, is referred as the MFL transport equation, and its eventual solution is the MFL distribution. Subsection 3.2 solves the MFL transport equation and proves the existence and uniqueness of the MFL distribution for system (4). Subsection 3.3 studies different aspects of the convergence towards the MFL distribution.

3.2 Study of the mean-field equation

This subsection gives a characterization of the MFL distribution as the unique solution of a non-local transport equation obtained as the limite case of transport equation (17). Let us assume informally that, for some metric over the space of probability measures (namely the Wasserstein distance, see subsection 3.3), the empirical measure $\mu[t,Z_{N}^{0}]$ has a limit $\mu[t]$ when $N\rightarrow+\infty$. Then it is reasonable to think that, for some other metric, the velocity field $\mathcal{G}_{N}(\mu[t,Z_{N}^{0}],\cdot)$ converges towards a velocity field depending on $g$ and $\mu[t]$. The only expression this velocity field can reasonably have, when $N\rightarrow+\infty$ in equation (18), is

[TABLE]

So if the MFL distribution exists, it has to be a measure solution of the transport equation of velocity field $\mathcal{G}$ and of initial condition $\mu_{0}$, as it is reminded in subsection 3.3 that $\mu[0,Z_{N}^{0}]$ converges towards $\mu_{0}$ for the Wasserstein distance. The MFL distribution is therefore an eventual trajectory $t\in\mathbb{R}_{+}\mapsto\mu[t]\in\mathcal{P}_{1}(\mathcal{Z})$, continuous for the metric $W_{1}$, such that for all test function $\varphi$ and for all time $t\in\mathbb{R}_{+}$

[TABLE]

Drawing largely on Golse [2013a], we would like to use the characteristic flow method to prove the existence and uniqueness of the solution to the equation (23). The characteristic flow method is a classical idea to study a transport equation : the transport PDE describes the dynamics, while the characteristic flow equation describes the motion of a single particle, subject to the same velocity field. Let $t\in\mathbb{R}_{+}\mapsto X_{\infty}(t,X,\theta)\in\mathcal{X}$ be the trajectory of a particle immersed in velocity field $\mathcal{G}$, of initial configuration $(X,\theta)\in\mathcal{Z}$.

[TABLE]

As $\mu[t]$ is unknown, we cannot evaluate the derivative $\displaystyle\frac{\mathrm{d}X_{\infty}}{\mathrm{d}t}$, except at $t=0$, when $\mu[0]=\mu_{0}$. However, there is a strong connection between the MFL distribution $\mu[t]$ and the flow $X_{\infty}(t,\cdot)$, as $\mu[t]$ describes the state of a population which is composed of infinite number of particles $(X_{\infty}(t,X^{\prime},\theta^{\prime}),\theta^{\prime})$, whose initial configuration is given by $\mu_{0}$. In other words, we can look at the interaction with the rest of the population not as an average over all the possible states at time $t$, as it is the case in equation (24), but as an average over all initial configurations.

[TABLE]

There is only a single unknown in the above functional equation, which is the flow $X_{\infty}:\mathbb{R}_{+}\times\mathcal{Z}\rightarrow\mathcal{X}$. The next theorem shows that this characteristic flow is well defined.

Theorem 1**.**

Let $\mu_{0}\in\mathcal{P}_{2}(\mathcal{Z})$ a probability measure having second order moment, i.e. $\displaystyle\int_{\mathcal{Z}}|z|^{2}\mu_{0}(\mathrm{d}z)<+\infty$, and $g:\mathcal{Z}^{2}\rightarrow\mathcal{X}$ satisfying assumptions (A1) and (A2). Then there exists an unique flow such that

$\forall t\in\mathbb{R}_{+},~{}\forall(X,\theta)\in\mathcal{Z},~{}\displaystyle\int_{\mathcal{Z}}|g(X_{\infty}(t,X,\theta),\theta,X_{\infty}(t,X^{\prime},\theta^{\prime}),\theta^{\prime})|\mu_{0}(\mathrm{d}X^{\prime},\mathrm{d}\theta^{\prime})<+\infty$** 2. 2.

$\forall(X,\theta)\in\mathcal{Z},~{}t\in\mathbb{R}_{+}\mapsto X_{\infty}(t,X,\theta)$* is continuously differentiable.* 3. 3.

$\forall(X,\theta)\in\mathcal{Z}$,

[TABLE]

The above functional equation can be seen as a continuous version of system (4). Formally, this equation is a differential equation with an initial condition being a probability measure and having trajectories in vector space $\mathcal{X}$. This theorem can therefore be proved by following exactly the same steps as for Cauchy-Lipschitz theorem with traditional differential equations. A common proof of Cauchy-Lipschitz theorem is based on fixed point theorem within a functional Banach space. The functional space where the flow solution $X_{\infty}$ lies must be complete, as the fixed point theorem recquires the convergence of any Cauchy sequence. The next lemma introduces the functional space used in the proof of theorem 1.

Lemma 1**.**

(Golse [2013a]) Let $\mathcal{Y}$ be the functional space defined by

[TABLE]

Then $\mathcal{Y}$ is a Banach space for the metric $f\in\mathcal{Y}\mapsto\|f\|_{\mathcal{Y}}=\displaystyle\sup_{z\in\mathcal{Z}}\frac{|f(z)|}{1+|z|}$.

The proof of the above theorem is given in appendix 8.1, and it follows the same steps than the proof of Cauchy-Lipschitz theorem for ordinary differential equations: local existence and uniqueness, existence of a maximal solution, uniqueness of the maximal solution and finally definition over $\mathbb{R}_{+}$ of the maximal solution. The assumption (A1) on $g$ is important to ensure that the flow solution is defined over $\mathbb{R}_{+}$ and also the stability within the functional space $\mathcal{Y}$. Besides, the control of the Lipschitz factor in assumption (A2) can be relaxed, as long as the Lipschitz factor is compensated by the initial distribution. For instance, if the Lipschitz factor in (A2) is $K_{2}(1+|X|^{n}+|X^{\prime}|^{n})$ for some $n>0$ instead of $K_{2}(1+|X|+|X^{\prime}|)$, then similar reasoning can be carried out to obtain the existence and uniqueness of $X_{\infty}$, if $\mu_{0}$ is chosen in $\mathcal{P}_{n+1}(\mathcal{Z})$, i.e. such that $\displaystyle\int_{\mathcal{Z}}|z|^{n+1}\mu_{0}(\mathrm{d}z)<+\infty$.

In the case of Schneider model, the flow solution $t\in\mathbb{R}_{+}\mapsto r_{\infty}(t,.)$ satisfies the following equation

[TABLE]

If $\mu_{0}$ is the distribution defined in equation (10), we have for all time $t\in\mathbb{R}_{+}$ and for all $(r,\theta)\in\mathcal{Z}$,

[TABLE]

where $r^{0}=\log\left(\frac{s^{0}}{s_{m}}\right)$. This relation holds because the marginal distribution of the initial state is a Dirac distribution centered at $r^{0}$.

The characteristic flow (26) leads to the unique solution of the mean-field transport equation (23), which appears as the pushforward probability measure of the initial distribution $\mu_{0}$ by the map $X_{\infty}(t,.)$. This result is used qualitatively to derive equation (26) from equation (24).

Corollary 1**.**

Let $\mu_{0}\in\mathcal{P}_{2}(\mathcal{Z})$, g satisfying assumptions (A1), (A2) and (A3), and $z^{0}=(X^{0},\theta)$ a random variable of distribution $\mu_{0}$. Then the unique measure-solution to the transport equation (23) is $t\in\mathbb{R}_{+}\mapsto\mu[t]\in\mathcal{P}_{2}(\mathcal{Z})$ where for all $t\in\mathbb{R}_{+}$, $\mu[t]$ is the probability distribution of $z^{t}=(X_{\infty}(t,z^{0}),\theta)$.

The proof of this corollary (appendix, section 8.2) is a generalization of the method of characteristic flows, classically used in the field of hyperbolic PDE. The proof that the pushforward measure is effectively a measure solution of the mean-field transport equation (23) is mainly based on the change-of-variable formula, stating that for every test function $\phi$, we have $\displaystyle\int_{\mathcal{Z}}\phi(t,z)\mu[t](\mathrm{d}z)=\int_{\mathcal{Z}}\phi(t,X_{\infty}(t,X,\theta),\theta)\mu_{0}(\mathrm{d}X,\mathrm{d}\theta)$. The continuity of the trajectory $t\mapsto\mu[t]$ is therefore implied by the continuity of $t\mapsto\|X_{\infty}(t,.)\|_{\mathcal{Y}}$. Besides, the proof of uniqueness requires additional assumptions on the regularity of the interaction function $g$, namely assumption (A3). It was added for the sake of brevity, but it seems that this assumption can be avoided by adding technical developments using an argument of density of the space of test functions. This regularity enables to prove that the flows $X^{\nu}$ associated to the velocity field $\mathcal{G}(\nu[t],z)$, where $t\mapsto\nu[t]$ is a fixed measure-trajectory, are continuously differentiable with respect to their arguments. This implies also the regularity of the solution $(t,z)\mapsto\varphi(t,z)$ to the following transport equation

[TABLE]

As $\varphi$ is regular, it can be used as a test function itself. If $t\mapsto\nu[t]$ is a measure-solution of (23), it can be shown that $t\mapsto\int_{\mathcal{Z}}\varphi(t,z)\nu[t](\mathrm{d}z)$ is constant. This implies in particular that $\nu[t]$ is the pushforward measure of $\mu_{0}$ by the map $X^{\nu}(t,.)$. It follows that $X^{\nu}$ and $X_{\infty}$ satisfies the same equation (26) and we conclude by uniqueness of the characteristic flow provided by theorem 1.

Alternative proofs of existence and uniqueness of the MFL distribution can be found in Lagoutière and Vauchelet [2017] or in Bolley et al. [2011], with weaker assumptions made on the velocity field. Lagoutière and Vauchelet [2017] uses Filipov characteristics and compactness arguments to solve the transport equation associated to a bounded velocity field with a finite set of discontinuities. The velocity field considered in Bolley et al. [2011] is not globally Lipschitz continuous, as in our case. The resolution of an equation similar to (24) follows an iterative procedure : the measure trajectory $t\mapsto\mu^{n}[t]$ is fixed at iteration $n$, and is used to compute the characteristic flow $X^{\mu^{n}}$, by solving a standard differential equation ; the distribution $t\mapsto\mu^{n+1}[t]$ chosen at the next iteration is the pushforward measure of $\mu_{0}$ by the characteristic flow $X^{\mu^{n}}$.

Let us now consider the case where the initial distribution $\mu_{0}$ is absolutely continuous with respect to the Lebesgue measure $\lambda^{\otimes d_{z}}$. We denote by $p_{0}:\mathcal{Z}\rightarrow\mathbb{R}_{+}$ its associated probability density. $p_{0}$ can be factorized in two terms using the chain rule.

[TABLE]

It follows from proposition 4 that, in this case, the MFL distribution $\mu[t]$ is absolutely continuous for all time $t\in\mathbb{R}_{+}$, and that the associated density $p_{t}:\mathcal{Z}\rightarrow\mathbb{R}_{+}$ is given by change-of-variable formula.

[TABLE]

In the above equation, for all $t\in\mathbb{R}_{+}$ and $\theta\in\Theta$, $X\in\mathcal{X}\mapsto X_{\infty}^{-1}(t,X,\theta)$ is the inverse function of $X\in\mathcal{X}\mapsto X_{\infty}(t,X,\theta)$ and $\displaystyle\frac{\partial X_{\infty}^{-1}}{\partial X}(t,X,\theta)$ is the Jacobian matrix of this function.

The fact that $X_{\infty}(t,.,\theta)$ is a one-to-one map from $\mathcal{X}$ to $\mathcal{X}$ is a consequence of the flow property. This property considers a variation of the initial time in equation (26). There exists an unique map $(t,t_{0},z)\mapsto X_{\infty}(t,t_{0},z)$ satisfying

[TABLE]

By unicity, we have that for all $t,t_{0}\in\mathbb{R}_{+}$ and $(X,\theta)\in\mathcal{Z}$, $X_{\infty}(t,t_{0},X_{\infty}(t_{0},t,X,\theta),\theta)=X$. It follows that $X_{\infty}^{-1}(t,X,\theta)=X_{\infty}(0,t,X,\theta)$. The differentiability of the map $X\in\mathcal{X}\mapsto X_{\infty}^{-1}(t,X,\theta)$ is a consequence of assumption (A3).

In the case of Schneider model, the MFL distribution $\mu[t]$ is the law of the random variable $(r_{\infty}(t,r^{0},\theta),\theta)$ with $\theta\sim p_{0}^{\theta}$. Therefore $\mu[t]$ is entirely determined by the map $(t,\theta)\in\mathbb{R}_{+}\times\theta\mapsto r_{\infty}(t,r^{0},\theta)$. This is due to the fact that the initial state is constant over the whole population, equal to $r^{0}$. We can notice that this situation seems much simpler than the case where the marginal density of the initial state is absolutely continuous : here, we only need to compute the characteristic flow $r_{\infty}$, and we do not need to compute its inverse function and its derivative. Section 4 describes a methodology to approximate the characterstic flow $r_{\infty}$ and therefore to sample the MFL distribution in the specific case of Schneider model.

3.3 Dobrushin stability and propagation of chaos

The relation between the microscopic level, represented by the empirical measure of the population, and the mean-field level, represented by solution of (23), is mainly based on a convergence of the initial empirical distribution $\mu[0,Z_{N}^{0}]$ towards the initial distribution $\mu_{0}$ as N $\rightarrow+\infty$ and on the fact that this results can be extended at all time $t\in\mathbb{R}_{+}$, i.e. the same type of convergence is verified by $\mu[t,Z_{N}^{0}]$ towards $\mu[t]$. The convergence discussed here is the one associated to the metric $W_{1}$, whose expression is recalled in definition 1, and which metrizes the weak convergence in the space of probability distribution (see corollary 6.13 in Villani [2008]). According to Varadarajan [1958], if $(z_{n}^{0})_{n\in\mathbb{N}}$ is a sequence of independent random variables of distribution $\mu_{0}$, and if for all $N>1$ $Z_{N}^{0}=(z_{1}^{0},...,z_{N}^{0})$, we have that $\displaystyle\lim_{N\rightarrow+\infty}W_{1}\left(\mu[0,Z_{N}^{0}],\mu_{0}\right)=0$ almost surely in $\mathbb{P}$. This means that there exists $\Omega^{*}\in\mathcal{F}$ such that $\mathbb{P}(\Omega^{*})=1$ and for all $\omega\in\Omega^{*}$, we have for all $\phi\in C_{L}(\mathcal{Z})$ any Lipschitz continuous function that

[TABLE]

A concise proof of this result can also be found in Golse [2013a] (theorem 3.3.5). This result is a consequence of the strong law of large numbers and of the fact that the space of continuous functions with compact support over $\mathbb{R}^{d_{z}}$ is separable. The rate of convergence of the random variable $W_{1}(\mu[0,Z_{N}^{0}],\mu_{0})$, along with Wasserstein distance of higher orders, is a well-documented topic in the literature. Dudley [1969] stated that in the case where $\mu_{0}$ is absolutely continuous with respect to the Lebesgue measure, i.e. can be associated to a probability density $f_{0}:\mathcal{Z}\rightarrow\mathbb{R}_{+}$, and if $d_{\mathcal{Z}}\geq 2$, then there exists a constant $C(\mu_{0})>0$ such that for all $N\in\mathbb{N}^{*}$, $\mathbb{E}_{Z_{N}^{0}\sim\mu_{0}^{\otimes N}}W_{1}\left(\mu[0,Z_{N}^{0}],\mu_{0}\right)\leq C(\mu_{0})N^{-1/d_{\mathcal{Z}}}$. In the example of the Schneider model, the chosen initial density described in equation (10) is not absolutely continuous with respect to the Lebesgue measure over $\mathcal{Z}$, since the marginal distribution $\mu_{0}^{r}$ of the variable $r^{0}$ is reduced to a Dirac distribution $\delta_{r^{0}}$, representing the fact hat all plants have the same size $s^{0}=s_{m}e^{r^{0}}$ initially. However, the marginal $\mu_{0}^{\theta}$ of variable $\theta$ is associated to a density over $\Theta$, so the upper-bound of Dudley can be rewritten as $\mathbb{E}_{Z_{N}^{0}\sim\mu_{0}^{\otimes N}}W_{1}\left(\mu[0,Z_{N}^{0}],\mu_{0}\right)\leq C(\mu_{0}^{\theta})N^{-1/d_{\Theta}}=C(\mu_{0}^{\theta})N^{-1/4}$. Faster convergence rates can be obtained in the case where probability measure $\mu_{0}$ is less regular. We can quote notably Weed and Bach [2017] in the case where $\mu_{0}$ is compactly-supported, and Lei [2018] for a generalization to unbounded metric spaces.

If the random variable $z^{0}\sim\mu_{0}$ has at least one of its component with a probability density, then the weak and almost sure convergence of $\mu[0,Z_{N}^{0}]$ towards $\mu_{0}$ means visually that the point cloud $(z_{i}^{0})_{1\leq i\leq N}$ is more and more alike the continuous set represented by measure $\mu_{0}$. In what follows, the argument of Dobrushin stability (see proposition 4 in Dobrushin [1979] or theorem 3.3.3 in Golse [2013a]) is used to prove that for all time $t\in\mathbb{R}_{+}$ $\displaystyle\lim_{N\rightarrow+\infty}W_{1}(\mu[t,Z_{N}^{0}],\mu[t])=0$ almost surely.

Theorem 2**.**

Let $\mu_{0}\in\mathcal{P}(\mathcal{Z})$ be a probability measure having a compact support, such that the support is included in $\mathcal{B}(0,R_{0})=\{z\in\mathcal{Z}||z|\leq R_{0}\}$ for some $R_{0}>0$. Let $(z_{n}^{0})_{n\in\mathbb{N}^{*}}$ be a sequence of independent random variables of distribution $\mu_{0}$ and $\forall N>1,Z_{N}^{0}=(z_{1}^{0},...,z_{N}^{0})$. There exists $\Omega^{*}\in\mathcal{F}$ such that $\mathbb{P}(\Omega^{*})=1$ and such that $\forall\omega\in\Omega^{*},~{}\displaystyle\lim_{N\rightarrow+\infty}W_{1}(\mu[0,Z_{N}^{0}(\omega)],\mu_{0})=0$. We introduce $\mu[t]$ the solution of problem (23) for an interaction function $g$ satisfying assumptions (A1), (A2), (A3) and (A4). Then we have

[TABLE]

The argument of Dobrushin consists in deriving an upper bound of $W_{1}(\mu[t,Z_{N}^{0}],\mu[t])$ depending on $W_{1}(\mu[0,Z_{N}^{0}],\mu_{0})$, which holds for all initial configuration $Z_{N}^{0}\in\mathcal{Z}^{N}$. The derivation of the upper bound is exactly a generalization of Grönwall lemma to characteristic flows of the type of $X_{\infty}$, i.e. solutions of differential equations taking values in space $\mathcal{X}$ and having as initial condition a probability measure over $\mathcal{Z}$. Indeed, we can prove that for all initial configuration $Z_{N}^{0}\in\mathcal{Z}^{N}$, we have

[TABLE]

The functions $E_{\mu_{0}}$ and $F_{\mu_{0}}$ appearing in the previous inequality depend on the first and second moments $M^{1}_{\mu_{0}}=\displaystyle\int_{\mathcal{Z}}|z|\mu_{0}(\mathrm{d}z)$ and $M^{2}_{\mu_{0}}=\displaystyle\int_{\mathcal{Z}}|z|^{2}\mu_{0}(\mathrm{d}z)$ of the distribution $\mu_{0}$. The other functions appearing in the inequality (32) are such that $\forall t\in\mathbb{R}_{+},\forall\omega\in\Omega^{*},~{}\displaystyle\lim_{N\rightarrow+\infty}\epsilon_{1}(t,\mu[t,Z_{N}^{0}])=\lim_{N\rightarrow+\infty}\epsilon_{2}(t,\mu[t,Z_{N}^{0}])=0$.

Historically, Dobrushin [1979] introduced this methodology to obtain uniqueness results on solutions of Vlasov equations. In this article, the studied interaction functions are globally Lipschitz continuous, and the author does not resort to Grönwall lemma. With the same assumptions, a proof of Dobrushin stability was suggested by Golse [2013a], theorem 1.4.3, making clear use of Grönwall lemma. In Lagoutière and Vauchelet [2017], the proposition 1 gives a quite similar contraction estimate, in the case where the transition function is expressed as a convolution product with the mean-field limit measure. In all aforementionned works, the transport functions $\mathcal{G}_{N}$ at microscopic level have the same expression as the transport function $\mathcal{G}$ at macroscopic level, and the physical models do not recquire to exclude the interaction of a particle with itself, notably thanks to a property of anti-symmetry of the underlying potential. In our case, the transition function is only assumed to be locally Lipschitz continuous, but this difficulty is bypassed by assuming that the i-nitial distribution $\mu_{0}$ has a compact support. The obtained upper-bound of $W_{1}\left(\mu[t,Z_{N}^{0}],\mu[t]\right)$ in (32) is a much faster increasing function than in Golse [2013a]. The assumption on global Lipschitz continuity of the function $g$ leads to a factor of order $e^{Kt}$ for some constant $K$, whereas the assumptions on quadratic variations of the functions, namely (A2) and (A4), leads to a factor of order $\exp\left(e^{Kt}\right)$ for some constant $K$, because of two subsequent applications of Grönwall lemma (see the proof in appendix 9.2). Needless to say that the upper bound in (32) seems far from being optimal.

The next corollary uses the argument of Dobrushin stability to show the relation between the solution of the microscopic system (4) and the MFL characteristic flow.

Corollary 2**.**

With the same assumptions as in theorem 31, we consider the sequence of random variables $(z_{n}^{0})_{n\in\mathbb{N}^{*}}$ independent and of distribution $\mu_{0}$. For all $N>1$, we define $Z_{N}^{0}=(z_{1}^{0},z_{2}^{0},...,z_{N}^{0})\in\mathcal{Z}^{N}$ the initial configuration of the system (4) and $t\in\mathbb{R}_{+}\mapsto(X_{1}(t,Z_{N}^{0}),...,X_{N}(t,Z_{N}^{0}))$ the solution of the system (4). Then we have

[TABLE]

The above results provides a more visual intuition of the asymptotic link between the microscopic level of system (4) and the mean-field limit. The trajectories obtained by solving system (4) are more and more alike the trajectories given by the MFL characteristic flow $X_{\infty}$. A generalization to any sub-group of fixed size within the population can also be obtained. Indeed, for $k\in\mathbb{N}^{*}$ and for $N>k$, any sub-group of size $k$ $(X_{i_{1}}(t,Z_{N}^{0}),...,X_{i_{k}}(t,Z_{N}^{0}))$, with $i_{1},...,i_{k}$ being distinct integers in $\llbracket 1;N\rrbracket$, has the same distribution as $(X_{1}(t,Z_{N}^{0}),...,X_{k}(t,Z_{N}^{0}))$ by symmetry. According to the previous corollary, the almost sure convergence for a single individual can be generalized to any sub-group of size $k$.

[TABLE]

The limit distribution of the sequence of random variables $((X_{1}(t,Z_{N}^{0}),...,X_{k}(t,Z_{N}^{0}))_{N>k}$ is factorized and is exactly $\mu[t]^{\otimes k}$, as $(X_{\infty}(t,z_{1}^{0}),...,X_{\infty}(t,z_{k}^{0}))\sim\mu[t]^{\otimes k}$. For finite $N$ and for $t>0$, the random variables $X_{1}(t,Z_{N}^{0}),...,X_{k}(t,Z_{N}^{0})$ are strongly interdependent. At the limit $N\rightarrow+\infty$, the individuals are independent. More accurately, if one focuses on a finite group of individuals, while the rest of the population is increasing towards infinity, then these observed individuals have independent trajectories in the probabilistic sense. Their distribution is said to be asymptotically factorized. An alternative proof of the phenomenon of chaos propagation is given in Golse [2013b], section 1.6. This proof is based on a characterization of asymptotically factorized sequence of probability measures (see theorem 1.6.2 in Golse [2013b]).

The phenomenon of chaos propagation may have applications for statistical inference, paving the way for methodologies based on variational Bayes approximation. Let us consider the following example : we aim at studying the dynamics of an heterogeneous crop from the observation of the growth of few dozens of plants. Their growth is assumed to be well represented by a model of the form of Schneider et al. [2006], but some parameters of the interaction function $g$ are unknown. In general, we do not know accurately the exact number $N$ of individuals in the population, but we know that $N$ is much larger than the number of observed individuals. In a Bayesian setting, i.e. when we want to compare prior knowledge and assumptions with field observations, the resulting inference problem is of great difficulty. Among other things, it requires to determine the posterior distribution of the number of individuals in the population444A possible prior for the random variable $N$ would be a Poisson distribution., but also the posterior distributions of all the unobserved individuals, i.e. of their positions and of their characteristics $\gamma$ and $S$. This is clearly intractable for a population having the dimension of a crop. Otherwise, if we make the approximation that the observed individuals are in interaction with an infinity of individuals, which is quite a relevant approximation after all, and that this continuum of individuals is represented by the MFL distribution $\mu[t]$, then the inference problem is significantly simplified : the observed individuals are then mutually independent, and there is no need to extract the information of all the unobserved individuals. Of course, the difficulty is elsewhere : how to simulate the MFL distribution efficiently, so that it can be used within a statistical inference process. The next section gives a first attempt to answer this issue.

4 Simulation of the MFL distribution using Gaussian process regression

In this section, we present a preliminary work on the numerical approximation of the MFL distribution $t\mapsto\mu[t]$, which is defined as the measure-solution of variational problem (23). So it boils down to solve numerically a hyperbolic PDE with non-local velocity. The simulation of solutions of kinetic equations is a well-documented in the literature. Amongst others, we can quote the upwind scheme introduced by Lagoutière and Vauchelet [2017], which consists in a reconstruction of the solution using finite volumes. The reconstruction is piecewise constant over a discretization of the phase space $\mathcal{Z}$. In the case of Schneider heterogeneous population model, the space is of dimension higher than 3, and this makes the discretization of the space a too expensive task on the computational view point. This constraint of the dimension leads rather towards mesh-free methods.

The method suggested here consists in approximating by regression a consistent sequence of reconstructions of the exact characteristic flow. It is therefore a semi-Lagrangian method with an interpolation step. The family of functions used for the interpolation is defined from the interaction function $g$, and takes the form of linear combinations of reproducing kernels. The proof of the consistency of the scheme is an on-going work. However, some numerical tests seems to confirm that this approach is relevant.

For the sake of simplicity, the method is presented through the simulation of the Schneider model. In this case, the MFL distribution is the law of the random variable $(r_{\infty}(t,r^{0},\theta),\theta)$, where $r^{0}\sim\delta_{r^{0}}$ is a constant, where $\theta\sim p_{0}^{\theta}$ is defined in equation (10), and $r_{\infty}$ is the characteristic flow defined by equation (29). By change of variable, we can consider the characteristic flow associated to the size variable $s$, which is defined as the solution of the equation

[TABLE]

So our aim is to approximate the function $(t,\theta)\in\mathbb{R}_{+}\times\Theta\mapsto s_{\infty}(t,s^{0},\theta)$.

A direct resolution of equation (33) using an explicit Euler method, with time discretization $\Delta t>0$, chosen small enough, would lead to a sequence of functions $(s_{n})_{n\in\mathbb{N}}$ defined an induction equation.

[TABLE]

This sequence of functions cannot be computed exactly, as the integral is not analytical. This integral is in fact an expectation with respect to the density $p_{0}^{\theta}$. Let $\omega_{M}=(\theta_{i}^{\omega})_{1\leq i\leq M}$ be sample of the distribution $\mu_{0}^{\theta}$ of density $p_{0}^{\theta}$. We consider the sequence of functions $(s_{n}(.,\omega_{M}))$ defined as the empirical approximation of the sequence $(s_{n})_{n\in\mathbb{N}}$ using the sample $\omega_{M}$.

[TABLE]

It is quite straightforward to prove that for any fixed $n\in\mathbb{N}$, the sequence $(s_{n}(.,\omega_{M}))_{M\in\mathbb{N}^{*}}$ is an almost sure approximation of the characteristic flow at time $n\Delta t$.

[TABLE]

Indeed, the sequence of functions $(s_{n}(.,\omega_{M}))$ is stochastic because of its dependency with respect to the sample $\omega_{M}$. The above convergence is mainly based on the law of large numbers, enabling to prove a uniform almost sure convergence over the space $\Theta$. So $(s_{n}(.,\omega_{M}))_{n\in\mathbb{N}}$ constitutes a simple approximation of the characteristic flow, but it has some limitations. It can only give a local estimation of the function $s_{n}$. Indeed, to compute $s_{100}(\theta,\omega_{M})$ at a given point $\theta_{0}\in\Theta$, then it requires the computation of $s_{99}(\theta_{0},\omega_{M})$, and in turn the computation of $s_{98}(\theta,\omega_{M})$, etc… We cannot know the values of the function $s_{n}(.,\omega_{M})$ outside of the set of points we have decided to observe a priori, from the very initial time $n=0$, this set of observation points including also the sample $\omega_{M}$. For a global approximation of the function $s_{n}$, a grid covering the whole space $\theta$ has to be build, so this boils down exactly to the construction of a mesh, which is to be avoided in our case. An interpolation method is used at this point so that the local information given by some values of $s_{n}(.,\omega_{M})$ could be extended to the whole space $\Theta$.

The basis of functions used for interpolation has been chosen from a qualitative estimation of the correlation between the values of the function $s_{n}(.,\omega_{M})$. According to central limit theorem, the asymptotic covariance matrix of the random variables $s_{1}(\theta_{1},\omega_{M})$ and $s_{1}(\theta_{2},\omega_{M})$ for $\theta_{1}$ and $\theta_{2}$ in $\Theta$ when $M\rightarrow+\infty$ depends on the interaction $g$.

[TABLE]

$\forall d\in\mathbb{N}^{*},$ $\mathcal{N}_{d}(\mu,\Sigma)$ is the normal distribution of mean $\mu$ and of covariance matrix $\Sigma$. In other words, the random vector $(s_{1}(\theta_{1},\omega_{M}),s_{1}(\theta_{2},\omega_{M}))^{\textsf{T}}$ behaves approximately like a Gaussian vector. From this result, we make the approximation that this property holds for all $n\in\mathbb{N}^{*}$555this approximation may not be justified theoretically..

[TABLE]

This reasonable expression of the covariance leads to a choice of interpolation functions being defined from the covariance function, which is by construction a positive kernel.

[TABLE]

This kernel cannot be used per se as the sequence $(s_{n})_{n\in\mathbb{N}}$ is unknown. Another kernel is therefore chosen, but still largely inspired from the above expression. As the sequence $(s_{n})_{n\in\mathbb{N}}$, they are replaced in the above expression by parametric functions that reproduce roughly their variations over the space $\Theta$. More specifically, polynomial functions $(m^{s}_{n})_{n\in\mathbb{N}}$ of degree 2 were chosen to approximate the sequence $(s_{n})_{n\in\mathbb{N}}$.

[TABLE]

where $v_{16}:\mathcal{M}_{4}(\mathbb{R})\rightarrow\mathbb{R}^{16}$ is the canonical bijection between the square matrices $4\times 4$ and the vector of 16 components. The coefficients $(a_{n},b_{n},c_{n})$ are chosen so that the parametric function is close to function $s_{n}$ is the $L^{2}$ sense.

[TABLE]

Equivalently, $(a_{n},b_{n},c_{n})$ is the solution of a linear system expressed with expectations with respect to the density $p_{0}^{\theta}$.

[TABLE]

In this linear system, the functions $s_{n}$ can be replaced by their stochastic approximations $s_{n}(.,\omega_{M})$, and the theoretical mean can be replaced by an empirical mean over the set $\omega_{M}$.

The final expression for the kernel used for the interpolation depends on the sample $\omega_{M}$.

[TABLE]

The theoretical covariance is replaced by an empirical covariance over the sample $\omega_{M}$ and the characteristic flows $s_{n-1}$ are replaced by either the polynomial functions $m^{s}_{n-1}$ either the stochastic approximation $s_{n-1}(.,\omega_{M})$. If the parametric model $m^{s}_{n}$ is not too rough and if $M$ is large, then the above covariance function $k_{n}$ is consistent with the stochastic behaviour of $s_{n}(.,\omega_{M})$, and $k_{n}$ is easy to evaluate over the whole space.

In addition to the values $s_{n}(\omega_{M},\omega_{M})=(s_{n}(\theta_{i}^{\omega},\omega_{M}))_{1\leq i\leq M}$, the function $s_{n}(.,\omega_{M})$ is evaluated over another set of points $\Theta_{1:K}=(\theta_{j})_{1\leq j\leq K}$, called training set, that can also be taken as a sample from the density $p_{0}^{\theta}$. For all $n\in\mathbb{N}$, we extend the values of $s_{n}(\Theta_{1:K},\omega_{M})$ by making the approximation that the values of $s_{n}(.,\omega_{M})$ is a Gaussian process of mean function $\theta\mapsto m^{s}_{n}(\theta)$ and of covariance function $(\theta_{1},\theta_{2})\in\Theta^{2}\mapsto k_{n}(\theta_{1},\theta_{2})$ (cf. Rasmussen [2004] for an introduction to Gaussian processes). In particular, under this approximation, for all $\theta\in\Theta$

[TABLE]

The distribution of $s_{n}(\theta,\omega_{M})$ is given by conditioning with respect to the observed, or rather computed, values of $s_{n}(\Theta_{1:K},\omega_{M})$.

[TABLE]

Therefore, under this approximation, the most probable value of $s_{n}(\theta,\omega_{M})$ knowing the values of $s_{n}(\Theta_{1:K},\omega_{M})$ is given by the mode of the above conditional distribution. This is the reconstruction of the characteristic flow we use to estimate it over the whole space $\Theta$.

[TABLE]

One can notice that more information could have been used to compute the conditional distribution, as we also know the values of $s_{n}(\omega_{M},\omega_{M})$. The reason why the sample $\omega_{M}$ is omitted is just that inverting a matrix of dimension $K$ is cheaper than inverting a matrix of dimension $M+K$.

The relevancy of the reconstruction can be assessed using a test set $\Theta^{t}_{1:K}=(\theta^{t}_{j})_{1\leq j\leq K}$, that can also be a sample drawn from density $p_{0}^{\theta}$. A mean square error is used for this purpose.

[TABLE]

If $J_{n}$ remains relatively small during the iterations, then the reconstruction of $s_{n}(.,\omega_{M})$ is likely to be relevant.

The different steps of the simulation process are summarized in the following algorithm.

The algorithm was run with the parameters values given in table 1 and for $n_{\max}=100$ iterations, with a sample size $M=1000$ and size of training / reconstruction set of $K=100$. Figure 3 displays the evolution of the test error of the reconstruction $\hat{s}_{n}(.,\omega_{M})$, along with the test error associated with the polynomial approximation $m^{s}_{n}$.

Figure 3 shows that both the polynomial approximation and Gaussian process (GP) reconstruction seem to provide a good estimate of the function $\hat{s}_{n}(.,\omega_{M})$, with a relative remaining lower than $12\%$ for the polynomial approximation, and lower than $2\%$ for the GP reconstruction. The error made by the polynomial function increases almost linearly with the iterations, meaning that the shapes of the functions $(\hat{s}_{n}(.,\omega_{M}))_{n\in\mathbb{N}^{*}}$ become more and more complex for large $n$, and the approximation by parabolic functions become more and more rough. As a matter of fact, the GP reconstruction has also an increasing test error, but it still provides a significant improvement with respect to the polynomial approximation.

Once $\hat{s}_{n}(.,\omega_{M})$ is computed, an approximate sample of the marginal distribution of random variable $s^{n}\sim\mu^{s}[n\Delta t]$ can be drawn. The sample is obtained by drawing independent samples $(\theta_{i}^{\prime})_{1\leq i\leq M}$ from density $p_{0}^{\theta}$ and compute the values of the characteristic flow over this sample $(\hat{s}_{n}(\theta_{i}^{\prime},\omega_{M}))_{1\leq i\leq M}$. For $n>0$, the marginal distribution $\mu^{s}[n\Delta t]$ is absolutely continuous with respect to the Lebesgue measure $\lambda^{\otimes 4}$, and the associated density can be estimated by non-parametric kernel regression. We define $p_{n}^{s}:s\in\mathbb{R}\mapsto\displaystyle\frac{\partial\mu^{s}[n\Delta t]}{\partial\lambda^{\otimes 4}}(s)\in\mathbb{R}_{+}$ the associated density. Figure 4 illustrates the evolution of the marginal density of $s^{n}$ with the time.

In figure 4, we can observe the distribution of the sizes is in the beginning above the initial size $s^{0}$ : this is the first stage of the growth in the population, when all plants have their sizes increasing. This corresponds to the densities at times $t=1$ day and $t=5$ days. At some point, the competition becomes too important, mainly at the center of the domain $[0;L]^{2}$ and part of the plants decay, leading to the appearance of plant with sizes lower than $s^{0}$ at time $t=10$ days. Besides, the plants that keep on increasing are the ones that are located in close to the edge $x=L,y=L$, which have faster growth rates $\gamma$ and taller asymptotic sizes $S$. These plants therefore their equilibrium size faster than in the rest domain, so that there are very little change between density at $t=5$ days and density at time $t=10$ days for the plants of size higher than $1.5~{}s^{0}$. This result is consistent with the simulations of the differential system (9) displayed in figure 2.

A clearer visualization of the global behaviour of the MFL distribution can be made by computing the surface corresponding to the averaged size over the domain $[0;L]^{2}$, i.e. the expectation $(x,y)\in[0;L]^{2}\mapsto\hat{e}_{n}(x,y)=\mathbb{E}_{\mu[n\Delta t]}(\hat{s}_{n}(\theta,\omega_{M})|x,y)$, which is obtained by marginalizing growth parameters $\gamma$ and $S$.

[TABLE]

[TABLE]

where $(u_{i})_{1\leq i\leq M}$ and $(u^{\prime}_{i})_{1\leq i\leq M}$ are independent samples from the uniform distribution $\mathcal{U}([0;1])$.

As expected, the surface has its maximal value at the point $x=L,y=L$.The line $x\mapsto\hat{e}_{n}(x,L)$ does not change much from $n=50$ to $n=100$, because plants along this line are already close to their equilibrium size (with competition) at $n=50$, whereas the line $x\mapsto\hat{e}_{n}(x,0)$ has not converged yet for $n=50$, due to its slower growth rate. As $n$ tends towards infinity, we can expect the surface to be more and more invariant by translation along $y-$axis.

5 Conclusion and perspectives

Heterogeneous population models can be approximated by the MFL distribution when the population is large enough and when the interaction function describing the dynamics satisfies a set of assumptions. The phenomenon of chaos propagation, implied by Dobrushin stability, seems to provide an interesting research direction to circumvent the problem of fully-correlated individuals, that arises when the inference of the model is carried out. The suggested methodology for the simulation of the MFL distribution gives promising results, although a theoretical analysis of its consistency still needs to be conducted (on-going work). Our next step is to apply the MFL approximation to real experimental data, to study the impact of competition on the development of plants in mixed stands.

The MFL distribution is appealing because of the property of chaos propagation. But it is clear that MFL approximation might not be relevant for populations having relatively small sizes, as it can be expected when looking at figure 2. A limit seems to be reached for $N>100$, as there are very few changes in the dynamics above this threshold. For smaller $N$ however, the trajectory of a single individual is noisy, and the approximation of the population distribution by a factorized distribution might be too rough. In general terms, the critical size of the population $N_{c}$ is a function of the tolerance $\epsilon$ on some metric quantifying the discrepancy between the microscopic distribution and the MFL distribution, of the length $T$ of the time interval during which the system is observed and finally, of course, of the transition function $g$666An obvious example is given by a transition function $g$ having no dependency with respect to $X^{\prime},\theta^{\prime}$. In that case, the individuals defined by system 4 are in fact independent already, and the critical size is then $N_{c}=1$ for all $\epsilon,T$. The metric, over which a tolerance $\epsilon$ is defined, has to be chosen according to the objectives of the inference. For instance, if our aim is to compute the posterior distribution of some parameter of the model given a set of observations, then we need to find an estimate of the discrepancy between the result that would be obtained by direct inference and the result obtained under MFL approximation. This task seems rather unfeasible, as both of these distributions are either too difficult to compute or simulated with a procedure having a yet uncontrolled error. The upper bound provided by Dobrushin stability in inequality (32) is far too rough to be used for the estimation of the critical size of the population $N_{c}$.

For some systems however, MFL approximation is without doubt relevant. This is the case, for instance, for systems studied in statistical physics, systems that are constituted by a number of particles near or beyond the Avogadro constant ($\approx 6\times 10^{23}$). Even in this favourable case, the use of MFL approximation within a process of Bayesian inference is not well set yet. It would require the coupling of a numerical scheme, similar to the one presented in the previous section, with a time-consuming MH or MHWG algorithms. In machine learning and signal processing literature (Marnissi et al. [2016]), the distribution of the variational Bayes approximation is mainly chosen for its conjugation property with the prior distribution. This often leads to analytical posterior distributions of the parameters, and it may spare a lot of computation time. In our case, our choice is motivated by the behaviour of the dynamical system when it becomes infinite. There is no chance that applying Bayes rule in this context might lead to known or tractable posterior distributions. The solution to this issue may consist in a trade-off between traditional variational Bayes techniques, that are efficient but biased by convenience-motivated choices, and the simulation of the MFL distribution associated to the system, which may require a significant amount of computation time but which is asymptotically unbiased.

This paper has focused on a quite specific class of interaction models, namely the ones that can be decomposed into a sum of pairwise potentials. In the case of more realistic plant models, such decomposition cannot be obtained. The competition is not considered as being additive, and does not even have a closed-form expression in some cases. The necessay complexity of these models leads to the question of the derivation of MFL distribution associated to more generic dynamical systems. In our case, the velocity field $\mathcal{G}_{N}$ at the microscopic level has a linear dependency with respect to the population empirical measure $\mu[0,Z_{N}^{0}]$. The convergence of $\mathcal{G}_{N}$ towards a MFL velocity field may still be preserved when $\mathcal{G}_{N}$ is only continuous with respect to $\mu[0,Z_{N}^{0}]$ for some Wasserstein metric. Such theoretical study, in a more general setting than the one presented in this paper, may enable to study the asymptotic behaviour of more realistic plant population models, incorporating not only competition but also beneficial interactions, which constitute the main interest of mixed cropping.

6 Proofs of the subsection 2.1

6.1 Proof of proposition 1

Proof.

Let us set $\forall i\in\llbracket 1;N\rrbracket,~{}R_{i}=\displaystyle\log\left(\frac{S_{i}}{s_{m}}\right)$. We use the notation $\theta=(x,y,S,\gamma)\in\Theta=\mathbb{R}^{2}\times(\mathbb{R}_{+}^{*})^{2}$. Let us consider the functions

[TABLE]

We set $\theta_{1:N}=(x_{i},y_{i},S_{i},\gamma_{i})_{1\leq i\leq N}$ and the function

[TABLE]

For all $i,j\in\llbracket 1;N\rrbracket$ and $r_{1:N}\in\mathbb{R}^{N}$, we have

[TABLE]

We consider the following norm over $\mathbb{R}^{N}$, defined by $\forall r_{1:N}\in\mathbb{R}^{N},~{}|r_{1:N}|=\displaystyle\sum_{i=1}^{N}|r_{i}|$, known as the norm 1. We have for all $r_{1:N}\in\mathbb{R}^{N}$

[TABLE]

This inequality, along with the fact $G_{r}(\circ,\theta_{1:N})$ is a locally Lipschitz continuous map, proves that the differential system

[TABLE]

has an unique global solution defined over $\mathbb{R}_{+}$. Then the function $t\in\mathbb{R}_{+},~{}s_{1:N}:t\in\mathbb{R}_{+}\mapsto\left(s_{m}e^{r_{i}(t)}\right)_{1\leq i\leq N}$ is the unique solution of system (9). ∎

6.2 Proof of proposition 2

Proof.

We set $\forall t\in\mathbb{R}_{+},~{}Z_{N}(t,Z_{N}^{0})=(s_{i}(t,Z_{N}^{0}),x_{i},y_{i},S_{i},\gamma_{i})_{1\leq i\leq N}$. We consider the random intervall

[TABLE]

The intervall $T[{Z_{N}^{0}}]$ is almost surely an intervall not reduced to singleton $\{0\}$, as $Z_{N}(0,Z_{N}^{0})$ is in $\mathring{\mathcal{D}}^{N}$ almost surely and as $t\in\mathbb{R}_{+}\mapsto Z_{N}(t,Z_{N}^{0})$ is a continuous mapping. $t^{*}(Z_{N}^{0})=\sup(T[{Z_{N}^{0}}])$ is therefore positive random variable almost surely. Let $\Omega^{*}=\{\omega\in\Omega|t^{*}(Z_{N}^{0}(\omega))>0\}$ which is such that $\mathbb{P}(\Omega^{*})=1$. Let $\omega\in\Omega^{*}$. Then for all $t\in[0;t^{*}(Z_{N}^{0}(\omega))[$, we have for all $i\in\llbracket 1;N\rrbracket$

[TABLE]

According to Grönwall lemma, the latest inequalities lead to

[TABLE]

If $t^{*}(Z_{N}^{0}(\omega))<+\infty$, we can use inequalities (45) to obtain that $Z_{N}(t^{*}(Z_{N}^{0}(\omega)),Z_{N}^{0}(\omega))\in\mathring{\mathcal{D}}^{N}$. By continuity and by the fact that $\mathring{\mathcal{D}}^{N}$ is a non-empty open set, we can find $\epsilon(\omega)>0$ such that $Z_{N}(t^{*}(Z_{N}^{0}(\omega))+\epsilon(\omega),Z_{N}^{0}(\omega))\in\mathring{\mathcal{D}}^{N}$, which is in contradiction with the definition of $t^{*}(Z_{N}^{0}(\omega))$. So $\forall\omega\in\Omega^{*},~{}t^{*}(Z_{N}^{0}(\omega))=+\infty$. ∎

7 Proof of the subsection 3.1 : proposition 3

Proof.

We only have to check that the trajectory $t\in\mathbb{R}_{+}\mapsto\mu[t,Z_{N}^{0}]\in\mathcal{P}_{1}(\mathcal{Z})$ is continuous for the metric $W_{1}$. Let $\phi\in\mathcal{C}_{L}(\mathcal{Z})$ a Lipschitz continuous function such that $\mathrm{Lip}(\phi)\leq 1$ and let $t_{1},t_{2}\in\mathbb{R}_{+}$.

[TABLE]

It follows that $t\in\mathbb{R}_{+}\mapsto\mu[t,Z_{N}^{0}]$ is continuous for the metric $W_{1}$, by continuity of the solution of the system (4). The other recquirement for $t\in\mathbb{R}_{+}\mapsto\mu[t,Z_{N}^{0}]$ to be a measure solution is given by equation (17). ∎

8 Proofs of the subsection 3.2

8.1 Proof of theorem 1

Proof.

(of theorem 1) Let us start by proving the local existence of functions satisfying the characteristic flow equation. Let $\alpha>0$. We introduce the following functional space $\mathcal{Y}_{\alpha}=\mathcal{C}^{0}([-\alpha;\alpha]\rightarrow\mathcal{Y})$ endowed with the functional metric $\displaystyle f\in\mathcal{Y}_{\alpha}\mapsto\|f\|_{\mathcal{Y}_{\alpha}}=\sup_{t\in[-\alpha;\alpha]}\|f(t,.)\|_{\mathcal{Y}}$. $\mathcal{Y}_{\alpha}$ is a Banach space for this metric. Over the functional space $\mathcal{Y}_{\alpha}$, we define the following map

[TABLE]

For all $f\in\mathcal{Y}_{\alpha}$, $\Phi_{\alpha}(f,\cdot,\cdot)\in\mathcal{Y}_{\alpha}$. Let $R>1$ and $\mathcal{Y}_{\alpha,R}=\{f\in\mathcal{Y}_{\alpha}|\|f\|_{\mathcal{Y}_{\alpha}}\leq R\}$. We have for all $f\in\mathcal{Y}_{\alpha,R}$, for all $(X,\theta)\in\mathcal{Z}$ and for all $t\in[-\alpha;\alpha]$

[TABLE]

we set $\displaystyle M^{1}_{\mu_{0}}=\int_{\mathcal{Z}}|z^{\prime}|\mu_{0}(\mathrm{d}z^{\prime})$the first order moment, then

[TABLE]

So if we choose $\alpha$ such that $1+K_{1}\alpha(1+(2+M^{1}_{\mu_{0}})R)\leq R$, i.e.

[TABLE]

we have that for all $f\in\mathcal{Y}_{\alpha,R}$, $\Phi_{\alpha}(f,\cdot,\cdot)\in\mathcal{Y}_{\alpha,R}$.

Let $f_{1},f_{2}\in\mathcal{Y}_{\alpha,R}$. We have for all $z=(X,\theta)\in\mathcal{Z}$ and for all $t\in[-\alpha;\alpha]$

[TABLE]

[TABLE]

If $\alpha$ is chosen such that $K_{2}\alpha(2+4R+(1+6R)M_{\mu_{0}}^{1}+2RM_{\mu_{0}}^{2})<1$, and such that it satisfies the inequality (47), i.e.

[TABLE]

then $\Phi_{\alpha}$ is a contractive map over $\mathcal{Y}_{\alpha,R}$. According to fixed-point theorem, there exists an unique $f_{\alpha,R}\in\mathcal{Y}_{\alpha,R}$ such that $\Phi_{\alpha}(f_{\alpha,R},\cdot,\cdot)=f_{\alpha,R}$, i.e. for all $t\in[-\alpha;\alpha]$ and $z=(X,\theta)\in\mathcal{Z}$, we have

[TABLE]

Let us prove now that any function satisfying the equation on a sub-interval of $[-\alpha;\alpha]$ is the restriction of the previous function $f_{\alpha,R}$ to this sub-interval. Without loss of generality, we work on sub-intervals of type $[-\beta;\beta]$ with $\beta\leq\alpha$.

Let $f_{\beta}:[-\beta;\beta]\times\mathcal{Z}\rightarrow\mathcal{X}$ be such that

[TABLE]

We can then distinguish two cases :

Either $\displaystyle\sup_{-\beta\leq t\leq\beta}\sup_{z\in\mathcal{Z}}\frac{|f_{\beta}(t,z)|}{1+|z|}\leq R$. Then, by following the same reasoning as previously, $f_{\beta}$ is the unique fixed point of the map $\Phi_{\beta}$ over the set $\mathcal{Y}_{\beta,R}$. Since the restriction of $f_{\alpha,R}$ to the interval $[-\beta;\beta]$ is also a fixed point of $\Phi_{\beta}$, then we have that $(f_{\alpha,R})_{[-\beta;\beta]}=f_{\beta}$. 2. 2.

Either $\displaystyle\sup_{-\beta\leq t\leq\beta}\sup_{z\in\mathcal{Z}}\frac{|f_{\beta}(t,z)|}{1+|z|}>R$. Let us introduce the following time $\beta_{R}=\sup\{\delta\in[0;\beta]|\forall t\in[-\delta;\delta],~{}\|f_{\beta}(t,.)\|_{\mathcal{Y}}\leq R\}$. Then $\beta_{R}>0$ necessarily, since $\|f_{\beta}(0,.)\|_{\mathcal{Y}}=1<R$. For all $\delta\in[0;\beta_{R}[$, we have, by deriving the same inequalities as in (46),

[TABLE]

By continuity, the previous inequality is also valid for $\delta=\beta_{R}$. Since $\displaystyle\sup_{-\beta\leq t\leq\beta}\sup_{z\in\mathcal{Z}}\frac{|f_{\beta}(t,z)|}{1+|z|}>R$, we have that $\beta_{R}<\beta$. By reinjecting this inequality in (48), we have in fact in that

[TABLE]

which is in contradiction with the definition of $\beta_{R}$. So the current case 2 is absurd.

We can extend the following reasoning to any interval of $\mathbb{R}$ containing 0. We define the following set of tuples

[TABLE]

This set is non-empty as it contains at least $([-\alpha,\alpha],f_{\alpha,R})$ and all its restriction to sub-intervals. $\mathcal{S}_{0,\mu_{0}}$ is partially ordered by the following relationship

[TABLE]

Let us consider the set $\bar{\mathcal{S}}_{0,\mu_{0}}$ of maximal elements of $\mathcal{S}_{0,\mu_{0}}$, i.e.

[TABLE]

We prove now that the set of maximal elements $\bar{\mathcal{S}}_{0,\mu_{0}}$ is reduced to a singleton.

Let $(J_{1},f_{J_{1}}),(J_{2},f_{J_{2}})$ be two maximal elements of $\bar{\mathcal{S}}_{0,\mu_{0}}$. We consider $J=J_{1}\cap J_{2}$ and $T_{+}=\{t\in J|t\geq 0,\forall s\in[0;t],~{}\forall z\in\mathcal{Z},~{}f_{J_{1}}(s,z)=f_{J_{2}}(s,z)\}$. Let us assume by contradiction that $T_{+}\neq J\cap\mathbb{R}_{+}$. If $t^{*}=\sup(T_{+})$, we can exclude two cases :

If $t^{*}=+\infty$, then $T_{+}=\mathbb{R}_{+}$, so $\mathbb{R}_{+}\subset J\cap\mathbb{R}_{+}$, leading to $J=\mathbb{R}_{+}=T_{+}$, which is a contradiction. So $t^{*}$ must be finite, $t^{*}<+\infty$. 2. 2.

If $t^{*}\in\partial J$, i.e. the boundary of interval $J$, then $t^{*}=\sup(J)$, and therefore $T_{+}\cap\mathbb{R}_{+}=J\cap\mathbb{R}_{+}$, which is a contradiction. So under our assumptions, $t^{*}$ must be in the interior of the interval $J$.

For all $t\in[0;t^{*})$, we have $\|f_{J_{1}}(t,.)-f_{J_{2}}(t,.)\|_{\mathcal{Y}}=0$, so by continuity $\|f_{J_{1}}(t^{*},.)-f_{J_{2}}(t^{*},.)\|_{\mathcal{Y}}=0$. Let $\delta>0$ such that $t^{*}+\delta\in J$ and such that $\forall t\in[t^{*};t^{*}+\delta]$, $\max(\|f_{J_{1}}(t,.)\|_{\mathcal{Y}},\|f_{J_{2}}(t,.)\|_{\mathcal{Y}})\leq R^{*}=\|f_{J_{1}}(t^{*},.)\|_{\mathcal{Y}}+1$.

Let $z=(X,\theta)$ and $t\in[t^{*};t^{*}+\delta]$

[TABLE]

The last inequality implies that for all $t\in[t^{*};t^{*}+\delta]$, $\|f_{J_{1}}(t,.)-f_{J_{2}}(t,.)\|_{\mathcal{Y}}=0$ by Grönwall lemma, which is a contradiction with the definition of $t^{*}$. So we have necessarily that $T_{+}=J\cap\mathbb{R}_{+}$. We can conduct the same reasoning to prove that $T_{-}=\{t\in J|t\leq 0,~{}\forall s\in[t;0],~{}f_{J_{1}}(s,.)=f_{J_{2}}(s,.)\}$ is equal to $J\cap\mathbb{R}_{-}$. So the functions $f_{J_{1}}$ and $f_{J_{2}}$ coincide on $J=J_{1}\cap J_{2}$. If $J_{1}\cap J_{2}\subsetneq J_{1}\cup J_{2}$, we could construct the following function

[TABLE]

Then we would have that $(J_{1}\cup J_{2},f_{J_{1}\cup J_{2}})\in\mathcal{S}_{0,\mu_{0}}$ and that $(J_{1},f_{J_{1}})\prec(J_{1}\cup J_{2},f_{J_{1}\cup J_{2}})$ and $(J_{2},f_{J_{2}})\prec(J_{1}\cup J_{2},f_{J_{1}\cup J_{2}})$, which would be in contradiction with the maximality of $(J_{1},f_{J_{1}}),(J_{2},f_{J_{2}})$. So $J_{1}=J_{2}$ and $f_{J_{1}}=f_{J_{2}}$, and $\bar{\mathcal{S}}_{0,\mu_{0}}$ is reduced to a singleton.

Let us prove now that the unique maximal element $(J,f_{J})$ is in fact defined over $\mathbb{R}_{+}$, i.e. $\mathbb{R}_{+}\subset J$. We consider $t^{*}=\sup(J)$. Let us assume by contradiction that $t^{*}<+\infty$. Then we have necessarily that $t^{*}\notin J$. Otherwise, we could apply the same reasoning as for the local existence in the beginning of the proof, to the initial time $t^{*}$ and to the initial distribution $\mu_{t^{*}}$ the probability distribution of $(f_{J}(t^{*},{z_{0}}),{\theta_{0}})$ where ${z_{0}}=({X_{0}},{\theta_{0}})$ is a random variable of distribution $\mu_{0}$. So we would be able to extend the interval of definition $J$, which would be in contradiction with the maximality of $(J,f_{J})$.

Let $t\in[0;t^{*}[$ and $(X,\theta)\in\mathcal{Z}$, we have

[TABLE]

We use the last inequality to show that the derivative $\displaystyle t\in[0;t^{*}[\mapsto\frac{\partial f_{J}}{\partial t}(t,z)$ is bounded for all $z\in\mathcal{Z}$. Let $t\in[0;t^{*}[$, $(X,\theta)\in\mathcal{Z}$

[TABLE]

So $\displaystyle\lim_{t\rightarrow t^{*}}\int_{0}^{t}\left|\frac{\partial f_{J}}{\partial t}(s,z)\right|\mathrm{d}s$ is finite, and $\displaystyle\lim_{t\rightarrow t^{*}}f_{J}(t,z)=X+\lim_{t\rightarrow t^{*}}\int_{0}^{t}\frac{\partial f_{J}}{\partial t}(s,z)\mathrm{d}s$ exists. Then we can define the following function

[TABLE]

We have then that $(J\cup\{t^{*}\},f_{J\cup\{t^{*}\}})\in\mathcal{S}_{0,\mu_{0}}$ and that $(J,f_{J})\prec(J\cup\{t^{*}\},f_{J\cup\{t^{*}\}})$, which is in contradiction with the maximality of $(J,f_{J})$. So $t^{*}=+\infty$ and the maximal element is defined over $\mathbb{R}_{+}$. ∎

8.2 proof of corollary 1

Lemma 2**.**

Let $\mu_{0}\in\mathcal{P}_{2}(\mathcal{Z})$, $g$ satisfying assumptions (A1) and (A2), $z^{0}=(X^{0},\theta)$ a random variable of distribution $\mu_{0}$ and $X_{\infty}:\mathbb{R}_{+}\times\mathcal{Z}\rightarrow\mathcal{X}$ the flow solution of equation (26). For all time $t\in\mathbb{R}_{+}$, we denote by $\mu[t]$ the distribution of the random variable $z^{t}=(X_{\infty}(t,z^{0}),\theta)$. Then $\mu[t]$ is a measure solution of the transport equation (23).

Proof.

Let $\phi\in\mathcal{C}_{L}(\mathcal{Z})$ be a Lipschitz continuous function such that $\mathrm{Lip}(\phi)\leq 1$, and $t_{1},t_{2}\in\mathbb{R}_{+}$.

[TABLE]

The continuity of $t\in\mathbb{R}_{+}\mapsto\mu[t]\in\mathcal{P}_{1}(\mathcal{Z})$ for the metric $W_{1}$ is therefore implied by the continuity of $t\in\mathbb{R}_{+}\mapsto X_{\infty}(t,.)\in\mathcal{Y}$ for the metric $\|.\|_{\mathcal{Y}}$. Let $\varphi$ be a test function. For the initial time $t=0$, $\mu[0]$ is the distribution of $(X_{\infty}(0,z^{0}),\theta)=(X^{0},\theta)$, which is $\mu_{0}$ by definition. So we have

[TABLE]

$t\in\mathbb{R}_{+}\mapsto\displaystyle\int_{\mathcal{Z}}\varphi(t,z)\mu[t](\mathrm{d}z)$ is continuously differentiable and

[TABLE]

∎

Before proving that this measure-solution is in fact the unique one, we need to establish some auxiliary results.

Lemma 3**.**

Let $\mu_{0}\in\mathcal{P}_{2}(\mathcal{Z})$ and $g$ satisfying assumptions (A1) and (A2). Let $t\in\mathbb{R}_{+}\mapsto\nu[t]\in\mathcal{P}_{1}(\mathcal{Z})$ be a trajectory in the space of probability measures continuous for the metric $W_{1}$. Then there exists an unique flow to the (ordinary) differential equation

[TABLE]

where $\mathcal{G}$ is non-local velocity field defined in equation (22).

Proof.

Let $T>0$, $t\in[0;T]$, $z_{1},z_{2}\in\mathcal{Z}$

[TABLE]

It follows that for all $\theta\in\Theta$, the map $(t,X)\in\mathbb{R}_{+}\times\mathcal{X}\mapsto\mathcal{G}(\nu[t],X,\theta)$ is globally Lipschitz continuous over the intervall $[0;T]$, for any $T>0$. The proof is concluded by Cauchy-Lipschitz theorem. ∎

The two following lemmas are classical results from dynamical systems theory and transport equations.

Lemma 4**.**

*(Golse [2013a], theorem 2.2.3)

Let $a:(t,X)\in\mathbb{R}_{+}\times\mathcal{X}\mapsto a(t,x)\in\mathcal{X}$ such that $a\in C(\mathbb{R}_{+}\times\mathcal{X}\rightarrow\mathcal{X})$ and $\displaystyle(t,X)\mapsto\frac{\partial a}{\partial X}(t,X)$ is defined and continuous over $\mathbb{R}_{+}\times\mathcal{X}$. We assume that there exists $K>0$ such that for all $t\in\mathbb{R}_{+}$ and for all $X\in\mathcal{X}$, $|a(t,X)|\leq K(1+|X|)$, and we consider the flow associated to the differential equation*

[TABLE]

Then the flow $X^{a}$ is continuously differentiable with respect to its three arguments, i.e. $X^{a}\in\mathcal{C}^{1}(\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathcal{X}\rightarrow\mathcal{X})$.

Lemma 5**.**

*(Golse [2013a], theorem 2.2.4)

Let $\varphi_{0}\in C^{1}(\mathcal{X}\rightarrow\mathbb{R})$ and $(t,X)\in\mathbb{R}_{+}\times\mathcal{X}\mapsto a(t,X)\in\mathcal{X}$ be such that $a\in C(\mathbb{R}_{+}\times\mathcal{X}\rightarrow\mathcal{X})$ and $\displaystyle\frac{\partial a}{\partial X}\in C(\mathbb{R}_{+}\times\mathcal{X}\rightarrow\mathcal{M}_{d_{\mathcal{X}}\times d_{\mathcal{X}}}(\mathbb{R}))$. We assume that for some $T>0$ there exists $K>0$ such that for all $t\in[0;T]$, $|a(t,X)|\leq K(1+|X|)$. Then there exists an unique solution $\varphi\in C^{1}([0;T]\times\mathcal{X}\rightarrow\mathbb{R})$ to the partial differential equation*

[TABLE]

The solution $\varphi$ has the following expression

[TABLE]

where $(t,t_{0},X)\in[0;T]\times[0;T]\times\mathcal{X}\rightarrow X^{a}(t,t_{0},X)$ is the flow of equation (51).

Lemma 6**.**

Let $\mu_{0}\in\mathcal{P}_{2}(\mathcal{Z})$, g satisfying assumptions (A1), (A2) and (A3), and $(X^{0},\theta)$ a random variable of distribution $\mu_{0}$. Then the unique measure-solution to the transport equation (23) is $t\in\mathbb{R}_{+}\mapsto\mu[t]\in\mathcal{P}_{2}(\mathcal{Z})$ where for all $t\in\mathbb{R}_{+}$, $\mu[t]$ is the probability distribution of ${z}^{t}=(X_{\infty}(t,X^{0},\theta),\theta)$.

Proof.

Let $t\in\mathbb{R}_{+}\mapsto\nu[t]\in\mathcal{P}_{1}(\mathcal{Z})$ be a measure solution to equation (23). Let us consider the flow $X^{\nu}$ associated to the differential equation (50). Thanks to assumption (A3), we have by Leibniz integral rule

[TABLE]

According to lemma 4, for all $\theta\in\Theta$, the map $(t,t_{0},X)\in\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathcal{X}\mapsto X^{\nu}(t,t_{0},X,\theta)$ is continuously differentiable with respect to $t,t_{0}$ and $X$.

Let $\varphi_{0}\in\mathcal{C}^{1}_{0}(\mathcal{Z}\rightarrow\mathbb{R})$, i.e. continuously differentiable and such that $\displaystyle\lim_{|z|\rightarrow+\infty}|\varphi_{0}(z)|+\left|\frac{\partial\varphi_{0}(z)}{\partial z}\right|=0$. We consider the linear transport equation of unknown $\varphi$

[TABLE]

Then, using lemma 5, the unique solution of above equation is

[TABLE]

From previously, we have that $\varphi\in\mathcal{C}^{1}_{0}(\mathbb{R}_{+}\times\mathcal{Z}\rightarrow\mathbb{R})$. As $t\mapsto\nu[t]$ is a measure solution, we can write for all time $t\in\mathbb{R}_{+}$

[TABLE]

If we introduce for all time $t\in\mathbb{R}_{+}$, a random variable $z^{t}=(X^{t},\theta^{t})$ of distribution $\nu[t]$, we can rewrite the above equation as $\mathbb{E}_{\nu[t]}(\varphi(t,X^{t},\theta^{t}))=\mathbb{E}_{\mu_{0}}(\varphi_{0}(X^{0},\theta))$. Let us introduce the random variable $z^{-t}=(X^{-t},\theta^{t})=(X^{\nu}(0,t,X^{t},\theta^{t}),\theta^{t})$ and $\nu_{-t}[t]$ its probability distribution. We have then

[TABLE]

Hence, as the last equality holds for any $\varphi_{0}$ verifying $\displaystyle\lim_{|z|\rightarrow+\infty}|\varphi_{0}(z)|+\left|\frac{\partial\varphi_{0}(z)}{\partial z}\right|=0$, the distributions $\nu_{-t}[t]$ and $\mu_{0}$ are equal for all time $t\in\mathbb{R}_{+}$. $z^{t}$ has the same distribution as the random variable $(X^{\nu}(t,0,X^{0},\theta),\theta)$ and therefore for all $t\in\mathbb{R}_{+}$ and for all $(X,\theta)\in\mathcal{Z}$, we have

[TABLE]

By unicity of the characteristic flow, it follows that $\forall z\in\mathcal{Z},~{}\forall t\in\mathbb{R}_{+},~{}X^{\nu}(t,0,z)=X_{\infty}(t,0,z)$. ∎

9 Proofs of the subsection 3.3

9.1 Proof of lemma 7

Lemma 7**.**

Let $\mu_{0}\in\mathcal{P}_{1}(\mathcal{Z})$ and $g:\mathcal{Z}^{2}\rightarrow\mathcal{X}$ a transition function satisfying assumptions (A1), (A2). For any initial configuration of the population $Z_{N}^{0}\in\mathcal{Z}^{N}$, with $N>1$, there exists an unique function $\hat{X}(Z_{N}^{0},.,.):(t,z)\mapsto\hat{X}(Z_{N}^{0},t,z)\in\mathcal{X}$ such that

[TABLE]

Let $\theta\in\Theta$. We consider the velocity field $(t,X)\in\mathbb{R}_{+}\times\mathcal{X}\mapsto\mathcal{G}(\mu[t,Z_{N}^{0}],X,\theta)\in\mathcal{X}$, where $\mathcal{G}_{N}$ is the non local velocity field defined in equation (18). Let $t\in\mathbb{R}_{+}$ and $X_{1},X_{2}\in\mathcal{X}$.

[TABLE]

So $(t,X)\in\mathbb{R}_{+}\times\mathcal{X}\mapsto\mathcal{G}_{N}(\mu[t,Z_{N}^{0}],X)$ is locally Lipschitz continuous with respect to the variable $X$. Besides, for any time $T>0$ and $X\in\mathcal{X},t\in[0;T]$

[TABLE]

According to Cauchy-Lipschitz theorem, there exists an unique flow for the differential equation of velocity $(t,X)\in\mathbb{R}_{+}\times\mathcal{X}\mapsto\mathcal{G}_{N}(\mu[t,Z_{N}^{0}],X,\theta)$.

[TABLE]

Let $t\in\mathbb{R}_{+}\mapsto(X_{i}(t,Z_{N}^{0}))_{1\leq i\leq N}$ be the trajectories solution of the system (4). For all $i\in\llbracket 1;N\rrbracket$, $X_{i}(0,Z_{N}^{0})=\hat{X}(Z_{N}^{0},0,X_{i}^{0},\theta_{i})=X_{i}^{0}$ and $\forall t\in\mathbb{R}_{+},$ $\displaystyle\frac{\mathrm{d}X_{i}}{\mathrm{d}t}(t,Z_{N}^{0})=\mathcal{G}_{N}(\mu[t,Z_{N}^{0}],X_{i}(t,Z_{N}^{0}),\theta_{i})$. It follows by unicity that $\forall t\in\mathbb{R}_{+},$ $\hat{X}(Z_{N}^{0},t,X_{i}^{0},\theta_{i})=X_{i}(t,Z_{N}^{0})$. We finally obtain that

[TABLE]

9.2 proof of theorem 31

For the simplicity of notations, we introduce the the map $\hat{Z}(Z_{N}^{0},.,.):(t,X,\theta)\in\mathbb{R}_{+}\times\mathcal{Z}\mapsto(\hat{X}(Z_{N}^{0},t,X,\theta),\theta)\in\mathcal{Z}$ and the map $Z:(t,X,\theta)\in\mathcal{Z}\mapsto(X_{\infty}(t,X,\theta),\theta)\in\mathcal{Z}$. Let $\pi_{0}$ be a probability distribution in the set of couplings $\Pi(\mu[0,Z_{N}^{0}],\mu_{0})$ and let $(z_{1},z_{2})$ be a random variable of distribution $\pi_{0}$. We consider for all time $t\in\mathbb{R}_{+}$, the distribution $\pi_{t}$ of the random variable $(\hat{Z}(Z_{N}^{0},t,z_{1}),Z(t,z_{2}))$. Then it is straightforward that $\pi_{t}$ is in the set of couplings $\Pi(\mu[t,Z_{N}^{0}],\mu[t])$. We then have the following inequality.

[TABLE]

Let $z_{1},z_{2}\in\mathcal{Z}$.

[TABLE]

Let us consider the term depending on $A_{N}$.

[TABLE]

As $g$ satisfies assumptions (A2) and (A4), we have for all $z_{1},z_{1}^{\prime},z_{2},z_{2}^{\prime}\in\mathcal{Z}$

[TABLE]

We use also the following notation for all time $t\in\mathbb{R}_{+}$ $\|\hat{Z}(Z_{N}^{0},t,.)\|=\|\hat{X}(Z_{N}^{0},t,.)\|_{\mathcal{Y}}$

[TABLE]

We now look for an upper-bound of the function $t\in\mathbb{R}_{+}\mapsto\|\hat{Z}(Z_{N}^{0},t,.)\|$.

[TABLE]

Now let us consider the empirical characteristic $\hat{Z}(Z_{N}^{0},.,.)$.

[TABLE]

We use the notation $\forall t\in\mathbb{R}_{+},~{}M_{N}(t,Z_{N}^{0})=\displaystyle\frac{(3+2M^{1}(0,Z_{N}^{0}))\exp(2K_{N}(1+M^{1}(0,Z_{N}^{0}))t)-1}{2(1+M^{1}(0,Z_{N}^{0}))}$. Then we obtain that

[TABLE]

Now let us consider the term depending on $B_{N}^{\pi_{0}}$.

[TABLE]

The same reasonning as previously can be applied here to the function $t\in\mathbb{R}_{+}\mapsto\|Z(t,.)\|$ to show that $\forall t\in\mathbb{R}_{+},~{}\|Z(t,.)\|\leq\displaystyle\frac{(3+2M^{1}_{\mu_{0}})\exp(2K_{1}(1+M^{1}_{\mu_{0}})t)-1}{2(1+M^{1}_{\mu_{0}})}=M(t)$, which leads to

[TABLE]

We use the previous inequalities to find an upper-bound of the quantity $D^{\pi_{0}}_{N}(t)=\displaystyle\iint_{\mathcal{Z}^{2}}|\hat{Z}(Z_{N}^{0},t,z_{1})-Z(t,z_{2})|\pi_{0}(\mathrm{d}z_{1},\mathrm{d}z_{2})$.

[TABLE]

[TABLE]

By taking the infimum over $\Pi(\mu[0,Z_{N}^{0}],\mu_{0})$, we obtain the inequality

[TABLE]

Let us study the convergence of the upper bound when the sequence of random variables $(Z_{N}^{0})_{N>1}$ is evaluated at $\omega\in\Omega^{*}$ such that $\lim_{N\rightarrow+\infty}W_{1}(\mu[0,Z_{N}^{0}(\omega)],\mu_{0})=0$. The last convergence implies in particular that $M^{1}(0,Z_{N}^{0}(\omega))\underset{N\rightarrow+\infty}{\overset{}{\longrightarrow}}M^{1}_{\mu_{0}}$. As the distribution $\mu_{0}$ have a compact support of diameter upper bounded by $2R_{0}$, we can write

[TABLE]

So $W_{2}(\mu[0,Z_{N}^{0}(\omega)],\mu_{0})\underset{N\rightarrow+\infty}{\overset{}{\longrightarrow}}0$, and therefore $M^{2}(0,Z_{N}^{0}(\omega))\underset{N\rightarrow+\infty}{\overset{}{\longrightarrow}}M^{2}_{\mu_{0}}$. It follows that $\forall t\in\mathbb{R}_{+}$

[TABLE]

Then we obtain by dominated convergence that $W_{1}(\mu[t,Z_{N}^{0}(\omega)],\mu[t])\underset{N\rightarrow+\infty}{\overset{}{\longrightarrow}}0$.

9.3 Proof of corollary 2

Let us start by establishing an estimation of the quantity $|X_{1}(t,Z_{N}^{0})-X_{\infty}(t,z_{1}^{0})|$ for any time $t$ and for any $Z_{N}^{0}=(z_{1}^{0},...,z_{N}^{0})\in\mathcal{Z}^{N}$.

[TABLE]

where the functions $A_{N}$ and $B^{\pi_{0}}_{N}$ are defined in equation (53) with $\pi_{0}$ any coupling of $\Pi(\mu[0,Z_{N}^{0}],\mu_{0})$. In the proof of theorem 31 (cf appendix section 9.2), we have proved the following inequalities

[TABLE]

We use also as in the previous proof the notation $\displaystyle D^{\pi_{0}}_{N}(t,Z_{N}^{0})=\iint_{\mathcal{Z}^{2}}|z_{1}^{\prime}-z_{2}^{\prime}|\pi_{t}(\mathrm{d}z_{1}^{\prime},\mathrm{d}z_{2}^{\prime})$. The argument of Dobrushin leads to the following inequality

[TABLE]

By gathering the previous inequalities, we obtain finally

[TABLE]

As this inequality holds for any $\pi_{0}\in\Pi(\mu[0,Z_{N}^{0}],\mu_{0})$, we can take $\pi_{0}$ equal to the optimal plan, so that $D^{\pi_{0}}_{N}(t,Z_{N}^{0})=W_{1}(\mu[0,Z_{N}^{0}],\mu_{0})$. By setting $k_{N}(t,Z_{N}^{0})=\displaystyle\frac{e_{N}(t,Z_{N}^{0})}{N-1}+h_{N}(t,Z_{N}^{0})e^{F_{N}(t,Z_{N}^{0})}\left(W_{1}(\mu[0,Z_{N}^{0}],\mu_{0})+\frac{1}{N-1}\int_{0}^{t}E_{N}(s,Z_{N}^{0})e^{-F_{N}(s,Z_{N}^{0})}\mathrm{d}s\right)$, we obtain by Grönwall lemma

[TABLE]

It is clear that for all time $t\in\mathbb{R}_{+}$ and for all $\omega\in\Omega^{*}$, we have that $k_{N}(t,Z_{N}^{0}(\omega))\underset{N\rightarrow\infty}{\overset{}{\longrightarrow}}0$ and $h_{N}(t,Z_{N}^{0}(\omega))\underset{N\rightarrow\infty}{\overset{}{\longrightarrow}}K_{24}(1+2M(t)(2+R_{0}+|z_{1}^{0}(\omega)|))$. By an argument of dominated convergence, we can obtain that

[TABLE]