Simulation and analysis of an individual-based model for graph-structured plant dynamics

We propose a stochastic individual-based model for clonal plant dynamics in continuous time and space, focusing on the effects of the network structure of the plants on the reproductive strategy of ramets. This model is coupled with an explicit advection-diffusion dynamics for resources. After giving a partially exact simulation scheme of the model, the capacity of the model to reproduce specific features of clonal plants, such as their efficiency to forage resources over the field, is numerically studied. Next, we propose a large population approximation of the model for phalanx-type populations, taking the form of an advection-diffusion PDE for population densities, where the influence of the local graph structure of the plant takes the form of a nonlinear dependence in the gradient of resources. Finally, extensions of the model and other possible large population scalings are discussed.


Introduction
Individual-based models are in constant development in computational ecology (Huston et al., 1988;Grimm, 1999;Lomnicki, 1999;DeAngelis and Gross, 1992;Grimm and Railsback, 2005). These models aim to represent the dynamics of populations, and in contrast with conventional models where the population is represented as an aggregate state such as the population size or the total biomass, they explicitly describe each individual as well as each mechanism acting on these individuals. In this sense, conventional models correspond to a macroscopic approach and IBM's to a microscopic one. Most often individuals in IBM's are distributed in a physical (geographical) space or in more abstract spaces like phenotypic trait spaces. Such a description of an ecosystem at the scale of the individual usually rely on stochastic mechanisms reflecting the interactions of the individual with its neighbors and its environment.
We can essentially consider two categories of IBM's. On the one hand, "realistic IBM's" seek to replicate as faithfully as possible the dynamics of an ecosystem, taking into account a large number of mechanisms. These models are difficult or even impossible to analyze. They are used to develop simulation softwares that allow one to highlight emergent properties through simulation.
On the other hand, "simple IBM's" intend to account for the dynamics of simplified ecosystems allowing to focus ✩ This work received support by the French national research. agency (ANR) within the SYSCOMM project ANR-08-SYSC-012. * Corresponding author Email addresses: Fabien.Campillo@inria.fr (F. Campillo), Nicolas.Champagnat@inria.fr (N. Champagnat) on a specific feature of the studied ecosystem. It is possible to simulate such IBM's, but they are also amenable to mathematical analysis. One can for example derive the macroscopic behavior of such systems thanks to scaling up techniques. This last approach is parsimonious as it tries to account for the dynamics of a complex system through a minimum number of mechanisms, or to identify the global effect of a specific macroscopic mechanism at the population level.
One can distinguish several types of IBMs, depending on whether they are defined in continuous or discrete time and space. For clonal plant growth dynamics (e.g. grassland communities), as data are rarely available in continuous time and space, most of the IBM's considered in this domain are in discrete space and time (Dieckmann et al., 1999;Kun and Oborny, 2003;Oborny, 1994;Oborny et al., 2000Oborny et al., , 2001Oborny and Kun, 2002;Winkler and Klotz, 1997;Fischer, 1999, 2002;Winkler and Stöcklin, 2002), see for example the model proposed in (Herben and Suzuki, 2002). Another reason for this choice is that unrealistic spatial over-crowding cannot occur in such models, as only finitely many choices of directions for rhizomes or stolons growth are possible.
In this work, we consider continuous space and time Markovian IBMs, avoiding artifacts due to spatial and temporal grids. Similar models have been introduced in a series of papers for adaptive dynamics in evolutionary biology (Metz et al., 1996;Dieckmann and Law, 1996;Champagnat et al., 2006;Champagnat and Méléard, 2007;Méléard and Tran, 2009;Champagnat and Méléard, 2010) or for competition-colonization in population ecology Pacala, 1997, 1999;Dieckmann et al., 1999 links nodes Figure 1: The plant is represented as a set of nodes connected by links. The nodes can be seen as ramets and the links as rhizomes. Law and Dieckmann, 2000;Fournier and Méléard, 2004;Birch and Young, 2006).
Our goal is to construct and analyze a "simple IBM" putting emphasis on one of the main features of clonal plants, believed to confer them an advantage for resource exploitation and space colonization-their network structure and the transfer of information and resources through rhizomes or stolons. This graph structure allows clonal plants to localize favorable areas in terms of resources (Bell, 1984, foraging), and can favor the growth of new ramets in these areas (Roiloa and Retuerto, 2006). This both allows a better prospection of resources like water, light, organic or inorganic nutrients, and a better survival of descendants. In addition, ramets are able to share resources through links in the network (Stuefer et al., 2004;Wijesinghe and Hutchings, 1997;Hutchings, 1999, physiological integration). In particular, ramets located in areas with few resources can benefit from resources stored by ramets located in richer areas, transmitted through connections between ramets. In addition, translocation between ramets also allows clonal plants to resist to other stresses, like mechanical stress (Yu et al., 2008) or parasitism (D'Hertefeldt and van der Putten, 1998). The model we study here concentrates on these phenomena and their influence on the reproductive strategy. In particular, several other known features of clonal plants are not modeled, e.g. the specialization of ramets (Charpentier and Stuefer, 1999), the resource storage in rhizomes or stolons (Klimeš et al., 1997), or the effects of the graph structure on the growth of ramets and rhizomes or stolons.
The paper is organized as follows. In Section 2, we describe the stochastic individual-based model, coupled with the graph structure of the plant and a piecewise deterministic advection-diffusion model for resources. In Section 3, a partially "exact" Monte Carlo simulation scheme is described. Simulation results are discussed in Section 4. Large population approximations of the population are then given in Section 5. Finally, several extensions of our model are discussed in Section 6.

Dynamics of the plant and of the resources
At time t the clonal plant is represented as a set of nodes (ramets) that may be connected by links (rhizomes  Figure 2: The nodes evolves in a resource landscape (black/white: high/low resource availability). The nodes accessing high levels of resources are more likely to give birth to new nodes ( t 1 → t 2 → t 3 ). Simultaneously the node and the link between this node and the mother node are created. When a node disappears all the links connected with it simultaneously disappear ( t 3 → t 4 ).
or stolons), see Figure 1. The state of the nodes is described by the following finite measure: where x i t ∈ R 2 is position of the ith node and N t total number of nodes; δ x denotes the Dirac measure giving mass 1 to the point x. The measure ν t describes the distribution of nodes over the space D ⊂ R 2 of spatial positions. In this article, we consider for simplicity max ]. The measure ν t is a counting measure: the total mass ν t (D) is the number N t of individuals and ν t (B) is the number of individuals in the subdomain B of D.
For any node at position x we define the set of indices of the nodes connected to x: Note that the basic entities of our IBM are ramets and connections. In particular, we do not account for the growth of ramets and rhizomes or stolons (except through the birth and death mechanisms described below). The plant evolves in a resource landscape (see Figure  3). At each time t, this resource landscape is represented by r(t, x) ∈ [0, r max ] the available resources at position x ∈ D. The nodes accessing high levels of resources r(t, x) are more likely to give birth to new nodes.
First we describe the rates at which nodes disappear and give birth to new individuals, then we describe the dispersion kernel, and finally the dynamics of the resource.

Birth and death rates
Each node of ν t in position x may disappear at a death rate µ(t, x) and give birth to a new node at a birth rate λ(t, x). These rates are per capita rates. Global death and birth rates at population level are respectively: The global event rate is: Basically, the per capita rates depend on the local availability of resources: we suppose that the birth rate λ(t, x) is an increasing function of r(t, x) and the death rate is decreasing function of r(t, x). For example: (4) These rates may account for many other mechanisms. One can limit the number of connections per node with the following birth rate: where |J(t, x)| is the cardinal of the set J(t, x) and N max is the maximum number of connections per node. One can also account for the fact that resources can be translocated through connections. A possibility could be to include a dependence of λ with respect to the incoming flux of resources in the ramet, given by (See also the discussion below.) When a node is added to the population, it is always linked with the mother node, and the set of connections J(t, x i t ) corresponding to the mother node and the new node are modified accordingly. In addition, when a node x is removed from population, all connections to x are suppressed from all the sets J(t, x i t ) (see Figure 2).   (7) is the product of the angle probability density function f and the length probability density function g. When dt,x = 0 then we can for example choose f ( dt,x , θ) = 1 2π (uniform distribution).

Dispersion kernel
A node at position x at time t gives birth to a new node at position y = x + v according to the p.d.f. D t,x (v). We propose the following distribution: where (i) (d t,x , v) is the angle between a preferred direction of reference d t,x and the direction of the new shoot v, f (a, θ) is a p.d.f. on [−π, π) for all a ≥ 0; (ii) g( v ) is the p.d.f. on the length v of the connection (see Figure 4).
For f (a, θ) we can choose the Von Mises distribution with a parameter that may depend on a or, if no direction is preferred, the uniform distribution. In this last case, D t,x (v) ∝ g( v ). Note that, in the case where the preferred direction d t,x is 0, the angle (d t,x , v) is undefined, so one should take the uniform distribution, i.e. f (0, θ) = 1/2π. For g( v ), if we want to fix the length of the connection to ℓ 0 then we just let D t, For the preferred direction of reference d t,x , we need to account for the fact that the ramet can "perceive" the resource map from the connections with other ramets, for example because of resource translocations. One possible choice for d t,x is: This is an approximation of the resource gradient based on the values of r(t, x) at x and at the positions of all connected nodes (see Figure 5). The way d t,x approaches ∇r(t, x) can be made more precise as follows: using the approximation r(t, x i t )−r(t, x) ≈ , valid if the connections are short enough, we always have In addition, in the case where the position of each x i t for i ∈ J(t, x) is uniformly distributed on a circle centered at x (i.e. if f (a, θ) dθ is the uniform distribution on [−π, π)), a simple computation gives In order to keep the population inside the domain D, one should decide what happens at the boundary ∂D of the domain. Two choices are possible: either one considers ∂D as a natural boundary which cannot be crossed by the plants (like rivers, rocks or other areas where plants cannot live), or one thinks of the domain as part of a larger field, but where we neglect the influence of the rest of the field. In both cases, the convenient way to handle births outside of the domain D is simply by ignoring them and the corresponding new connections. In other words, when v is sampled so that y = x+v ∈ D, nothing happens to the population. Note that in resource landscapes like the one of Figure 3, the resource concentration is very small close to ∂D, so that in practice no plant crosses it in simulations.

Interactions between nodes and resources
The natural way to model resource concentration is as a density function r(t, x) over the domain D. Coupling (discrete) individual dynamics with resource density dynamics is a non-standard problem which requires a choice. We propose the following model: with r(0, x) = r 0 (x) and In the absence of plants, the resource concentration r(t, x) follows a classical advection/diffusion transport equation. The last term of (9a) represents the resource consumption of the ramets. We assume a rate of resource intake proportional to the local resource concentration during the whole lifetime of each ramet. In addition, we model with the function Γ the fact that resource consumption is not local. The parameter σ r can then be interpreted as the mean range of roots of a ramet in the species. Of course, other resource consumption models can be considered. For example, consumption can be assumed to occur also at birth, leading the the following resource update at a birth of a ramet at position y and time t: One could also take into account the degradation of plants at position y after its death at time t by updating the resource concentration according to which corresponds to an instantaneous degradation of the plant. More complicated degradation processes could also be considered.
The PDE (9) has to be coupled with appropriate boundary conditions. We assume here that the boundary of the domain D is a natural boundary that cannot be crossed by the plants, because of the absence of resources. This corresponds to the Dirichlet boundary condition If one thinks of D as part of a larger field, one would have to consider flux (Neumann) boundary conditions by giving a value to the gradient of r at x ∈ ∂D in the direction normal to ∂D at x. Local increase or decrease of resources (due for example to rain) can also be easily taken into account by adding a source term f (t, x) to the right-hand side of the PDE (9a).

Numerical approximation of the IBM
We now describe the simulation algorithm: starting from the state at last event time T k−1 , we first sample the time of the next event (birth or death): andλ T k−1 andμ T k−1 defined by (3a). The next event: • is a birth event with probabilitȳ Then sampleî according to and update the sets of connections J(T k , x) accordingly; • is a death event with probabilitȳ Then sampleî according to and update the sets of connections J(T k , x) accordingly.
Note that this algorithm is valid if the rates λ(t, x i t ) and If it is not the case, we should make use of an acceptance/rejection algorithm (Fournier and Méléard, 2004;Campillo and Joannides, 2009). In parallel, we should numerically integrate the PDE (9), for example with implicit or explicit finite-difference schemes. In practice, if T k − T k−1 is small enough, which is usually the case if the population is large, a single time step of the finitedifference scheme is sufficient. Note also that, in order to compute the birth and death rates λ and µ, one needs to interpolate the resource concentration at each ramet position from the resource concentrations on the discretization grid. The algorithm is presented in Algorithm 1.
Note that this algorithm can be very costly when the population size is large, as it requires to compute the sums λ andμ at each time step. Another possibility is to use an acceptance/rejection procedure (Fournier and Méléard, 2004) in order to replace this sum by a random sampling of the individual to which the next event will apply. The drawback is that some (and sometimes many) of these events may actually be void, leading to an increase of the number of time steps in the simulation.

Simulation
The proposed model can account for conventional phalanx or guerrilla dynamics. More specifically, in (7): (i) if the p.d.f. f (a, θ) of the shoot angle has a small variance and if the p.d.f. g( v ) favors large lengths, then the model will present the characteristics of a guerrilla plant; (see Figure 6).
In Figure 7 we focus on the local plant network in the guerrilla and in the phalanx cases. As there is no limitation in terms of number of connexions, the plant network does not appear very realistic. The architecture of the plant can be improved for example by limiting the number of connexions per node to N max = 3 (see (5)). Simulation results are given in Figure 9 and show a more realistic plant network.
One of the main advantages that network structure are believed to confer to clonal plants is their capacity to adapt their reproductive strategy to the local resources landscape, felt either by foraging or physiological integration. This capacity is given in our IBM by the dependency of the dispersion kernel-more precisely the shoot angle p.d.f. f -on the estimated resource gradient d t,x (8). The efficiency to exploit resources can then be compared with simulations of IBM involving different dispersion kernels. This is done in Figure 8 with the intial resource landscape given by Figure 3. In this landscape, resources are distributed in a fragmented way as several patches. In such a landscape, an efficient strategy requires (1) an ability to cross areas with small resources densities in order to visit new patches and (2) to be able to absorb resources within a patch efficiently. Pure guerrilla strategies are efficient for the first point (because the plant grows along straight lines), but inefficient for the second point. Pure Phalanx strategies are less efficient for the first point, but more efficient for the second point (as they grow in all directions). small variance) and if the p.d.f. g( v ) favors small lengths (resp. large lengths), then the model will present the characteristics of a phalanx growth strategy (resp. guerrilla growth strategy).
We explore a range of intermediate strategies simply by varying the variance of the shoot angle distribution f . A small variance means that the plant grows nearly along straight lines and has a very small probability to change its growth direction. A large variance corresponds to the case where the graph structure has no influence on the population dynamics, i.e. the case where f is the uniform distribution on [0, 2π). As expected, the simulations of Figure 8 show that there is an optimal tradeoff: when the shoot angle p.d.f. is not directive enough or too directive, the plant fails to reach areas with high level of resources; in intermediate cases, the plant reaches these areas and reaches them rapidly.

Large population approximation of the IBM
Individual-based models are a convenient tool for modeling small-scale ecological systems. However, their simu- guerilla phalanx Figure 7: Typical graph patterns for the guerrilla and the phalanx cases. In Figure 9 we present an example where the number of connexions per node is limited.
lation is often costly and can hardly provide relevant fieldscale informations in reasonable computational time. This is typically the case for phalanx-type clonal plants, where the connection length is often short and the plant architecture is dense. In order to understand the global interaction between plants and resources, the search for simpler approximation models is crucial. One then typically seeks for partial differential equations (PDE) governing the time dynamics of the population density over space, obtained in a limit of large population. This has been done under various scalings of the individual parameters for plants systems without network structure (Fournier and Méléard, 2004;Champagnat et al., 2006Champagnat et al., , 2008. Three main families of scalings were described in the second paper: in the first scaling, space is unscaled, leading to a non-local integro-differential equation for the population density; in the second one, space is scaled, and births and deaths are accelerated accordingly in order to obtain a limit, which takes the form of a local reactiondiffusion PDE; in the last scaling, births and deaths are even more accelerated, leading to a stochastic reactiondiffusion PDE. Because of the underlying network structure and the explicit coupling with resource dynamics, these results do not apply directly to our IBM. Similar results can be obtained in network-structured populations (Campillo and Champagnat, 2010), but without explicit dynamics for the resources. The result we present here is a first attempt to fill this gap. A full mathematical proof is likely to be very technical because our model couples a stochastic, discrete structure for the population and a deterministic, continuous structure for resource concentrations. We leave this proof as an open problem and only present and justify our macroscopic model as a reasonable large population approximation of the IBM.

The limit PDE
The more relevant scaling within the context of phalanxtype clonal plants is the space-scaling and acceleration of births and deaths which leads to a reaction-diffusion PDE for population densities. It is also the less technical, and  Figure 8: We assume that for all a > 0, f (a, θ) = f (θ) is the Von Mises distribution with variance parameter κ. Taking different values for this parameter (κ = 0, 0.1, 0.5, 2, 3) we plot, on the left, the corresponding shoot angle p.d.f. f (θ) and, on the right, the time evolution of the size of the population (red) and of the total resource (blue); we suppose that there is no advection of resources (i.e. b ≡ 0 in (9)). The simulation is done with the initial resource map of Figure 3 with high resources spot on the north of the domain, a less important spot in the south and 3 negligible spots (left, center and right), the initial plant is located on the left spot. In the first case, κ = 0, the plant explore all directions without any preference (uniform distribution of the shoot angle): by chance the plant reaches the north spot corresponding to the increase of population and the population subsequently decreases. In the second case, κ = 0.1, the plant first reaches the south spot and then the north spot. In the third case, κ = 0.5, the plant rapidly reaches the north spot and then the south spot. In the fourth case, κ = 2, the plant need more time to reach the two important spots. In the last case, the plant does not reach any resource spot and node population goes extinct rapidly. Guerrilla Phalanx n Figure 9: Here the birth rate is adjusted so that each node can have a maximum of 3 connections (i.e. Nmax = 3 in (5)).
we restrict ourselves to this case here. Concerning other possible scalings, we refer to the discussion below (see also Campillo and Champagnat. 2010). In this case, the PDE approximation of the IBM takes the following form: denoting by u(t, x) the population density at time t and position x in D, (11) This model is a diffusion-advection PDE with a term of local growth for the population density coupled with an diffusion-advection PDE for resources, where the resource consumption depends on the local population density u(t, x) and the advection of the population density depends on the local gradient of resources ∇r(t, x). Note that the PDE for u is linear, but depends non-linearly on r through F , λ and µ. Therefore, this is not strictly speaking a reaction diffusion PDE as in the competition-colonization plant model of Champagnat et al. (2006). Note also that the local growth rate of the population density is given as one could expect by λ(t, x) − µ(t, x).
The parameters β, γ, δ and F appearing in this equation are obtained as explicit functions of individual parameters. Before giving their expression, we first explain the scaling corresponding to this PDE and how it is obtained.

The parameter scalings
Fix a parameter η ∈]0, 1[ and introduce a constant K which will correspond to the order of the population size. This constant K is related to the carrying capacity of the system and is also often called "system size" (Metz et al., 1996). We denote by ν K and r K the population measure and resource concentration corresponding to the parameter K, and we scale ν K as First, in order to make the population larger, we need to increase the quantity of resources available, or equivalently to reduce the rate of resource consumption per individual: which amounts to replacing α by α/K in (9). We also apply a space scaling, which corresponds to the choice in the resource dynamics. This means that the typical size of a ramet (and its roots) is of order K −η/2 . Since the integral of Γ K y over space remains of order 1, the total amount of resources consumed per ramet and per unit of time remains of order 1/K (because α has been divided by K). This is consistent with a total population size of order K.
Concerning the population dynamics, we require that births and deaths are accelerated as follows: where λ and µ may depend on r as in (4). This means that, while keeping the growth rate per individual λ K − µ K = λ − µ constant, births and deaths are accelerated by the term K η γ, where γ is a function that may depend on x ∈ D. This corresponds to a time-scale where ramets have short life spans and reproduce rapidly, while the total population grows or declines on a slower time scale. In other words, we want to describe the global changes in the population on an intermediate time scale, after the initial (fast) invasion of space, but before it stabilizes to a global equilibrium. Finally, the space scaling must also be applied to the dispersion kernel, and the preferred dispersion direction must be consistent with the phalanx-type population we consider: we replace the dispersion distribution of (7) by The scaling of the p.d.f. g corresponds to the space scaling of K −η/2 which was used above for the resource consumption kernel Γ K . The scaling of the p.d.f. f amounts to a uniform dispersion direction with probability 1 − K −η/2 and otherwise to a dispersion direction around the preferred direction according to f . This corresponds to phalanxtype populations, and we will see below that this scaling is the one that leads to a non-trivial limit when K goes to infinity. Note that the PDE approximation (11) shows that even a slight effect of the preferred direction choice can lead to a non-trivial global effect at the field scale (expressed through the term F (x, ∇r(t, x))). We refer to the discussion below for other possible scalings of the dispersion kernel.

Approximation argument
The way Equation (11) can be derived from the previous scaling is as follows: we assume that the initial population state ν K 0 converges to some non-trivial measure u 0 (x)dx on D, i.e. that the initial population size is of order K and is distributed over space close to the population density u 0 (x). Now, assume that K is large and that at some time t, the population state ν K t is close to the density u(t, x), and consider some test function φ over D. Then, the quantity A t = D φ(x) ν K t (dx) changes on the time interval [t, t + dt] due to the following three effects: • Each individual close to x ∈ D gives birth to a number λ K (t, x) dt of individuals, and a proportion µ K (t, x)dt die. This gives a local change of the population density (λ K (t, x) − µ K (t, x)) u(t, x) dt, and a change of A t given by This gives the first-order change of A t . Higher-order changes (i.e. changes involving derivatives of φ) must also be taken into account.
• Among all the births from the individuals around x ∈ D, a proportion 1 − K −η/2 are located symmetrically around x at a distance K −η/2 ρ, where ρ is distributed as g(ρ). Making a Taylor expansion of φ around x, the corresponding change of A t is given by where v θ = (cos θ, sin θ) is the unit vector of direction θ. Note that the first-order term does not appear as it vanishes because of the symmetry of dispersions. Note also that the scaling of λ K in (13) has been chosen so that it exactly cancels the scaling of the dispersion distance K −η/2 ρ. Neglecting the smaller order terms, an easy computation gives a change of • Among all the births from the individuals around x ∈ D, a proportion K −η/2 are located at where v has distribution (7). Making again a Taylor expansion of φ, we obtain a local change of A t of where v t,x is the mean of v t,x when v has distribution (7), in which d t,x has the distribution of the preferred dispersion direction in the population around x. Of course, this last distribution is unknown and must be computed (see below).
Putting together all the previous results, we obtain This is the weak formulation of the PDE (11), where we put The fact that v t,x only depends on x and ∇r(t, x) will be justified below.
Before coming to the computation of F , we observe that the last term of the PDE (12) for r K may be written as Replacing ν K (dy) with its limit u(t, y)dy and passing to the limit K → +∞ would formally give the second equation of (11) with δ = α √ 2π σ r .
However, this requires a sufficient regularity for u(t, x) and a convergence of ν K t (dy) to u(t, y)dy such that space means of ν K t over balls of diameter of the order of K −η/2 also converge to u(t, x). This is the place where a full mathematical proof requires a careful analysis that we do not examine here.
Finally, we need to compute the distribution of the preferred dispersion direction over the population close to a fixed x ∈ D at a fixed time t ≥ 0. We consider one of these ramets chosen at random and call it RSI (Randomly Sampled Individual) at x, and we divide the computation into several steps.
Scaled age distribution of a RSI. Because of the scaling of the death rate µ K (t, x), the life span of the RSI is of order K −η , and the position of its mother ramet is at a distance of order K −η/2 because of the scaling of the p.d.f. g. Over these time and space scales, the population densities and resource concentration can be assumed constant equal to u(t, x) and r(t, x). We call scaled age of an individual at x at time t, its age multiplied by K η .
In particular, the ramets alive at time t close to x appeared at an (approximately) constant rate during a time length of order K −η . Since they die (approximately) at constant rate K η γ(x), the probability of survival of an individual born at time t − yK −η is e −yγ(x) . Since the birth rate of individuals close to x is approximately constant, we deduce that the scaled age of a RSI at x is asymptotically distributed when K → +∞ as an exponential random variable with parameter γ(x).
Number of living offspring of a RSI. Conditionally on the age K −η y (or the scaled age y) of a RSI at x, the birth rate at x over the time interval [t − K −η y, t] is constant (approximately) equal to K η γ(x). Therefore, the scaled age of all its offspring is distributed as a Poisson process of parameter γ(x) on the interval [0, y]. In particular, the total number of offspring is distributed as a Poisson random variable with parameter yγ(x). Among them, some are already dead at time t.
In order to compute the number of living neighbors of the RSI, we need to introduce an additional conditioning on the number k of offspring of a RSI conditioned to have scaled age y. We use the following feature of Poisson processes: conditionally of the number k of offspring of the RSI, their scaled age are distributed as k independent and identically distributed (i.i.d.) random variables with uniform distribution on [0, y]. In other words, conditionally on the age K −η y of the RSI and the number of its children k, the ages of its children are distributed as k i.i.d. random variable with uniform distribution on [0, K −η y].
Since they all die at (approximate) rate K η γ(x), the survival probability of each is We denote by p x,y the right-hand side of this equation. This leads to a number of living children of the RSI with binomial distribution with parameters k and p x,y .

Number of living neighbors of a RSI.
Similarly, conditionally on the scaled age y of the RSI, its mother ramet is still alive at time t with probability e −yγ(x) . Therefore, the number N of living neighbors of a RSI at x conditionally on A = y, where A is the scaled age of the RSI, is distributed as Denoting for convenience these quantities q x,y (n), we obtain the following distribution for the number N of living neighbors of a RSI at x: for all n ≥ 0, Positions of the neighbors of a RSI. Neglecting the smallerorder terms in (14), we have , which means that the direction of dispersal is (approximately) uniform for each birth. This is also the case for the births from the RSI at x, as well as its mother ramet. Therefore, the positions of all the neighbors of the RSI at x are i.i.d. and distributed as x + K −η/2 ρv θ , where ρ has p.d.f. g(ρ) and v θ = (cos θ, sin θ).
Asymptotic distribution of d t,x . Since the length of connections is of order K −η/2 , we can replace r(t, (8) for the preferred dispersion direction d t,x . Therefore, conditionally on N = n ≥ 1, the preferred dispersion direction d t,x of a RSI at x is distributed as where θ 1 , θ 2 , . . . are i.i.d. random variable uniformly distributed on [0, 2π). Introducing θ t,x such that ∇r(t, x) = ∇r(t, x) v θt,x , and applying the rotation of angle −θ t,x to the vector (16) where the θ ′ i are i.i.d. and uniformly distributed on [0, 2π). Denoting by µ(w)dw the distribution of cos(θ)v θ for θ uniformly distributed on [0, 2π), the distribution of the vector (17) is µ * n (nw) dw, where µ * n = µ * . . . * µ and * denotes the convolution operator. Therefore, a RSI at x at time t has the following distribution for its preferred dispersal direction d t,x rotated of the angle −θ t,x : where P(N = n) is given by (15). Denoting by π(dw) the previous distribution, we finally obtain Note that the right-hand side only depends on x through the distribution of N , and on ∇r(t, x) through the angle θ t,x .

Extensions and perspectives
We presented in this article an individual-based model for clonal plants which focuses on the network structure and its effects on the reproductive strategies of ramets. In order to specifically study this network structure and make possible the large population analysis of the model, we concentrated on foraging (Bell, 1984) and physiological integration (Stuefer et al., 2004), and ignored other features of clonal plants like specialization of ramets (Charpentier and Stuefer, 1999), resource storage (Klimeš et al., 1997) or the use of network structure to handle other stresses, e.g. mechanical or parasitic (Yu et al., 2008;D'Hertefeldt and van der Putten, 1998). In addition, we focused on the reproductive strategies and ignored effects of the network structure on the growth of ramets and spacers. In contrast to Campillo and Champagnat (2010), we coupled the IBM with an explicit resource dynamics in order to describe explicitly the effects of foraging and physiological integration of resources on the population dynamics.
We presented a partly explicit 1 simulation scheme, and some numerical simulations of the invasion of clonal plants in a given resource landscape. In particular we were able to reproduce the phalanx and guerirlla patterns specific of various clonal plant species. We also compared the efficiency of different invasion strategies of clonal plants in a particular example of a fragmented resource landscape.
Next, we proposed a PDE model which approaches the IBM when the population is large, the births and deaths are fast and the size of ramets and spaces is small. This particular scaling is relevant for clonal plants with guerrilla strategies. The justification of the PDE requires us to compute the local distribution of the network structure around a randomly sampled individual.
In order to keep the explanation short, we described a simple IBM, however our model can easily be extended to more realistic situations, which can be analyzed either numerically like in Section 4, or mathematically like in Section 5. Let us discuss in details how to extend our model and results in several directions.
• We describe the possible effect of the network structure on the choice of dispersion direction (foraging) by introducing the "estimated gradient" (8). The model and the large population results can be easily extended to the case where one takes into account physiological integration by assuming a dependance of the birth rate λ(t, x) on the "resource flux" felt the ramet at x, due to resource translocation via connections. This flux is given by the formula (6).
In particular, in a scaling of large population like the one studied in Section 5, the computation of the limit PDE can be easily adapted to this case, using our computation of the distribution of the number and positions of neighbors of a ramet.
• It is also possible to introduce more complex dependancies of the birth and death rates λ(t, x) and µ(t, x) on the graph, such as the one proposed in (5). Such extensions are desirable as they produce more realistic graph architectures for the plant (compare Figures 7 and 9). The simulation of the corresponding IBM would be as easy as in Section 3, but the computation of the limit PDE would be more complex. In principle, such a computation should be possible (although maybe not explicit) as long as the birth and death parameters of a ramet only depend on its nearest neighbors in the network.
• Our IBM is defined in continuous time and space, in contrast with classical IBMs of clonal plants (Herben and Suzuki, 2002, e.g.). This avoids artifacts due to the discrete grid structure, but can lead to unrealistic plant shapes, since a ramet can produce two very close new ramets. This can be avoided by assuming a dependency of the direction where an individual can produce a new ramet, on the positions of all living neighbors of this individual. Again, the simulation of the IBM would be easy, but here the mathematical analysis is likely to be much more complicated.
• Since connections are explicit in our IBM, it is also possible to include more complicated dynamics for spacers. For example, they could disconnect at some given rate which could depend on the spatial position. One could also model resource storage in spacers by specifying the amount of resources contained by each connection. These resources could then be used by neighboring ramets to reproduce with rates similar to λ(t, x). As far as we know, the mathematical analysis of such models is open.
• Concerning resources, advection-diffusion PDEs allow one to easily take into account flux of resources at the boundary of the domain D (Neumann boundary conditions), local increase or decrease of resource concentrations (by adding a source term to the PDE), or local loss or degradation of resources by the soil (by adding a negative zeroth-order term to the PDE).
The main modeling difficulty is to determine how the discrete individuals should interact with the continuous resource concentrations. In Section 2, we propose simple examples easy to analyze, but a more precise understanding of the range and rate of resource absorption during the lifetime of a plant is necessary to construct more realistic IBMs.
Although simple, our IBM allows one to recover and describe various specific features of clonal plants. The simulations of Section 4 show that guerrilla or phalanx-type dynamics are very easy to obtain. They can be made even more realistic by slightly changing the birth rates (Figure 9).
In addition, our IBM reproduces known facts about the advantage conferred by the reproductive strategies for the invasion of space and the consumption of resources. The simulations of Section 4 illustrate this fact in a specific resource landscape, where the optimal strategy appears to be intermediate between pure-guerrilla strategies (corresponding to clonal plants without network structure) and pure-phalanx strategies. The phalanx strategy allows the plant to cross low resource areas and to visit new patches of resources. The guerrilla strategy allows the plant to exploit the resources of a patch more efficiently.
The mathematical analysis of the IBM proposed in Section 5 is only done in a very specific case (large population for guerrilla strategies with time and space scaling). We actually restricted ourselves to the simplest possible scaling, but even in this case, the mathematical analysis is non-trivial, and the argument presented in Section 5 is not a complete proof, as several mathematical difficulties remain unsolved: • First, in order to prove the convergence of the resources, we need to pass to the limit K → +∞ in (12). But this requires that the population state ν K t converges to the solution u(t, x) of (11) in a stronger sense than the usual weak convergence of measures (see Section 5).
• Second, the limit PDE we obtain takes into account the local network structure of the plant through a dependency of the advection of population density on the gradient of resources. This kind of nonlinearity makes the analysis of the PDE much more complex. In particular, one needs to check the wellposedness of the PDE (11) to prove the convergence, and for practical issues, the convergence of approximation schemes of the solution should be analyzed. The non-linearity in the gradient of r makes these questions non-classical.
Several other scalings of our IBM can be considered, in the fashion of those obtained in clonal plants IBMs without explicit resource dynamics by Campillo and Champagnat (2010).
• The acceleration of births and deaths studied in Section 5 implicitly assumes that the global demographic changes in the population occur on a slower timescale than the individual births and deaths. This is the correct way to obtain a diffusion-advection PDE for the population density u(t, x). It is also possible to make a large population scaling without space scaling and birth and death acceleration (i.e. with the parameter η of Section 5 set to 0). As a result, the global changes in the population occur at the same time scale as individual births and deaths, and the resulting PDE is non-local as the network is unscaled. In particular, the effect of the preferred dispersion direction on the dynamics is much harder to compute. Still, it is possible to write a limit PDE similar to the one of Campillo and Champagnat (2010), by taking into account the age and parent position distribution of RSI at x. As the computation of F in Section 5 suggests, this information is crucial to characterize the distribution of the estimated gradient d t,x , which is itself needed to compute the distribution of the position of a new ramet born from a RSI at x. The limit PDE takes the form of an integro-differential equation with age and parent position structure. Note that this scaling applies both to guerrilla and phalanx IBMs.
• The space scaling and births and deaths acceleration of Section 5 can be also applied to phalanx-type strategies by taking as dispersal distribution Comparing with (14), we see that the direction of all dispersions is given by the p.d.f. f , and that space is scaled by K −η instead of K −η/2 . Assuming that the space is also scaled by K −η in (12), one can easily follow the argument of Section 5 to obtain a limit PDE for the population density without diffusion. This is a pure advection (or transport) PDE with non-linearity in ∇r. Therefore, when space is scaled and births and deaths are accelerated, the PDE models corresponding to phalanx and guerrilla strategies are drastically different. In the first case, diffusion and advection occur on the same time scale; in the second case, the dominant phenomenon is advection, and diffusion can only occur on a much longer time scale.
• Finally, the analysis of Section 5 can also be done in the case η = 1, corresponding to even faster births and deaths. In this case, some randomness remains in the limit of large population, in the fashion of the Fleming-Viot process. The limit process is a superprocess (a continuous measure-valued stochastic process), or a stochastic PDE (Etheridge, 2004;Champagnat et al., 2006). Again, the non-linearity in ∇r is going to lead to technical difficulties in the mathematical analysis of this convergence, as the theory of stochastic quasilinear PDEs is difficult and poorly developed today.