Vincent Bansaye
email: vincent.bansaye@upmc.fr

Cell contamination and branching process in random environment with immigration

Keywords: Branching processes in random environment with immigration (IBPRE). Markov chain indexed by a tree. Empirical measures. Renewal theorem 60J85, 60K37, 92C37, 92D25, 92D30

We consider a branching model for a population of dividing cells infected by parasites. Each cell receives parasites by inheritance from its mother cell and independent contamination from outside the cell population. Parasites multiply randomly inside the cell and are shared randomly between the two daughter cells when the cell divides. The law of the number of parasites which contaminate a given cell depends only on whether the cell is already infected or not. We determine first the asymptotic behavior of branching processes in random environment with state dependent immigration, which gives the convergence in distribution of the number of parasites in a cell line. We then derive a law of large numbers for the asymptotic proportions of cells with a given number of parasites. The main tools are branching processes in random environment and laws of large numbers for Markov tree.

Introduction

We consider the following model for cell division with parasite infection and state dependent contamination. The cell population starts from one single cell and divides in discrete time. At each generation, (i) the parasites multiply randomly inside the cells, (ii) each cell is contaminated by a random number of parasites which come from outside the cell population, (iii) each cell divides into two daughter cells and the parasites are shared randomly into the two daughter cells.

It is convenient to distinguish a first daughter cell called 0 and a second one called 1. We denote by T = ∪ n∈N {0, 1} n the binary genealogical tree of the cell population, by G n the set of cells in generation n and by Z i the number of parasites of cell i ∈ T. We write then i0 and i1 the two daughter cells of the cell i ∈ T.

First, we describe by a branching process the random multiplication and sharing of parasites in the cell, i.e. this branching process combines (i) and (iii). Second, we describe the random contamination (ii) by immigration. Finally, we combine both in an i.i.d. manner to fully describe the model. I Parasite infection and cell division For every cell, we choose randomly a mechanism for multiplication of the parasites inside and sharing of their offspring when the cell divides. This mechanism is independent and identically distributed for every cell. Its distribution is specified by a random couple probability generating function (p.g.f) f . This means that f is a.s. the p.g.f of a pair of random variables taking values in N.

More precisely let (f i ) i∈T be a sequence of i.i.d. couple p.g.f distributed as f . For each cell i, f i gives the reproduction law and sharing of the offspring of its parasites in the following way. For every i ∈ T, let (X

(0) k (i), X (1)
k (i)) k∈N be a sequence of r.v. such that conditionally on f i = g, (X

(0) k (i), X (1)
k (i)) k∈N are i.i.d. with common couple p.g.f g:

∀i ∈ T, ∀k ∈ N, ∀s, t ∈ [0, 1], E(s X (0) k (i) t X (1) k (i) | f i = g) = g(s, t).
Then, in each generation, each parasite k of the cell i gives birth to X (0)

k (i) + X (1)
k (i) children, X (0) k (i) of which go into the first daughter cell and X [START_REF] Asmussen | Branching processes[END_REF] k (i) of which go into the second one, when the cell divides. This is a more general model for parasite infection and cell division than the model studied in [START_REF] Bansaye | Proliferating parasites in dividing cells : Kimmel's branching model revisited[END_REF], where there was no random environment (f was deterministic) and the the total number of parasites was a Galton Watson process. See [START_REF] Kimmel | Quasistationarity in a branching model of division-withindivision[END_REF] for the original model in continuous time.

Our model includes also the two following natural models, with random binomial repartition of parasites. Let Z be a random variable in N and (P i ) i∈T be i.i.d. random variable in [0, 1]. In each generation, every parasite multiplies independently with the same reproduction law Z. Thus parasites follow a Galton Watson process. Moreover P i gives the mean fraction of parasites of the cell i which goes into the first daughter cell when the cell divides. More precisely, conditionally on P i = p, every parasite of the mother cell i chooses independently the first daughter cell with probability p (and the second one with probability 1 -p). It contains also the following model. Every parasite gives birth independently to a random cluster of parasites of size Z and conditionally on P i = p, every cluster of parasite goes independently into the first cell with probability p (and into the second one with probability 1 -p).

We want to take into account asymmetric repartition of parasites and do not make any assumption about f . Indeed unequal sharing have been observed when the cell divides, see e.g. experiments of M. de Paepe, G. Paul and F. Taddei at TaMaRa's Laboratory (Hpital Necker, Paris) who have infected the bacteria E. Coli with a lysogen bacteriophage M13 [START_REF] Stewart | Aging and Death in a Organism that reproduces by Morphologically Symmetric Division[END_REF]. In Section 6.1, we consider this model where a cell receive parasites only by inheritance from its mother cell. We determine when the number of infected cells becomes negligible compared to the number of cells when the generation tends to infinity. II State dependent contamination In each generation, each cell may be contaminated by a random number of parasites which also multiply randomly and are shared randomly between the two daughter cells. This contamination depends only on whether the cell already contains parasites or not.

More formally, if a cell i contains x parasites, the contamination brings Y (0) x parasites to the first daughter cell of i and Y [START_REF] Asmussen | Branching processes[END_REF] x to the second one, where

∀x ≥ 1, Y 1 : d = Y (0) x d = Y (1)
x , Y 0 :

d = Y (0) 0 d = Y (1)
0 .

Moreover we assume that contamination satisfies

0 < P(Y 0 = 0) < 1, 0 < P(Y 1 = 0), (1)
which means that each non-infected cell may be contaminated with a positive probability but the cells are not contaminated with probability one.

This model contains the case when the contamination is independent of the number of parasites in the cell (Y 0 and Y 1 are identically distributed). It also takes into account the case when only non infected cells can be contaminated (Y 1 = 0 a.s.) and the case when infected cells are 'weaker' and parasites contaminate them easier (Y 1 ≥ Y 0 a.s.). For biological and technical reasons, we dot make Y x depend on x ≥ 1. But the results given here could be generalized to the case when the contamination depends on the number of parasites x inside the cells soon as x is less than some fixed constant.

III Cell division with parasite infection and contamination

We describe now the whole model. We start with a single cell with k parasites and denote by P k the associated probability. Unless otherwise specified, we assume k = 0.

For every cell i ∈ T, conditionally on Z i = x and f i = g, the numbers of parasites (Z i0 , Z i1 ) of its two daughter cells is distributed as

x k=1 (X (0) k (i), X (1)
k (i)) + (Y (0)
x (i), Y (1) x (i)), where

(i) (X (0) k (i), X (1)
k (i)) k≥1 is an i.i.d. sequence with common couple p.g.f g.

(ii) (Y (0)

x (i), Y (1)
x (i)) is independent of (X

(0) k (i), X (1)
k (i)) k≥1 .
Moreover, (X

(0) k (i), X (1)
k (i)) k≥1 , (Y (0)
x (i), Y

x (i)) x≥0 are i.i.d. for i ∈ T. This model is a Markov chain indexed by a tree. This subject has been studied in the literature (see e.g. [START_REF] Athreya | Some limit theorems for positive recurrent branching Markov chains I[END_REF][START_REF] Athreya | Some limit theorems for positive recurrent branching Markov chains II[END_REF][START_REF] Benjamini | Markov chains indexed by trees[END_REF]) in the symmetric independent case. That is,

∀(i, k) ∈ T × N, P((Z i0 , Z i1 ) = (k 0 , k 1 ) | Z i = k) = P(Z i0 = k 0 | Z i = k)P(Z i0 = k 1 | Z i = k).
But this identity does not hold here since we are interested in unequal sharing of parasites. Guyon [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF] proves limit theorems for a Markov chain indexed by a binary tree where asymmetry and dependence are allowed. His theorem is the key argument to prove the convergence of asymptotic proportions of cells with a given number of parasites here. Indeed, contamination ensures that the process which counts the number of parasites along the random walk on the binary tree of the cell population is ergodic and non trivial (see Section 5). This is the fundamental assumption to use Guyon's law of large numbers. Let us then introduce more precisely this process which gives the number of parasites in a random cell line.

Let (a i ) i∈N be an i.i.d. sequence independent of (Z i ) i∈T such that

P(a 1 = 0) = P(a 1 = 1) = 1/2. ( 2
)
Denote by f (0) (resp f (1) ) the random p.g.f which gives the law of the size of the offspring of a parasite which goes in the first daughter cell (resp. in the second daughter cell):

f (0) (s) := f (s, 1) a.s., f (1) (t) := f (1, t) a.s., (s, t ∈ [0, 1]).
Let f be the mixed generating function of f (0) and f (1) , i.e.

P(f ∈ dg) = P(f (0) ∈ dg) + P(f (1) ∈ dg) 2 .

Then (Z n ) n∈N = (Z (a 1 ,a 2 ,..an) ) n∈N is a Branching Process in Random Environment with immigration depending on the state is zero or not: the reproduction law is given by its p.g.f f , the immigration law in zero is distributed as Y 0 , and the immigration law in k ≥ 1 as Y 1 . Thus, we first need to prove asymptotic results for this process.

Main results

Galton Watson processes with immigration are well known (see e.g. [START_REF] Asmussen | Branching processes[END_REF][START_REF] Lyons | Probability on Trees and Networks[END_REF]). If the process is subcritical and the expectation of the logarithm of the immigration is finite, then it converges in distribution to a finite random variable. Otherwise it tends to infinity in probability. Key [START_REF] Key | Limiting Distributions and Regeneration Times for Multitype Branching Processes with Immigration in a Random Environment[END_REF] has obtained the analogue result for Branching Processes in Random Environment with Immigration (IBPRE), in the subcritical case, with finite expectation of the logarithm. Actually he states results for multitype IBPRE, which have been complemented by Roitershtein [START_REF] Roithershtein | A note on multitype branching processes with immigration in a random environment[END_REF] who obtained a strong law of large numbers and a central limit theorem for the partial sum.

In Section 4, we give the asymptotic behavior of IBPRE in the different cases, which means that we also consider the critical or supercritical case and the case when the expectation of the logarithm of the immigration is infinite. To get these results, we use some general statements on Markov processes (Section 3.2), classical arguments for Galton Watson process with immigration (see [START_REF] Lyons | Probability on Trees and Networks[END_REF], which was inspired from [START_REF] Asmussen | Branching processes[END_REF]) and the tail of the time when IBPRE returns to 0 in the subcritical case, which is proved in [START_REF] Key | Limiting Distributions and Regeneration Times for Multitype Branching Processes with Immigration in a Random Environment[END_REF].

We can then state results about branching processes in random environment (Z n ) n∈N with immigration depending on the state is zero or not (Section 5) using coupling arguments and Section 3.2. This process gives the number of parasites along a random cell line. Recalling that immigration in state zero is distributed as Y 0 and immigration in state k ≥ 1 is distributed as Y 1 , we prove the following expected result.

Theorem. (i) If E log(f ′ (1)) < 0 and max(E(log + (Y i )) : i = 0, 1) < ∞, then there exists a finite r.v. Z ∞ such that for every k ∈ N, Z n starting from k converges in distribution to Z ∞ as n → ∞. (ii) If E log(f ′ (1)) ≥ 0 or max(E(log + (Y i )) : i = 0, 1) = ∞, then Z n converges in probability to ∞ as n → ∞.
With additional assumptions, we provide in Section 5 an estimate of the rate of convergence of (Z n ) n∈N depending on the initial state.

Then, in Section 6, we prove asymptotic results on the population of cells in generation n as n → ∞.

First, we consider the case when there is no contamination: Y 0 = Y 1 = 0 a.s. We determine when the organism recovers, meaning that the number of infected cells becomes negligible compared to the total number of cells. As stated in Proposition 2, the recovery occurs a.s. iff E(log(f ′ (1))) ≤ 0. Thus, we generalize results of Section 3 in [START_REF] Bansaye | Proliferating parasites in dividing cells : Kimmel's branching model revisited[END_REF] to random environment. Again, for any reproduction rate of parasites, we can find a necessary and sufficient condition on sharing of their offspring so that the organism recovers a.s.

As explained in introduction, a natural example is the random binomial repartition of parasites. If the reproduction of parasites is given by the r.v. Z and the random parameter of the binomial repartition is given P ∈ [0, 1], the a.s. recovery criterion becomes log(E(Z)) ≤ E(log(1/P )).

Second we take into account the contamination by parasites from outside the cell population with assumptions (1). We focus on proportions of cells in generation n with a given number of parasites:

F k (n) := #{i ∈ G n : Z i = k} 2 n (k ∈ N).
Using [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF] and the theorem above, we prove the following law of large numbers.

Theorem. If E(log(f ′ (1))) < 0 and max(E(log + (Y i )) : i = 0, 1) < ∞, then for every k ∈ N, F k (n) converges in probability to a deterministic number f k as n → ∞, such that f 0 > 0 and ∞ k=0 f k = 1. Otherwise, for every k ∈ N, F k (n) converges in probability to 0 as n → ∞.
Finally, in Section 7, we give the asymptotic behavior of the total number of parasites in generation n in the case when the growth of parasites follows a Galton Watson process and the contamination does not depend on the state of cell.

Preliminaries

We recall first some results about Branching Processes in Random Environment (BPRE) and then about Markov chains, which will be both useful to study BPRE with immigration (Z n ) n∈N . Recall that we denote by k the initial number of parasites and by P k the probability associated with.

Branching Processes in Random Environment (BPRE)

We consider here a BPRE (Z n ) n∈N specified by a sequence of i.i.d. generating functions (f n ) n∈N distributed as f [START_REF] Athreya | On branching processes with random environments, I : extinction probability Ann[END_REF][START_REF] Athreya | On branching processes with random environments, II : limit theorems[END_REF][START_REF] Smith | On branching processes in random environments[END_REF]. More precisely, conditionally on the environment (f n ) n∈N , particles at generation n reproduce independently of each other and their offspring has generating function f n . Then Z n is the number of particles at generation n and Z n+1 is the sum of Z n independent random variables with generating function f n . That is, for every n ∈ N,

E s Z n+1 |Z 0 , . . . , Z n ; f 0 , . . . , f n = f n (s) Zn (0 ≤ s ≤ 1).
Thus, denoting by

F n := f 0 • • • • • f n-1 , we have for every k ∈ N, E k (s Z n+1 | f 0 , ..., f n ) = E(s Z n+1 | Z 0 = k, f 0 , ..., f n ) = F n (s) k (0 ≤ s ≤ 1).
When the environments are deterministic (i.e. f is a deterministic generating function), this process is the Galton Watson process with reproduction law N , where f is the generating function of N .

The process (Z n ) n∈N is called subcritical, critical or supercritical if

E log(f ′ (1))
is negative, zero or positive respectively. This process becomes extinct a.s.:

P(∃n ∈ N : Z n = 0) = 1
iff it is subcritical or critical [START_REF] Athreya | On branching processes with random environments, I : extinction probability Ann[END_REF] (see [START_REF] Geiger | Limit Theorems for subcritical branching processes in random environment[END_REF] for finer results).

In the critical case, we make the following integrability assumption:

0 < E(log(f ′ 0 (1)) 2 ) < ∞, E [1 + log(f ′ 0 (1))]f ′′ 0 (1)/2f ′ 0 (1) < ∞,
so that there exist 0 < c 1 < c 2 < ∞ such that for every n ∈ N (see [START_REF] Kozlov | The asymptotic behavior of the probability of non-extinction of critical branching processes in a random environment[END_REF])

c 1 / √ n ≤ P(Z n > 0) ≤ c 2 / √ n. (3)
See [START_REF] Afanasyev | Criticality for branching processes in random environment[END_REF] for more general result in the critical case.

Markov chains

We consider now a Markov chain (Z n ) n∈N taking values in N and introduce the first time T 0 when (Z n ) n∈N visits 0 after time 0:

T 0 := inf{i > 0 : Z i = 0}.

Denote by

u n := P 0 (Z n = 0),

u ∞ := 1/E 0 (T 0 ) (1/∞ = 0).
By now, we assume 0 < P 0 (Z 1 = 0) < 1 and we give the asymptotic behavior of (Z n ) n∈N . The first part of (i) is the classical ergodic property for an aperiodic positive recurrent Markov chain and we provide an estimate of the speed of convergence depending on the initial state. Then (ii) gives the null recurrent case, which is also a classical result.

Lemma 1. (i) If for every k ∈ N, E k (T 0 ) < ∞, then Z n starting from k converges in distribution to a finite random variable Z ∞ ,

which does not depend on k and verifies

P(Z ∞ = 0) > 0.
Moreover there exists A > 0 such that for all n, k ∈ N,

l∈N |P k (Z n = l) -P(Z ∞ = l)| ≤ A sup n/2≤l≤n {|u l -u ∞ |} + E 0 (T 0 1l T 0 >n/4 ) + E k (T 0 1l T 0 >n/4 ) . ( 4
) (ii) If E 0 (T 0 ) = ∞ and for every l ∈ N, P l (T 0 < ∞) > 0, then for every k ∈ N, Z n → ∞ in P k -probability as n → ∞.
Proof of (i). First, note that by the Markov property, for every n ∈ N,

|P k (Z n = 0) -u ∞ | = | n j=1 P k (T 0 = j)P 0 (Z n-j = 0) -u ∞ | ≤ n j=1 P k (T 0 = j)|u n-j -u ∞ | + u ∞ P k (T 0 > n). (5)
On the event {T 0 ≤ n}, define R n as the last passage time of (Z n ) n∈N by 0 before time n:

R n := sup{i ≤ n : Z i = 0}.
For all 0 ≤ i ≤ n and l ∈ N, by the Markov property,

P k (Z n = l) = P k (T 0 > n, Z n = l) + n i=0 P k (T 0 ≤ n, R n = n -i, Z n = l) = P k (T 0 > n, Z n = l) + n i=0 P k (Z n-i = 0)P 0 (Z i = l, T 0 > i).
Define now

α l := u ∞ ∞ i=0 P 0 (Z i = l, T 0 > i).
We then have

|P k (Z n = l) -α l | ≤ P k (T 0 > n, Z n = l) + u ∞ ∞ i=n+1 P(Z i = l, T 0 > i) + n i=0 P(Z i = l, T 0 > i) u ∞ -P k (Z n-i = 0) . Summing over l leads to l∈N |P k (Z n = l) -α l | ≤ P k (T 0 > n) + u ∞ E 0 (T 0 1l T 0 >n+1 ) + n i=0 P(T 0 > i) u ∞ -P k (Z n-i = 0) . (6)
Moreover using (5), we have for all 0 ≤ n 0 ≤ n,

n i=0 P(T 0 > i) u ∞ -P k (Z n-i = 0) ≤ n i=0 P 0 (T 0 > i)   n-i j=1 P k (T 0 = j)|u n-i-j -u ∞ | + u ∞ P k (T 0 > n -i)   ≤ n i=0 P 0 (T 0 > i) n-i j=1 P k (T 0 = j)|u n-i-j -u ∞ | + u ∞ n i=0 P 0 (T 0 > i)P k (T 0 > n -i). (7)
Finally, denoting by M := sup n∈N {|u n -u ∞ |},

n i=0 P 0 (T 0 > i) n-i j=1 P k (T 0 = j)|u n-i-j -u ∞ | ≤ sup n 0 ≤l≤n {|u l -u ∞ |} n i=0 P 0 (T 0 > i) n-i j=1 P k (T 0 = j)1 n-i-j≥n 0 +M n i=0 P 0 (T 0 > i) n-i j=1 P k (T 0 = j)1 n-i-j<n 0 ≤ sup n 0 ≤l≤n {|u l -u ∞ |} n i=0 P 0 (T 0 > i) n-i j=1 P k (T 0 = j) +M n-n 0 i=0 P 0 (T 0 > i)P k (T 0 > n -n 0 -i). (8)
Combining ( 6), ( 7) and ( 8) and using that

n i=0 P 0 (T 0 > i)P k (T 0 > n -i) ≤ E 0 (T 0 1l T 0 >n/2 ) + E k (T 0 1l T 0 >n/2 ), n-n 0 i=0 P 0 (T 0 > i)P k (T 0 ≥ n -n 0 -i) ≤ E 0 (T 0 1l T 0 >(n-n 0 )/2 ) + E k (T 0 1l T 0 >n-n 0 )/2 ),
we get, for all 0 ≤ n 0 ≤ n,

l∈N |P k (Z n = l) -α l | ≤ P k (T 0 > n) + u ∞ E 0 (T 0 1l T 0 >n+1 ) + sup n 0 ≤l≤n {|u l -u ∞ |}E 0 (T 0 ) +[u ∞ + M ][E 0 (T 0 1l T 0 >(n-n 0 )/2 ) + E k (T 0 1l T 0 >(n-n 0 )/2 )]. ( 9
)
As P 0 (Z 1 = 0) > 0, by the renewal theorem [START_REF] Feller | An introduction to probability theory and its applications[END_REF],

u n n→∞ -→ u ∞ . Adding that E k (T 0 ) < ∞ and E 0 (T 0 ) < ∞ ensures that l∈N |P k (Z n = l) -α l | n→∞ -→ 0,
which proves that Z n starting from k converges in distribution to a r.v. Z ∞ which does not depend on k.

The inequality of (i) is obtained by letting n 0 = n/2 in (9).

Proof of (ii)

. If E 0 (T 0 ) = ∞, then by the renewal theorem again [START_REF] Feller | An introduction to probability theory and its applications[END_REF],

u n n→∞ -→ 0. So D n = inf{k -n : k ≥ n, Z k = 0} n→∞ -→ ∞, in probability.
Assume that there exist l ∈ N, ǫ > 0 and an increasing sequence of integers (u n ) n∈N such that P k (Z un = l) ≥ ǫ.

As P l (T 0 < ∞) > 0 by hypothesis, there exists N > 0 such that

P l (T 0 = N ) > 0.
Thus, by the Markov property,

P k (Z un+K = 0) ≥ P k (Z un = l)P l (T 0 = N ) ≥ ǫP l (T 0 = N ).
Then, for all n ∈ N,

P k (D un ≤ N ) ≥ ǫP l (T 0 = N ) > 0,
which is in contradiction with the fact that

D n → ∞ in P k as n → ∞. Then, P k (Z n = l) → 0 as n → ∞.
4 Branching processes in random environment with immigration (IBPRE)

We consider here a BPRE (Z n ) n∈N whose reproduction law is given by the random p.g.f f and we add at each generation n + 1 a random number of immigrants Y n independent and identically distributed as a r.v Y such that

P(Y = 0) > 0.
More precisely, for every n ∈ N,

Z n+1 = Y n + Zn i=1 X i , (10)
where (X i ) i∈N , Y n and Z n are independent and conditionally on f n = g, the (X i ) i∈N are i.i.d. with common probability generating function g.

Note that if the contamination does not dependent on the fact that this cell is already infected or not (i.e. Y 0 and Y 1 are identically distributed), then the number of parasites in a random cell line defined in Introduction is a IBPRE whose reproduction law given by f and immigration by

Y d = Y 0 d = Y 1 .
We give now the asymptotic behavior of this process. These results are classical for the Galton Watson process with immigration [START_REF] Asmussen | Branching processes[END_REF][START_REF] Lyons | Probability on Trees and Networks[END_REF]. We follow the same method in the case of random environment for the subcritical and supercritical cases. We give in (ii) the tail of the time

T 0 = inf{n > 0 : Z n = 0}
when the process returns to 0 in the subcritical case, which is proved in [START_REF] Key | Limiting Distributions and Regeneration Times for Multitype Branching Processes with Immigration in a Random Environment[END_REF] and we use Section 3.2 for the critical case.

Proposition 1. (i) If E log(f ′ (1)) < 0 and E(log + (Y )) < ∞, then Z n converges in distribution to a finite random variable as n → ∞ and lim n→∞ P(Z n = 0) > 0.

Otherwise Z n → ∞ as n → ∞.

(ii) If E log(f ′ (1)) < 0 and there exists q > 0 such that E(Y q ) < ∞, then there exist c, d > 0 such that for every n ∈ N,

P(T 0 > n) ≤ ce -dn . (iii) Assume E(f ′ (1) -1 ) < 1 and E(log + (Y )) < ∞, then there exists a finite r.v. W such that [Π n-1 i=0 f ′ i (1)] -1 Z n n→∞ -→ W, in P.
Note also that by the Borel-Cantelli lemma, if

E(log + (Y 1 )) = ∞, then for every c > 1, lim sup n→∞ c -n Z n = ∞ a.s.
since Z n ≥ Y n a.s. Moreover the proof of Section 5 provides an other approach to prove that (

Z n ) n∈N tends to ∞ if E(log + (Y )) = ∞.
Proof of (i) and (ii) in the subcritical case: E log(f ′ (1)) < 0. The subcritical case with assumption E(log + (Y )) < ∞ is handled in [START_REF] Key | Limiting Distributions and Regeneration Times for Multitype Branching Processes with Immigration in a Random Environment[END_REF]: First part of (i) is Theorem 3.3 and (ii) is a consequence of Theorem 4.2 of [START_REF] Key | Limiting Distributions and Regeneration Times for Multitype Branching Processes with Immigration in a Random Environment[END_REF].

We focus now on the case E(log + (Y )) = ∞ and prove that Z n converges in probability to ∞. The proof is close to the Galton Watson case (see [START_REF] Asmussen | Branching processes[END_REF] or [START_REF] Lyons | Probability on Trees and Networks[END_REF]). First, by Borel-Cantelli lemma, lim sup

k→∞ log + (Y k )/k = ∞ a.s.
Then, for every c ∈ (0, 1), lim sup

k→∞ c k Y k = ∞ a.s. (11)
Note that

Z n = n-1 k=0 Z k,n ,
where Z k,n is the number of descendants in generation n of immigrants in generation n -k. Thus, denoting by Y k,n the number of immigrants in generation n -k and X i (k, n) the number of descendants in generation n of immigrant i in generation n -k, we have

Z n = n-1 k=0 Y k,n i=1 X i (k, n).
This sum increases stochastically as n tends to infinity and converges in distribution to

Z ∞ = ∞ k=0 Y k i=1 X i (k),
where conditionally on (f i : i ∈ N), (X i (k) : i ∈ N, k ∈ N) are independent and the probability generating function of X i (k) is equal to f k-1 • ... • f 0 . Roughly speaking, X i (k) is the contribution of immigrant i which arrives k generations before 'final time' ∞. The integer X i (k) is the population in generation k of a BPRE without immigration starting from 1.

Assume now that Z ∞ < ∞ with a positive probability. As (X i (k) : k ∈ N, 1 ≤ i ≤ Y k ) are integers, then conditionally on Z ∞ < ∞, only a finite number of them are positive. Thus, by Borel-Cantelli lemma, conditionally on (

Z ∞ < ∞, Y k : k ∈ N, f i : i ∈ N), ∞ k=0 Y k P(X 1 (k) > 0) < ∞ a.s.
Moreover, by convexity, for all g p.g.f and s ∈ [0, 1],

1 -g(s) 1 -s = g(1) -g(s) 1 -s ≥ g(1) -g(0) 1 -0 = 1 -g(0), (0 ≤ s ≤ 1).
Then 1 -g(s) ≥ (1 -g(0))(1 -s) and by induction, we have for every k ∈ N,

P(X 1 (k) > 0 | f i : i ∈ N) = 1 -f k-1 • ... • f 0 (0) ≥ Π k-1 i=0 (1 -f i (0)) = exp(S k ), where S k := k-1 i=0 log(1 -f i (0)). Thus, conditionally on (Z ∞ < ∞, Y k : k ∈ N, f i : i ∈ N), ∞ k=0 Y k exp(S k ) < ∞ a.s.
Thus, on the event {Z ∞ < ∞} which has a positive probability, we get

∞ k=0 Y k exp(S k ) < ∞ a.s.
Moreover S n is a random walk with negative drift E(log(1 -f 0 (0))). So letting α < E(log(1 -f 0 (1))), P(S n < αn) decreases exponentially by classical large deviation results.

Then by Borel-Cantelli lemma, S n is less than αn for a finite number of n, and

L := inf n∈N {S n -αn} > -∞ a.s.
Using that for every k ∈ N, S k ≥ αk + L a.s., we get

∞ k=0 exp(αk)Y k < ∞,
with positive probability. This is in contradiction with [START_REF] Feller | An introduction to probability theory and its applications[END_REF]. Then Z ∞ = ∞ a.s. and Z n converges in probability to ∞ as n → ∞.

Proof of (i) in the critical and supercritical case: E log(f ′ (1)) ≥ 0. First, we focus on the critical case. Recall that T 0 = inf{i > 0 : Z i = 0} and consider ( Zn ) n∈N the BPRE associated with (Z n ) n∈N , that is the critical BPRE with reproduction law f and no immigration. Thanks to (3), there exists c 1 > 0 such that for ever n ∈ N,

P 1 ( Zn > 0) ≥ c 1 / √ n.
Adding that

P 1 (T 0 > n) = P 1 (Z n > 0) ≥ P 1 ( Zn > 0), ensures that E 1 (T 0 ) = ∞. Then E 0 (T 0 ) = ∞ since IBPRE (Z n ) n∈N starting from 1 is stochastically larger than (Z n ) n∈N starting from 0. Moreover ∀k ∈ N, P k (T 0 < ∞) > 0, since P k ( T0 < ∞) = 1 and P(Y = 0) > 0. Then Lemma 1 (ii) ensures that Z n → ∞ in P as n → ∞.
For the supercritical case, follow the proof in the critical case (or use the result with a coupling argument) to get that Z n → ∞ in probability as n → ∞ Proof of (iii). We follow again [START_REF] Lyons | Probability on Trees and Networks[END_REF].

If E(log + (Y )) < ∞, by Borel-Cantelli Lemma lim sup k→∞ log + (Y k )/k = 0. Then for every c > 1, ∞ k=0 c -k Y k < ∞ a.s. ( 12
) Define P n := [Π n-1 i=0 f ′ i (1)] -1
, and denote by F n the σ-field generated by (Z i : 0 ≤ i ≤ n), (P i : 0 ≤ i ≤ n) and (Y k : k ∈ N). Then using [START_REF] Geiger | Limit Theorems for subcritical branching processes in random environment[END_REF], we have

E(P n+1 Z n+1 | F n ) = E(P n+1 [ Zn i=1 X i + Y n ] | F n ) = P n E(f ′ n (1) -1 Zn i=1 X i | F n ) + P n E(f ′ (1) -1 )Y n = P n E(f ′ n (1) -1 Z n E(X 1 | f n ) | F n ) + P n E(f ′ (1) -1 )Y n = P n Z n + P n E(f ′ (1) -1 )Y n . So P n Z n is a submartingale. Moreover E(P n Z n | F 0 ) = Z 0 + n-1 i=0 E(f ′ (1) -1 ) i+1 Y i .
By [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF], if E(f ′ (1) -1 ) < 1, P n Z n has bounded expectations and then converges a.s. to a finite r.v.

Ergodicity and convergence for a random cell line

Recall that (Z n ) n∈N defined in Introduction is the number of parasites in a random cell line. The Markov chain (Z n ) n∈N is a BPRE with state dependent immigration. The reproduction law is given by the p.g.f f , immigration in state 0 is distributed as Y 0 and immigration in state k ≥ 1 is distributed as Y 1 . More precisely, for every n ∈ N, conditionally on

Z n = x, Z n+1 = Y (n) x + x i=1 X (n) i ,
where

(i) (X (n) i ) i∈N and Y (n) x are independent. (ii) Conditionally on f n = g, the (X (n) i ) i∈N are i.i.d. with common probability generating function g. (iii) For all x ≥ 1 and n ∈ N, Y (n) x d = Y 1 .
We have the following results, which generalize those of the previous section to the case when immigration depends on whether the state is zero or not.

Theorem 1. (i) If E log(f ′ (1)) < 0 and max(E(log + (Y i )) : i = 0, 1) < ∞, then there exists a finite random variable Z ∞ such that for every k ∈ N, Z n starting from k converges in distribution to Z ∞ as n → ∞.
Moreover, if there exists q > 0 such that max(E(Y q i ) : i = 0, 1) < ∞, then for every ǫ > 0, there exist 0 < r < 1 and C > 0 such that for all n ∈ N and k ∈ N,

∞ l=0 |P k (Z n = l) -P(Z ∞ = l)| ≤ Ck ǫ r n . (ii) If E log(f ′ (1)) ≥ 0 or max(E(log + (Y i )) : i = 0, 1) = ∞, Z n converges in probability to infinity as n → ∞.

Note again that by Borel

-Cantelli lemma, if E(log + (Y 1 )) = ∞, then for every c > 1, lim sup n→∞ c -n Z n = ∞ a.s., since Z n ≥ Y n a.s.
The proof of (ii) in the critical or supercritical case (E log(f ′ (1)) ≥ 0) is directly derived from Proposition 1 and we focus now on the subcritical case:

E log(f ′ (1)) < 0.
Recall that T 0 is the first time after 0 when (Z n ) n∈N visits 0. Using IBPRE (see Section 4), we prove the following result in the subcritical case.

Lemma 2. If max(E(log

+ (Y i )) : i = 0, 1) < ∞, then for every kN, P k (T 0 < ∞) = 1 and sup n∈N {P k (Z n ≥ l)} l→∞ -→ 0.
Moreover if there exists q > 0 such that max(E(Y q i ) : i = 0, 1) < ∞, then for every ǫ > 0, there exist r > 0 and C > 0 such that for all n ∈ N, k ≥ 1:

P 0 (T 0 ≥ n) ≤ Cr n , P k (T 0 ≥ n) ≤ Ck ǫ r n .
Proof. We couple (Z n ) n∈N with an IBPRE ( Z n ) n∈N with reproduction law given by the random p.g.f f (such as (Z n ) n∈N ) and immigration Y defined by

Y := max(Y 0 , Y 1 , Y ),
where Y 0 , Y 1 , and Y are independent and Y is defined by

P( Y = 0) = 1/2; ∀n ∈ N * , P( Y = n) = αn -1-ǫ , α := [2 ∞ i=1 i -1-ǫ ] -1 .
Thus immigration Y for Z n is stochastically larger than immigration for Z n (whereas reproduction law is the same), so that coupling gives

∀n ∈ N, Z n ≤ Z n a.s.
Moreover, Z n is still subcritical. Recalling that min(P(Y i = 0) : i = 0, 1) > 0, P( Y = 0) = 1/2, and that the expectation of the logarithm of every r.v. is finite, we have

E(log + (Y )) < ∞, P(Y = 0) > 0.
Then Proposition 1 (i) ensures that Z n converges in distribution to a finite random variable, so that sup

n∈N {P k (Z n ≥ l)} ≤ sup n∈N {P k ( Z n ≥ l)} l→∞ -→ 0.
Proposition 1 (i) ensures also that for every k ∈ N, lim n→∞ P k (Z n = 0) > 0. Thus, for every k ∈ N, P k ( T 0 < ∞) = 1 and then P k (T 0 < ∞) = 1. This completes the first part of the lemma.

We assume now that there exists q > 0 such that max(E(Y q i ) : i = 0, 1) < ∞. Moreover E( Y ǫ/2 ) < ∞, so letting q ′ = min(ǫ/2, q), we have

E(Y q ′ ) < ∞.
We can then apply Proposition 1 (ii) to IBPRE ( Z n ) n∈N , so that there exist c, d > 0 such that for every n ∈ N, P 0 ( T 0 > n) ≤ ce -dn .

Recalling that for all k, n ∈ N,

P k (T 0 ≥ n) ≤ P k ( T 0 ≥ n),
we get P 0 (T 0 ≥ n) ≤ ce -dn . Moreover for every k ∈ N,

P 0 ( T 0 > n) ≥ P(Y ≥ k)P k ( T 0 ≥ n).
By definition of Y , there exists β > 0 such that for every n ∈ N,

P(Y ≥ n) ≥ βn -ǫ .
Using these inequalities gives

P k (T 0 ≥ n) ≤ P k ( T 0 ≥ n) ≤ β -1 k ǫ P 0 ( T 0 > n) ≤ β -1 ck ǫ e -dn .
This completes the proof.

Proof of Theorem 1 (i) and (ii) in the subcritical case: E log(f ′ (1)) < 0. We split the proof into 4 cases:

CASE 1: max(E(log + (Y i )) : i = 0, 1) < ∞.
CASE 2: There exists q > 0 such that max(E(Y q i ) : i = 0, 1) < ∞.

CASE 3: E(log + (Y 1 )) = ∞. CASE 4: E(log + (Y 0 )) = ∞.
First, note that P(Y 0 = 0) > 0 ensures that P 0 (Z 1 = 0) > 0 and we can use results of Section 3.2. CASE 1. In this case, by Lemma 2, (Z n ) n∈N is bounded in distribution:

sup n∈N {P 0 (Z n ≥ l)} l→∞ -→ 0.
If E 0 (T 0 ) = ∞, then Z n → ∞ in P 0 by Lemma 1 (ii), which is in contradiction with the previous limit. Then E 0 (T 0 ) < ∞. We prove now that ∀k ≥ 1, E k (T 0 ) < ∞ by a coupling argument. Let k ≥ 1 and change only immigration to get a Markov process ( Z n ) n∈N which is larger than

(Z n ) n∈N : ∀n ∈ N, Z n ≥ Z n a.s.
Its immigrations Y 0 and Y 1 satisfy

Y 1 d = Y 1 , ∀n ∈ N, P( Y 0 ≥ n) ≥ P(Y 0 ≥ n), P( Y 0 ≥ k) > 0, max(E(log( Y i ) : i = 0, 1) < ∞.
Then, we have again E 0 ( T 0 ) < ∞, which entails that E k ( T 0 ) < ∞ since P( Y 0 ≥ k) > 0.

As for every n ∈ N, Z n ≥ Z n a.s., we have

E k (T 0 ) ≤ E k ( T 0 ) < ∞.
Then Lemma 1 (i) ensures that for every k ∈ N, (Z n ) n∈N converges in distribution to a finite random variable Z ∞ , which does not depend on k and verifies P(Z ∞ = 0) > 0.

CASE 2: By Lemma 1 (i), we have

l∈N |P k (Z n = l) -P(Z ∞ = l)| ≤ A sup n/2≤l≤n
{|u l -u ∞ |} + E 0 (T 0 1l T 0 >n/4 ) + E k (T 0 1l T 0 >n/4 ) .

(13) Moreover by Lemma 2, for every ǫ > 0, there exists C > 0 such that

P k (T 0 ≥ n) ≤ Ck ǫ r n , P 0 (T 0 ≥ n) ≤ Cr n . ( 14
)
So for every r ′ ∈ (r, 1), E 0 (exp(-log(r)T 0 )) < ∞. Then, by Kendall renewal theorem [START_REF] Kendall | Unitary dilations of Markov transition operators, and the corresponding integral representations for transition-probability matrices[END_REF], there exists ρ ∈ (0, 1) and c > 0 such that for every n ∈ N,

|u n -u ∞ | ≤ cρ n . (15)
Finally, [START_REF] Kendall | Unitary dilations of Markov transition operators, and the corresponding integral representations for transition-probability matrices[END_REF] ensures that there exists D > 0 such that for every n ∈ N,

E 0 (T 0 1l T 0 >n/4 ) ≤ Dnr n/4 , E k (T 0 1l T 0 >n/4 ) ≤ Dnk ǫ r n/4 .
Combining these two inequalities with ( 13) and ( 15), we get

l∈N |P k (Z n = l) -P(Z ∞ = l)| ≤ A cρ n + Dnr n/4 + Dnk ǫ r n/4 ,
which ends the proof in CASE 2. CASE 3. Change immigration of (Z n ) n∈N to get an IBPRE ( Z n ) n∈N whose immigration is distributed as Y 1 and whose reproduction law is still given by f . Then Proposition 1 (i) and E(log + (Y 1 )) = ∞ ensures that ( Z n ) n∈N starting from 0 tends in distribution to ∞. Then Lemma 1 (i) entails that E 0 ( T 0 ) = ∞, so that for every k ≥ 1,

E k ( T 0 ) ≥ E 0 ( T 0 ) = ∞, since the IBPRE ( Z n ) n∈N starting from k ≥ 1 is stochastically larger than ( Z n ) n∈N starting from 0. Moreover, under P k , (Z n ) n∈N is equal to ( Z n ) n∈N until time T 0 = T 0 . So E k (T 0 ) = ∞. Let k ≥ 1 such that P 0 (Z 1 = k) > 0, then E 0 (T 0 ) ≥ P 0 (Z 1 = k)E k (T 0 -1). This entails that E 0 (T 0 ) = ∞.
By Lemma 1 (ii), (Z n ) n∈N starting from any k ∈ N tends to ∞ in probability. CASE 4. Denote by

X i := P(Z i > 0 | Z i-1 = 1, f i-1
), (i ≥ 1), the survival probability in environment f i-1 and introduce the following random walk

S n = n i=1 log(X i ). Then P 1 (Z n > 0 | (f 0 , f 1 , ..., f n-1 )) ≥ Π n 1 X i = exp(S n ) a.s., so that P k (Z n > 0 | (f 0 , f 1 , ..., f n-1 )) = 1 -P k (Z n = 0 | (f 0 , f 1 , ..., f n-1 )) = 1 -[1 -P 1 (Z n > 0 | (f 0 , f 1 , ..., f n-1 ))] k ≥ 1 -[1 -exp(S n )] k a.s. Thus P k (Z n > 0) ≥ E(1 -[1 -exp(S n )] k ).
Using the Markov property we have

E 0 (T 0 + 1) ≥ ∞ k=1 P(Y 0 = k)E k (T 0 ) = ∞ k=1 P(Y 0 = k) ∞ n=1 P k (T 0 ≥ n) ≥ ∞ k=1 P(Y 0 = k) ∞ n=1 P k (Z n > 0) ≥ ∞ k=1 P(Y 0 = k) ∞ n=1 E(1 -[1 -exp(S n )] k ).
Moreover for all x ∈ [0, 1[ and k ≥ 0, exp(k log(1 -x)) ≤ exp(-kx), and by the law of large numbers, S n /n tends a.s. to E(X 1 ) < 0 so that there exists n 0 ≥ 1 such that for every n ≥ n 0 , P(S n /n ≥ 3E(X 1 )/2) ≥ 1/2.

We get then

E 0 (T 0 + 1) ≥ ∞ k=1 P(Y 0 = k) ∞ n=1 E(1 -exp(-k exp(S n ))) ≥ [1 -e -1 ] ∞ n=1 ∞ k=1 P(k exp(S n ) ≥ 1)P(Y 0 = k) ≥ [1 -e -1 ] ∞ n=n 0 P(S n /n ≥ 3E(X 1 )/2) ∞ k≥exp(-3nE(X 1 )/2) P(Y 0 = k) ≥ 2 -1 [1 -e -1 ] ∞ n=n 0 P(Y 0 ≥ exp(-3nE(X 1 )/2)) ≥ 2 -1 [1 -e -1 ] ∞ n=n 0 P(β log(Y 0 ) ≥ n), where β := [-3E(X 1 )/2] -1 > 0. Then E(log(Y 0 )) = ∞ ensures that E 0 (T 0 + 1) = ∞, so E 0 (T 0 ) = ∞.
Conclude that (Z n ) n∈N tends to ∞ in P k using Lemma 1 (ii).

6 Asympotics for proportions of cells with a given number of parasites

Asymptotics without contamination

Here there is no contamination, i.e. Y 0 = Y 1 = 0 a.s. and we determine when the organism recovers, meaning that the number of contaminated cells becomes negligible compared to the total number of cells. We get the same result as Theorem 1 in [START_REF] Bansaye | Proliferating parasites in dividing cells : Kimmel's branching model revisited[END_REF] for the more general model considered here. Denote by N n the number of contaminated cells.

Proposition 2. N n /2 n decreases as n grows. If E(log(f ′ (1))) ≤ 0, then N n /2 n → 0 a.s. as n → ∞.
Otherwise, N n /2 n → 0 as n → ∞ iff all parasites die out, which happens with a probability less than 1.

Example. Consider the case of the random binomial repartition of parasites mentioned in Introduction. Let Z ∈ N be a r.v and (P i ) i∈T be an i.i.d. sequence distributed as a r.v. P ∈ [0, 1], such that P d = 1 -P . In every generation, each parasite gives birth independently to a random number of parasites distributed as Z. When the cell i divides, conditionally on P i = p, each parasite of the cell i goes independently in the first daughter cell with probability p (or it goes in the second daughter cell, which happens with probability 1 -p). Then,

P(f ′ (1) ∈ dx) = P(E(Z)P ∈ dx).
Thus, the organism recovers a.s. (i.e. N n /2 n tends a.s. to 0) iff log(E(Z)) ≤ E(log(1/P )). This is the same criteria in the case when the offspring of each parasite goes a.s. is the same daughter cell (there, p is the probability that this offspring goes in the first daughter cell.)

Proof. Note that N n /2 n decreases to L as n → ∞, since one infected cell has at most two 6.2 Asymptotics with contamination in the case E(log(f ′ (1))) < 0 and max(E(log + (Y i )) : i = 0, 1) < ∞.

Define F k (n) the proportion of cells with k parasites in generation n:

F k (n) := #{i ∈ G n : Z i = k} 2 n (k ∈ N).
We introduce the Banach space l 1 (N) and the subset of frequencies S 1 (N) which we endow with the norm . 1 defined by:

l 1 (N) := {(x i ) i∈N : ∞ i=0 |x i | < ∞}, (x i ) i∈N 1 = ∞ i=0 |x i |, S 1 (N) := {(f i ) i∈N : ∀ i ∈ N, f i ∈ R + , ∞ i=0 f i = 1}.
The main argument here is the law of large number proved by Guyon [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF] for asymmetric Markov chains indexed by a tree.

Theorem 2. If E(log(f ′ (1))) < 0 and max(E(log + (Y i )) : i = 0, 1) < ∞, then (F k (n)) k∈N converges in probability in S 1 (N) to a deterministic sequence (f k ) k∈N as n → ∞, such that f 0 > 0 and ∞ k=0 f k = 1. Moreover, for every k ∈ N, f k = P(Z ∞ = k).

Proof. Recall that (Z i ) i∈T is a Markov chain indexed by a tree and we are in the framework of bifurcating Markov chain studied in [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF]. Thanks to the ergodicity of the number of parasites in a random cell line proved in the previous section (Theorem 1 (i)), we can directly apply Theorem 8 in [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF] to get the convergence of proportions of cells with a given number of parasites.

But it seems that we can't apply Theorem 14 or Corollary 15 in [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF] to get a.s. convergence of proportions, because of the term k ǫ in estimation of Theorem 1. For examples, we refer to the previous proposition.

Using again [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF], we can prove also a law of large numbers and a central limit theorem for the proportions of cells with given number of parasites before generation n. Define, for every n ∈ N, Proof. Use again Theorem 1 (i) and Theorem 8 in [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF] to prove the law of large numbers.

P k (n) := #{i ∈ ∪ 0≤i≤n G i : Z i = k} 2 n+1 (k ∈ N).
For the central limit theorem, use Theorem 19 in [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF] by letting F be the set of continuous functions taking values in [0, 1].

6.3 Asymptotics with contamination in the case E(log(f ′ (1))) ≥ 0 or max(E(log + (Y i )) : i = 0, 1) = ∞.

In this case, cells become infinitely infected as the generation tends to infinity. By Theorem 1, P(Z n ≥ K) tends to 1, then 1 -#{i ∈ G n : Z i ≥ K}/2 n converges to 0 in L 1 , which gives the result.

Asymptotics for the number of parasites

We assume here that parasites multiply following a Galton Watson process with deterministic mean m, independently of the cell they belong to. That is, s → f (s, s) is deterministic and every parasite multiply independently with the reproduction law whose p.g.f. is equal to g : s → f (s, s). Moreover we assume that contamination of a cell does not depend on the number of enclosed parasites. That is

Y d = Y 0 d = Y 1 .
Set P n the number of parasites in generation n. Without contamination, in the supercritical case m > 1, it is well know that either P n becomes extinct or P n /m n converges to a positive finite random variable. In the presence of contamination, we have the following result. Proposition 3. If E(Y ) < ∞ and P(Y 0 > 0) > 0, then log(P n )/n converges in P to log(max(2, m)).

If m = 2, then E(P n ) = E(Y )nm n . This gives the upper bound by Markov inequality and completes the proof.

Figure 1 .

 1 Figure 1. Cell division with multiplication of parasites, random sharing and contamination. Each parasite gives birth to a random number of light parasites and dark parasites. Light parasites go into the first daughter cell, dark parasites go into the second daughter cell and square parasites contaminate the cells from outside the cell population. But light/ dark/ square parasites then behave in the same way.

Theorem 3 .

 3 If E(log(f ′ (1))) < 0 and max(E(log + (Y i )) : i = 0, 1) < ∞, then (P k (n)) k∈N converges in probability in S 1 (N) to the deterministic sequence (f k ) k∈N as n → ∞. Moreover for every k ∈ N, √ n(P k (n) -f k )converges in distribution to a centered normal law as n → ∞, with a non explicit variance.

Theorem 4 .

 4 If E(log(f ′ (1))) ≥ 0 or max(E(log + (Y i )) : i = 0, 1) = ∞, for every k ∈ N, then F k (n) tends to zero as n → ∞. That is, for very K ∈ N, lim n→∞ #{i ∈ G n : Z i ≥ K}/2 n P = 1.Proof. By Fubini's theorem, we haveE #{i ∈ G n : Z i ≥ K}/2 n = i∈Gn P(Z i ≥ K)/2 n = i∈Gn P((a 0 , ...a n-1 ) = i)P(Z i ≥ K) = P(Z n ≥ K).

daughter cells which are infected. Moreover, for every n ∈ N,

= i∈Gn P((a 0 , ..., a n-1 ) = i)P(Z i > 0)

= P(Z n > 0).

If E(log(f ′ (1))) ≤ 0 (subcritical or critical case), then P(Z n > 0) tends to 0 as n → ∞ (see Section 3.1). Thus, E(L) = 0 and N n /2 n tends to 0 a.s. as n → ∞.

If E(log(f ′ (1))) > 0 (supercritical case), then P(Z n > 0) tends to a positive value, which is equal to P(L > 0)0. We complete the proof with the following lemma.

Let us prove the following zero one law, where P n is the total number of parasites in generation n.

Proof. First, we prove that conditionally on non-extinction of parasites, for every K ∈ N, there exists a.s. a generation n such that N n ≥ K. Letting K ∈ N, we fix p as the first integer such that 2 p ≥ K. Then q := P 1 (N p ≥ K) > 0 since P(N 1 = 2) > 0. Either the number of infected cells in generation p is more than K, which happens with probability q, or we can choose in generation p an infected cell i(1), since parasites have not died out. Then, with probability larger than q , the number of infected cells in generation p of the subtree rooted in this cell i(1) contains more than K parasites. Note that this probability is exactly equal to q iff the infected cell i(1) contains one single parasite. Recursively, we find a.s. a generation n such than N n ≥ K.

Then, recalling that we still work conditionally on non-extinction of parasites, the stopping time T := inf{n ∈ N : N n ≥ K} < ∞ a.s. We now also condition by T = n and N T = k. We can then choose one parasite in every infected cell in generation n, which we label by 1 ≤ i ≤ k and we denote by

p the number of cells in generation n + p infected by parasites whose ancestor in generation n is the parasite i. By branching property, the integers (N p (i) :

As sup(L (i) : i ∈ N) = ∞ a.s., letting K → ∞ ensures that a.s. N p /2 p does not tend to 0.

Proof. First, we prove the lower bound. This is a consequence of the fact that P n is larger than (i) the total number of parasites P 1 n which contaminate cells of generation n, (ii) the number of parasites P n (p) in generation n with the same given parasite ancestor in generation p.

Indeed, first P 1 n is the sum of 2 n i.i.d. random variables with mean E(Y ), so law of large numbers ensures that

Moreover for every p < n,

with P(W > 0) > 0. Let now P 2 n be the sum of the number of descendants in generation n of each parasite of generation p. We get then the sum of P p i.i.d. quantities distributed as P n (p). Then ( 16) and ( 17) ensure that we can choose p such that

Using that N n is larger than P 1 n and P 2 n ensures that for every ǫ > 0, lim sup n→∞ P(log(P n )/n ≤ log(max(2, m))) < ǫ.

Letting ǫ → 0 gives the lower bound.

Second, we prove the upper bound. Note that the total number of parasites in generation n can be written as

where Y i,j is the number of parasites which contaminate the jth cell of generation i, and labeling by 1 ≤ k ≤ Y i,j these parasites, Z i,j k is the number of descendants in generation n of the kth parasites. Moreover (Y i,j : i ∈ N, j ∈ N) are identically distributed and independent of (Z i,j p (k), i ∈ N, j ∈ N, k ∈ N), (Z i,j k , i ∈ N, j ∈ N, k ∈ N) are independent and Z i,j p (k) is the population of a Galton Watson process in generation n -i with offspring probability generation function equal to g. Thus


Florence Le Ber

Claire Lavigne

Katarzyna Adamczyk

Frédérique Angevin

Nathalie Colbach

Jean-François Mari

Hervé Monod

Neutral modelling of agricultural

Keywords: Neutral landscape models, agricultural landscapes, gene flow, spatial point process, Voronoi diagrams, rectangular tessellation, field pattern, GenExP-LandSiTes

published or not. The documents may come

Abstract. Neutral landscape models are not frequently used in the agronomical domain, whereas they would be very useful for studying given agro-ecological or physical processes. Contrary to ecological neutral landscape models, agricultural models have to represent and manage geometrical patches and thus should rely on tessellation methods. We present a three steps approach that aimed at simulating such landscapes. Firstly, we characterized the geometry of three real field patterns; secondly, we generated simulated field patterns with two tessellation methods attempting to control the value of some of the observed characteristics and, thirdly, we evaluated the simulated field patterns. For this evaluation, we considered that good simulated field patterns should capture characteristics of real landscapes that are important for the targeted agro-ecological process. Real landscapes and landscapes simulated using either a Voronoi or a rectangular tessellation were thus compared when used as input data within a gene flow model. The results showed that neither tessellation method captured field shapes correctly, thus leading to over or (small) under estimation of gene flow. The

Introduction

Landscape simulation has been a major line of investigation in the fields of ecology and agronomy over the last few years. In landscape ecology, landscape simulation is used for studying the influence of landscape patterns on ecological phenomena [START_REF]Quantitative Methods in Landscape Ecology[END_REF]. Simulation methods fall into three groups: geostatistical models, neutral models, and process-explicit models [START_REF] Saura | Landscape patterns simulation with a modified random clusters method[END_REF]. Geostatistical models focus on interpolation methods of spatial data [START_REF] Gotway | The components of geostatistical simulation[END_REF]. Processexplicit models include the description of spatial processes and related ecological phenomena which underlie the landscape dynamics (e.g., [START_REF] Wu | A patch-based spatial modeling approach: conceptual framework and simulation scheme[END_REF][START_REF] Liu | FORMOSAIC: an individual-based spatially explicit model for simulating forest dynamics in landscape mosaics[END_REF]). In contrast, neutral landscape models (NLM, sensu [START_REF] Gardner | Neutral models for the analysis of broad-scale landscape pattern[END_REF]) provide random landscape structures as a baseline for comparison with real landscape patterns, or for an evaluation of landscape structure effects on ecological processes [START_REF] Gardner | Neutral models for the analysis of broad-scale landscape pattern[END_REF]16,[START_REF] With | The use and misuse of neutral landscape models in ecology[END_REF]. Such models are termed neutral as they model no explicit process giving rise to the landscape pattern, they do model totally random or somewhat constrained covering of the area. Neutral models are mainly based on raster approaches, land-use being allocated randomly to pixels and then clustered by various methods (e.g., [START_REF] Saura | Landscape patterns simulation with a modified random clusters method[END_REF]). Fractal models have also been used [START_REF] Hargrove | A fractal landscape realizer for generating synthetic maps[END_REF]. More recently, polygonal approaches, using geometrical tessellations, have been proposed to simulate patchy landscapes [START_REF] Gaucherel | Neutral models for patchy landscapes[END_REF][START_REF] Gaucherel | Neutral models for polygonal landscapes with linear networks[END_REF]. Landscape models have been widely used in forest management (e.g., [START_REF] Liu | FORMOSAIC: an individual-based spatially explicit model for simulating forest dynamics in landscape mosaics[END_REF][START_REF] Kurz | Telsa: the tool for exploratory landscape scenario analyses[END_REF]). In land-use planning and agronomy, such models are less frequent. Spatial models have been used for the simulation of crop allocations on a farm or over a regional territory (e.g., [START_REF] Ber | Modelling the spatial organisation of land use in a farming territory. Example of a village in the Plateau Lorrain[END_REF][START_REF] Carsjens | Strategic land-use allocation: dealing with spatial relationships and fragmentation of agriculture[END_REF]), in order to understand or to plan the organisation of agricultural landscapes according to environmental problems (e.g., water pollution). Closer to ecological approaches, other studies have focused on the links between farming practices and ecological processes with process-explicit models [START_REF] Gaucherel | A domain-specific language for patchy landscape modelling: The brittany agricultural mosaic as a case study[END_REF][START_REF] Aviron | Conservation of butterfly populations in dynamic landscapes: The role of farming practices and landscape mosaic[END_REF]. In both cases, the models rely on real spatial data, such as field patterns, landscape characteristics (i.e., soil, slope) or the current agricultural mosaic of a region. As far as we know, there are almost no neutral models, such as those developed in landscape ecology, though these models would be very useful when dealing with models of agro-ecological processes. Several reasons for using neutral models can be cited. First, real data are not always available or are too specific and thus reduce the scope of application of the model results. Second, in the case of anthropogenic landscapes, it is necessary to prospect new configurations in order to forecast their effects or to find the best configuration with respect to the given agro-ecological process. Finally, neutral landscapes can be used to test the sensitivity of process models to the spatial variability of agricultural landscapes.

We propose to simulate agricultural landscapes, including both their configuration (the field pattern) and their composition (occurrence of categories of land-use) [START_REF] Li | A simulation experiment to quantify spatial heterogeneity in categorical maps[END_REF]. However, as mentioned by [START_REF] Gaucherel | A domain-specific language for patchy landscape modelling: The brittany agricultural mosaic as a case study[END_REF], agricultural landscapes cannot be easily simulated with traditional ecological neutral landscape models for two reasons: (1) they are mainly geometrical, contrarily to less anthropogenic landscapes, and (2) in traditional landscape models, the basic unit is the pixel, and the land-use mosaic (i.e., both the configuration and the composition) only emerges as a result of the simulation (e.g., [START_REF] Wu | A patch-based spatial modeling approach: conceptual framework and simulation scheme[END_REF]), whereas in agricultural landscapes, the field pattern is generally stable and "precedes" the land-use. For handling geometric patches as basic units, a straightforward way is to use tessellation methods, such as proposed by [START_REF] Barrett | Voronoï tessellation methods to delineate harvests units for spatial forest planning[END_REF] or [START_REF] Gaucherel | Neutral models for patchy landscapes[END_REF]. However, various tessellation methods are available, that result in somewhat different geometries. Good tessellation methods should have the ability to correctly simulate characteristics of real landscapes that are of importance for the targeted agro-ecological process. In the following, we simulate agricultural landscapes using two different tesselation methods. These landscapes are used as input for a gene flow model that handles pollen dispersal as an agro-ecological process. More precisely, we investigate how the values of simulated cross-pollination rates from GM to non-GM fields depend on the tessellation method used to simulate the landscapes, with the rationale that if cross-pollination rates are similar between real and simulated landscapes, this would mean that simulated landscapes capture essential features of real landscapes. To measure and compare these features is necessary to interpret safely the results of gene flow models on simulated landscapes. Whereas various tessellation methods have been used as a basis for landscape generation (e.g., [START_REF] Barrett | Voronoï tessellation methods to delineate harvests units for spatial forest planning[END_REF][START_REF] Gaucherel | Neutral models for polygonal landscapes with linear networks[END_REF]), as far as we know, very few studies focused on their evaluation via their impact on the outcome of ecological processes (e.g., [START_REF] Haydon | Metapopulation theory, landscape models, and species diversity[END_REF]).

To this end, we followed a two-step procedure: we first generated simulated landscapes with basic properties of actual landscapes; then, we as-sessed the impact of the tessellation methods on gene flow. In Section 2, we describe tessellation methods and an innovative point-process model which can be used for generating tessellation seeds. Section 3 focuses on the combination of tessellation methods and the gene flow model. Results are given in Section 4 and they are discussed in Section 5. The conclusion proposes some promising leads for improving the realism of tessellation-based simulated landscapes.

Neutral models of agricultural landscapes

In this section, we focus on the simulation of agricultural landscapes in terms of configuration (field shape, size, and location) and composition (the land-use mosaic). The configuration has been designed with the help of tessellation methods. The composition can be obtained in a purely random manner, or on the basis of expert knowledge or real data. Furthermore, our objective was to explore the variability of landscape patterns on the basis of some given, controlable parameters describing the basic properties of actual agricultural landscapes. We therefore sought a correspondence between geometrical tessellations and spatial patterns of existing landscapes.

Simulating the geometry of agricultural landscapes

Tessellation entails filling a defined area with geometrical figures (hereafter refered to as cells) without any overlapping and leaving no holes. Tessellation methods are based on spatial statistics and computational geometry [START_REF] Okabe | Spatial tessellations: concepts and applications of Voronoi diagrams[END_REF]. They are used in several fields of study, as in image recognition or geography for example, but more for spatial analysis than for simulation (e.g., [START_REF] Michael | Geospatial Analysis -A Comprehensive Guide to Principles, Techniques and Software Tools[END_REF]). Purely random or regular tessellations were not taken into consideration, and we have focused on approaches that allow the handling of tessellation in terms of given characteristics of real landscapes. These approaches are based on a tessellation algorithm and a given set of points or segments, called the tessellation seeds.

The ability of a tessellation method to preserve the features observed in a real landscape can be evaluated, on the one hand, on the basis of the geometry of the tessellation cells (size, shape, orientation), and, on the other hand, on the basis of their spatial distribution (clustered or regular patterns). In both tessellation methods investigated here, the geometry of the cells is completely determined by the tessellation algorithm. In contrast, the spatial distribution of tessellation cells depends on the distribution of the tessellation seeds. Consequently, a tessellation which is based on the "seeds" of a real landscape will preserve the spatial configuration of the fields in that landscape.

Simulation of tessellation seeds

If tessellation seeds are to be generated in a stochastic manner, they can be considered as a realisation of a spatial (2D) point-process, within a bounded window. Furthermore, if the seeds correspond to specific points of a landscape pattern, then the point-process should fit the spatial distribution of these specific points. Eventually, if examples of such specific points are available, it is possible to estimate the point-process parameters in order to simulate similar configurations of seeds. We describe here an example of a point-process fitted to centroid patterns (see [START_REF] Adamczyk | GenExP, un logiciel simulateur de paysages agricoles pour l'étude de la diffusion de transgènes[END_REF] for more details). We chose a point-process that could yield an aggregated or regular centroid pattern and account for the fact that the distance between a pair of real centroids is greater than a certain threshold, proportional to the minimal field area. The Poisson process, commonly used as the reference point-process, is not appropriate for modelling centroid distribution. The pairwise-interaction process (see for example [START_REF] Van Lieshout | Markov point processes and their applications[END_REF]) is more suitable for our purposes. The model controls both the mean number of points per unit area (intensity) and the distribution of pairwise distances between the points. It is defined by the following density function:

f (x; λ, θ) ∝ λ n i<j γ θ ( x i -x j ) (1)
where x 1 , . . . , x n are the coordinates of n centroids, λ is the intensity of the process and γ θ (•) is the interaction function depending on the Euclidean distance between the pairs of points and the interaction parameter θ. For γ θ ≡ 1 the model is reduced to a homogeneous Poisson process with intensity λ: no interaction between the points is observed. The values of γ θ > 1 introduce an aggregated point pattern and the values of γ θ < 1 introduce a regular point pattern at a given distance.

Two tessellation methods

The Voronoi tessellation Given a set of points in the Euclidean plan, the set of seeds, a Voronoi tessellation -or Voronoi diagram-is a covering of the Euclidean plan with non-overlapping convex polygons, each surrounding a seed. The points in a polygon are closer to its seed than to any other seed located in any other polygon. The edges therefore consist in the points located at an equal distance from two seeds (Figure 1(a)). The Delaunay triangulation [START_REF] Okabe | Spatial tessellations: concepts and applications of Voronoi diagrams[END_REF] can be used to determine the polygons. The major advantage of this tessellation is the correspondence between the pair polygon/seed and the pair field/centroid. For example, equally distanced seeds yield regular pattern of cells, aggregated seeds yield clusters of small cells and aligned seeds result in range-structured cells. Similarly, the configuration of field centroids characterises the spatial distribution of the fields: for instance, an aggregated pattern of centroids reveals clusters of small fields. Therefore, a simulation based on statistic parameters that are estimated on centroids of real landscapes is possible. The drawbacks of this tessellation are related to the shape of the polygons which are convex and often have a higher number of vertices than actual fields.

A rectangular tessellation In a rectangular tessellation, the area is filled in with non-overlapping rectangles. We eliminated the simple case when the rectangles are defined by two orthogonal sets of parallel lines and considered only those rectangles sharing T-vertices (Figure 1(b)). In such tessellations, a vertex of a rectangle cannot belong to more than three rectangles (no X-vertex). A method for building such tessellations has been described in [START_REF] Mackisack | Homogeneous rectangular tessellations[END_REF]. The basic principles of the algorithm is to generate the edges of the rectangles by crossing two-directional rays starting from a set of points.

In this case, there is no obvious correspondence between the rectangles and the fields, since the basic geometric unit is a segment rather than a polygon and a seed is located anywhere along an edge. We found no specific points representative of a real landscape to be used as seeds for the rectangular tessellation. Nevertheless, the geometry of the rectangular cells yields a priori more realistic field shapes (for at least some regions, see figure 3(b)) than the Voronoi cells, and certain overall characteristics of a landscape pattern can be conserved: for example, the number of the polygons is equal to the number of seeds plus one.

Simulating the land-use mosaic of agricultural landscapes

The simulation of the land-use mosaic can be handled in several ways. Basically, in most landscape models, crops are randomly allocated according to a probability distribution: e.g., in binary patchy landscapes, each plot is grown with a crop according to a probability p (vs. 1-p) and its own surface area. When several crops are under consideration, each crop is associated with a probability value p c , so that c p c = 1.

Methods have been proposed to improve the purely random approach. In [START_REF] Gaucherel | Neutral models for patchy landscapes[END_REF], a Gibbs process is introduced to model the interactions between pairs of neighbouring patches. According to the parameter values of this model, the landscape mosaic can be more or less heterogeneous. Other processes could be used, based on the statistical study of spatial and temporal patterns of real field mosaics [START_REF] Castellazzi | New measures and tests of temporal and spatial pattern of crops in agricultural landscapes[END_REF][START_REF] Castellazzi | A systematic representation of crop rotations[END_REF]. Furthermore, stochastic neighbourhood rules can be learnt from a real dataset, as proposed in [START_REF] Mari | Temporal and Spatial Data Mining with Second-Order Hidden Markov Models[END_REF], where Hidden Markov Models are applied on land-use temporal and spatial data.

The GenExP-LandSiTes software

GenExP-LandSiTes1 is a software simulating two-dimensional agricultural landscapes. It is written in Java, and it is freely accessible through a Gnu Public Licence [START_REF] Le Ber | GenExP, un logiciel pour simuler des paysages agricoles, en vue de l'étude de la diffusion de transgènes[END_REF]. Its main features fall into three categories: the generation of the field patterns based on the simulation of spatial pointprocesses and on tessellation methods; the post-processing that improves the realistic aspect of the simulated landscape; and the computation of landscape statistics.

The two tessellation methods described above have been implemented: the Voronoi tessellation, using the "3D Hull" algorithm for Delaunay tri-angulation [START_REF] O'rourke | Computational Geometry in C[END_REF] and the rectangular tessellation, using the Mackisack algorithm [START_REF] Mackisack | Homogeneous rectangular tessellations[END_REF]. The land-use mosaic is built randomly according to a chosen distribution of land-uses. GenExP-LandSiTes was also coupled with the statistics R software to have high-level tools for the spatial point-process simulation [START_REF] Baddeley | Spatstat: an R package for analyzing spatial point patterns[END_REF]. Finally, GenExP-LandSiTes provides a library to calculate basic landscape descriptors (field area, perimeter, number of vertices, centroid and shape).

Application to the study of gene flow

In the case of coexisting GM, conventional and organic crops, two agronomic aspects have to be taken into account. The first aspect concerns the cropping of several agricultural species bearing the same transgene (for instance a tolerance to the same herbicide). Interactions between these different crops in a cropping system must be studied to define the best management of volunteers and avoid selection or breakdown of this particular gene.

The second aspect concerns the cultivation of GM and non-GM varieties of the same species for which, in Europe, coexistence rules between GM and non-GM crops (e.g., minimum distances between crops) have to be implemented to respect harvest impurity thresholds of non-GM crops, e.g., less than 0.9% of GM seeds in non-GM harvests (2003/556/EC, 1830/2003/EC). Several coexistence studies have been carried out to address this issue (see for instance: [START_REF] Bock | Scenarios for co-existence of genetically modified, conventional and organic crops in european agriculture[END_REF][START_REF] Tolstrup | Report from the Danish working group on the coexistence of genetically modified crops with conventional and organic crops[END_REF][START_REF] Sanvido | Konzept für die Koexistenz verschiedener landwirtschaftlicher Anbausysteme mit und ohne Gentechnik in der Schweiz[END_REF][START_REF] Messéan | New case studies on the coexistence of gm and crops in european agriculture[END_REF]). Few of them have dealt with the impact of the distribution of GM and non-GM crops in the landscape [START_REF] Messéan | New case studies on the coexistence of gm and crops in european agriculture[END_REF][START_REF] Lécroart | Coexistence of gm and non-gm maize: effect of regional structural variables on gm dissemination risk[END_REF]. Moreover, they were based on real agricultural landscape which limits the number of studied cases (due to the accessibility of spatialised data). Only a few recent studies were based on simulated landscapes [START_REF] Ceddia | Landscape gene flow, coexistence and threshold effect: The case of genetically modified herbicide tolerant oilseed rape (Brassicus napus)[END_REF][START_REF] Lavigne | How do genetically modified (GM) crops contribute to background levels of GM pollen in an agricultural landscape[END_REF][START_REF] Colbach | A simulation study of the medium-term effects of field patterns on cross-pollination rates in oilseed rape (brassica napus l.)[END_REF].

Combining landscape models with gene flow models

The GenExP-LandSiTes software was combined with Mapod-maize-a spatially explicit pollen dispersal model [START_REF] Angevin | Modelling impacts of cropping sys[END_REF]-and we assessed the influence of landscape structures on maize gene flow. The investigation entailed the following four-step process (figure 2).

Generation of field patterns

Three French landscapes were chosen for this study (Figure 3). The landscapes, restricted to a 1.5 km×1.5 km window, differed mainly in mean R (spatial statistics) original landscapes For each region, the pairwise-interaction model (section 2.1.1) was fitted to the pattern of field centroids (see [START_REF] Adamczyk | GenExP, un logiciel simulateur de paysages agricoles pour l'étude de la diffusion de transgènes[END_REF] for details). The parameters of the pairwise-interaction model were estimated using the maximum pseudolikelihood method proposed by [START_REF] Baddeley | Practical maximum pseudolikelihood for spatial point patterns[END_REF]. The seeds of Voronoi tessellations were generated according to the fitted pairwise-interaction process 2 whereas the seeds of rectangular tessellations were generated as independent unifor-mally distributed random points. 3 From each original landscape, ten virtual field patterns were designed using the GenExP-LandSiTes software, five based on the Voronoi tessellation and five based on the rectangular tessellation. Overall, thirty virtual field patterns were generated. 4 Taking into account the original patterns, there were 33 field patterns. Landscapes were simulated within a 1.5 km×1.5 km window without any predefined hole (no non-agricultural zone, see Figure 1).

Mapod-maize maps of GM impurities

Crop allocation to fields

The total proportion of field area allocated to maize was set at 70%, thus simulating French production areas where maize is a major crop. The proportion of GM maize within this area was set to either 10% or 50%, simulating low and high acceptance of GM maize, respectively. GM and non-GM maize were randomly attributed to the fields in each of the 1.5 x 1.5 maps until the desired proportions were met. Three replicated allocations of crops were performed per map. Overall, we thus produced 198 crop allocations (33 field patterns x 2 proportions of GM maize x 3 replicates).

Pollen dispersal

Mapod-maize simulates pollen dispersal through a quasi-mechanistic approach that takes into account crop phenology and various agronomic and climatic conditions [START_REF] Angevin | Modelling impacts of cropping sys[END_REF]. The sub-model for pollen dispersal was adapted from [START_REF] Klein | Corn pollen dispersal: quasi-mechanistic models and field experiments[END_REF] and it accounts for distance between source and receptor plants as well as wind direction and intensity, and as plant height. To focus on the impact of spatial characteristics of landscapes, we kept agronomic and climatic inputs constant and identical for the GM and the conventional varieties. In particular, both varieties were assumed homozygous with the same pollen production and synchronous flowering. Wind was supposed moderate (3m/s) and orientated from west to east to simulate common French situations.

GM impurity rates

GM impurity rates were calculated in each non-GM field as the proportion of seeds containing at least one copy of the transgene at harvest. For each simulated crop allocation, we also calculated the proportion of fields with impurity rates exceeding either 0.1% or 1%.

Statistical analyses

3.2.1 Assessing the impact of the tessellation method on landscapes characteristics [START_REF] Riitters | A factor analysis of landscape pattern and structure metrics[END_REF] showed that six indices are sufficient to characterize landscape patterns. Since we were not interested in fractal landscapes, land-use or scaling, only two remained: the average patch compaction (i.e., perimeter/area), and the average patch shape (perimeter/ √ area). In the case of simple convex shapes, as produced by the Voronoi and rectangular tessellations, the average patch compaction is highly correlated to the average patch area. The areas (in hectares) and the shapes (I S = perimeter/4

√ area) of fields were thus calculated over each of the 33 field patterns. To assess the within-landscape variability of field areas and shapes, we also calculated the standard deviation of these measurements within field patterns. Although statistical tests on simulated data should be carefully considered, means and standard deviations of area and shape index at the field pattern level were submitted to analyses of variance using original landscape (qualitative, three levels), tessellation method (three levels: none=original landscapes, Voronoi, and rectangular) and their interaction as factors. Field areas and standard deviations of the shape index within field patterns were log-transformed to improve normality of residuals.

Besides the geometry of the simulated fields, we were also interested in their spatial distribution. We evaluated this distribution by examining the distances between field centroids, using the nearest neighbour distance distribution function G. For a given distance r the value G(r) gives the probability that the nearest-neighbour of an arbitrary centroid lies within distance.

Assessing the impact of tessellation method on GM impurity rates

Impacts on GM impurity rates for landscape type and tessellation method were assessed through analyses of variance and comparisons of least-squares means. Fixed effect used in the model were original landscape, tessellation method, proportion of GM crop (qualitative, two levels) and their two-way interactions. To account for correlations among the five simulations performed with each tessellation for each original landscape, a repet factor was declared as random and nested in (region, tessellation, GM proportion).

Similarly, to account for correlations among fields within a single crop pattern, allocation nested in (repet, region, tessellation, GM proportion) was also declared as random. All GM impurity rates were log-transformed prior to analyses to improve the normality of residuals. The same model was applied to the proportions of fields exceeding the 0.1% and 1% threshold except that the allocation term was not needed because proportions are response variables defined at the landscape level. Prior to analyses, square roots of proportions were arcsine-transformed to homogenise variances.

Results

Comparison of original and simulated field patterns

The number of fields (Table 1) was fixed for the rectangular tessellation, since we used a fixed number of seeds. On the contrary, for the Voronoi tessellation, the number of seeds is a random variable whose mean value should be close to the number of observed centroids: in our case it appeared to be slightly lower. One explanation is that the Voronoi tessellation produced infinite cells on the landscape borders that were eliminated (cf. Figure 1). areas and shapes were not affected by the same factors: whereas field areas and their within-landscape variability were mostly explained by the original landscape, field shapes and their variability were mostly explained by the tessellation method. Variability due to small random variations of seed distribution among replicate patterns was small (variance estimates in Table 2). Average field areas (Table 1) with the rectangular tessellations were close to the ratio landscape area/number of original fields. Actually, the rectangular tessellation is a total covering of the initial window. For the same number of fields, the original landscapes had a lower average field area since they contained non-agricultural areas. The results of the Voronoi tessellations were slightly more variable. In contrast, the variability of field areas was lower for the Voronoi tessellations than for both real landscapes and rectangular tessellations. Rectangular tessellations had the highest variability. Regarding field shapes (Table 1), as expected, each tessellation gave homogeneous results with respect to the three original landscapes, but somewhat different from the the original field shapes. Indeed, the shape index of rectangular fields was too high (elongated shapes) while the Voronoi fields were too compact. Furthermore, the within-landscape variability of shapes was very high for the rectangular tessellation and somewhat low for the Voronoi tessellation.

Regarding the spatial distribution of fields, results for region A1 are given in Figure 4. The estimates of the G function for the Voronoi simulations were close to the original pattern. Indeed, the seed generation in Voronoi 1: For each original landscape (P1, A1, S4) and each tessellation method: average number of fields, field areas, and field shapes (mean ± standard error); variabilities of field areas and field shapes measured as within-landscape standard deviation (mean ± standard error).

simulations was based on a pairwise-interaction process fitted to the real centroids, and thanks to the tessellation method, the centroids of tessellation cells are close to the seeds. On the contrary, the estimates of the G function for rectangular tessellation differed from the original pattern. Indeed, the rectangular tessellation did not allow us to link the centroids of the simulated rectangles to the centroids of the actual landscapes.

Impact of the tessellation method on GM impurity rates

As expected, the proportions of GM maize over the landscape had the highest effect on impurity rates (Tables 3 and4). The impact of the tessellation method was comparable to that of the original landscape (close F-values) when either individual field cross-pollination rate or proportion of fields above the 0.1% threshold were considered. Differences between the tessellation methods decreased largely when the 1% threshold value was considered and became much smaller than differences between original landscapes (Table 4). The impact of the tessellation method did not differ much between regions (small F-value for the Tessellation x Region interactions).

The large effect of the tessellation method was mainly due to differences between Voronoi and rectangular simulated landscapes. Values for real landscapes were generally close to those on simulated Voronoi landscapes: pairwise comparisons of least-squares means indicate that Voronoi tessellations resulted in impurity rates closer to those simulated on original landscapes than rectangular tessellations. 5 The same result was observed in terms of the proportions of fields exceeding the 0.1% threshold. 6 As expected from the small F-value, outputs of simulations were very close for the two tessellation methods and the original landscape when the 1% threshold was considered. 7

Discussion

Methods have been investigated for modelling 2D-geometrical landscapes, such as a distribution of fields over an area, maintaining certain characteristics of the original landscapes. The rationale for doing so lies in the recent increasing interest in spatially explicit modelling that, to some extent, is due to the increased availability of land-use maps through aerial or satellite pictures. Relating landscape characteristics (e.g., hedgerow density) to ecological (e.g., spread of an invasive species) or physical (e.g., water pollution) processes has been increasingly frequent. However, understanding how landscape structures influence the observed spatial patterns requires more than data from a few maps, because each map is unique and differs from others in many ways. Coupling landscape models to process-explicit models (such as the here-used model for pollen dispersal) provides a means to explore the impact on a given process of small random variations in landscape characteristics (replication) as well as large systematic differences of these characteristics.

Whereas neutral landscape models are frequent in ecology [START_REF]Quantitative Methods in Landscape Ecology[END_REF][START_REF] With | The use and misuse of neutral landscape models in ecology[END_REF], tessellation methods have been rarely used. Furthermore, we know of no other landscape model exploring variation in the geometry of fields except that of [START_REF] Gaucherel | Neutral models for polygonal landscapes with linear networks[END_REF]. Such an approach means that the impact of the geometry of agricultural landscapes on agro-ecological processes can be investigated. It also allows crop allocations over similar (but not identical) geometrical field patterns to be varied. The difficulty is in controlling important landscape characteristics, here for example field areas and shapes or distances among field centroids, while allowing other characteristics to vary.

Comparing tessellation methods w.r.t. original landscapes

Neither simulation method really succeeded in providing replicates of real landscapes that would allow studying the impact of small random variations in landscape geometry. While both tessellation methods allowed simulation of landscapes with numbers of fields and average field areas that resembled those of the original landscapes, neither correctly reflected the shape and most of all the within-landscape variability of this shape: the Voronoi tessellation produced somewhat too compact fields with low variability while the rectangular tessellation produced overly elongated fields with an overly high variability. Systematic post-processing, e.g., deleting too small/elongated cells (in the case of the rectangular tessellation) or deleting vertices in order to deform shapes (in the case of the Voronoi tessellation) could be a solution to improve suitability of tessellation results with the original landscapes.

As mentioned in Section 2.1, the a priori major advantage of the Voronoi tessellation over the rectangular tessellation for simulating agricultural landscapes is that a correspondence can be established between the pair polygon/seed and the pair field/centroid. This property is not exact, especially for concave fields, but it proved efficient to simulate patterns where field centroids respect certain constraints, e.g., a minimum distance between two centroids. Indeed, using the pairwise-interaction point-process (see Section 2.1.1) we produced Voronoi tessellations with similar distributions of nearest-neighbour distances between centroids as the original landscapes. Other characteristics of landscapes have to be investigated. For example, in the rectangular tessellation, we could take advantage of the growing directions of the segments in order to simulate the average orientations of fields.

However, both the Voronoi and the rectangular tessellations start with the distribution of a set of points. This step is completely artificial in the context of real landscape formation and makes a model fit rather difficult. One possible improvement consists on studying the tessellation models based on the geometrical units naturally present in an agricultural landscape: polygons or lines. In a recent work, [START_REF] Gaucherel | Neutral models for polygonal landscapes with linear networks[END_REF] proposed to use the Delaunay triangles (rather than the Voronoi polygons) and to aggregate them according to linear networks. Furthermore, line-based tessellations, such as proposed by [START_REF] Nagel | The crack tessellation -characterization of the stationary random tessellations which are stable with respect to iteration[END_REF], can also be explored to simulate the geometry of agricultural landscapes.

Comparing pollen flows w.r.t the characteristics of landscapes

A recent spatial sensitivity of GM impurity rates that also used the Mapod-maize software has concluded that field geometry was a minor factor affecting impurity rates in comparison to maize frequency and spatial ar-rangement [START_REF] Viaud | Spatial sensivity of maize gene-flow from genetically modified to conventional varieties to landscape pattern: a simulation approach[END_REF]. Our results confirm this: GM maize proportion was the main factor explaining impurity rates. However, our simulations highlighted some influence of field geometry: GM impurity rates obtained using the Voronoi tessellations tended to be a little under-estimated, while those from the rectangular tessellation tended to be over-estimated, in particular more fields exhibited impurity rates that exceeded 0.1%. This difference might be explained both by the larger variability in field area as well as by the elongated field shapes in the rectangular tessellation. Indeed, in simulations that only took into account one source and one recipient field, large compact recipient fields such as those simulated by the Voronoi tessellation have been shown to be the least susceptible to incoming gene flow [START_REF] Damgaard | Gene flow of oilseed rape (Brassica napus) according to isolation distance and buffer zone[END_REF][START_REF] Klein | Pollen dispersal of oilseed rape: estimation of the dispersal function and effects of field dimension[END_REF]. This property is likely due to a larger average distance of points to field edges, distance being maximum within Voronoi polygons. general, impurity rates obtained with Voronoi tessellations were close to those simulated on real landscapes. This is to some extent due to the choice of real landscapes, two of them being covered by regular and mostly large fields, but it is also likely to be a general result as small elongated fields are less and less frequent in current landscapes due to mechanisation and land consolidation. On the contrary, landscapes with regular compact fields are more and more the case. As a first step, Voronoi tessellation thus appears a good basis for simulating landscapes for empirical gene flow such as Mapod-maize.

Simulating land-use

In this paper, we only simulated three land-uses (conventional maize, GM maize, others) that were allocated randomly within the landscapes. Characteristics of the fields, such as the area or the neighbourood, were not taken into account. Furthermore, we only needed a 1-year crop simulation, since maize pollen or seeds do not persist under European climatic conditions. To investigate the effect of landscape structure and land-use on other crops such as rapeseed [START_REF] Colbach | How to model and simulate the effects of cropping systems on population dynamics and gene flow at the landscape level. Example of oilseed rape volunteers and their role for co-existence of GM and non-GM crops[END_REF], temporal (succession of land-uses) and spatial (surface or neighbouring crops) constraints should be taken into account. This is possible thanks to the underway coupling of GenExP-LandSiTes with the data-mining software CarrotAge [START_REF] Ber | Studying crop sequences with CarrotAge, a HMM-based data mining software[END_REF].

Conclusion

When studying agro-ecological or physical processes such as gene flow, erosion, and spore dispersal, in an agricultural landscape, a large number of fields and their spatial organisation must be considered, and therefore, experimentation under field conditions is complicated and expensive, if not impossible. Computer simulations are less costly and dramatically increase the number of configurations that can be investigated. Furthermore, there is still little landscape data available, and the extrapolation of the results obtained in a few specific contexts is not easy. Thus, it is necessary to simulate neutral agricultural landscapes in order to explore the variability of landscape characteristics from which more general conclusions can be drawn. In this paper, we tested two tessellation methods with respect to gene flow. Granted, neither of the methods fully captured field shapes of real landscapes and therefore they only partly reached the given objectives. Indeed, two methods either over-estimated (for the rectangular tessellation) or (slightly) gene flow (for the Voronoi tessellation). However, the results have provided promising leads for further developments. We shall have to experiment other process models and improve tessellation methods for a better control of their behaviour with respect to the characteristics of actual landscapes. Possible research directions should include the use of post-processes (i.e., deletion of small cells, fusion of vertices), the control of cell orientations (rectangular tessellation), or the implementation of other tessellation techniques (e.g., line-based tessellations).

Figure 1 :

 1 Figure 1: Two different tessellations based on sets of seeds (region A1, scale 1/5): the Voronoi tessellation (a) and the rectangular tessellation with only T-vertices (b). Infinite Voronoi cells were eliminated.

Figure 2 :

 2 Figure 2: Combination process of the GenExP-LandSiTes and Mapodmaize softwares: (1) generation of field patterns, (2) crop allocation, (3) pollen dispersal, and (4) GM impurity rates.

Figure 3 :

 3 Figure 3: Three agricultural landscapes with plot centroids.

Figure 4 :

 4 Figure 4: Estimate of G function for the field centroids of pattern A1. Red line: original pattern. Gray lines: five Voronoi simulations. Blue lines: five rectangular simulations. X-axis: distance (m) between two nearest neighbours.

  tems and climate on maize cross pollination in agricultural landscapes: The MAPOD model. European Journal of Agronomy, 28(3):471-484, apr 2008.

  Neutral modelling of agricultural landscapes by tessellation methods -Application for gene flow simulation Texte paru dans Ecological Modelling, doi:10.1016/j.ecolmodel.2009.06.019 F. Le Ber * , C. Lavigne † , K. Adamczyk ‡ , F. Angevin § , N. Colbach ¶ , J.-F. Mari , and H. Monod * *

Table 2 :

 2 Analysis of variance on characteristics of original and simulated

	Source	DF	Type 3 Sum Sq (F)	Type 3 Sum Sq (F)
				Area	St dev Area
	Region	2 4.129	(90.79 *** ) 8.91	(92.69 *** )
	Tessellation	2 0.028	(2.67 ns )	7.95	( 82.67 *** )
	Tess. x Region	4 0.072	(3.40 * )	1.05	(5.45 ** )
	error	24 0.983		1.154
				Shape	St dev Shape
	Region	2 4.77 10 -3 (0.21 ns )	0.11	(2.03 ns )
	Tessellation	2 1.20	(53.27 *** ) 33.99 (642.6 *** )
	Tess. x region	4 2.38 10 -3 (0.05 ns )	0.21	(1.96 ns )
	error	24 0.272		0.635

field patterns (DF: degrees of freedom; Sum Sq: sum of squares; F: Fisher statistics; St dev: standard deviation; ns : P > 0.05; * : P > 0.01; ** : P > 0.001; *** : P < 0.001). The F tests were performed assuming the mixed model described in the main text.

Table 3 :

 3 Analysis of variance with fixed and random effects (Proc Mixed in SAS 8.01, SAS Institute, Cary, NC, USA) on individual field (log 10 ) GM impurity rates (DF: degrees of freedom; F: Fisher statistics; St dev: standard deviation). All Fisher statistics are calculated using the repet in (tess. x region x GM prop.) variance as error variance (52 degrees of freedom).

	Effect	DF	F
	Tessellation	2	7.20
	Region	2	8.53
	GM proportion	1 319.46
	Tessellation x region	4	1.54
	Tessellation x GM prop.	2	2.30
	Region x GM prop.	2	0.47
	Error variances	Estimation St dev
	repet in (tess. x region x GM prop.)	0.01	0.04
	alloc in (repet. x tess. x region x GM prop.)	0.38	0.06
	Residual	4.20	0.06

Table 4 :

 4 Analysis of variance with fixed and random effects (Proc Mixed in SAS 8.01, SAS Institute, Cary, NC, USA) on proportions of fields above a given threshold of impurity rate (DF: degrees of freedom; F: Fisher statistics; dev: standard deviation). All Fisher statistics are calculated using the repet in (tess. x region x GM prop.) variance as error variance (52 degrees of freedom).

LANDscape SImulation based on TESsellation.

We used the Metropolis-Hasting algorithm, implemented in the R-spatstat library[START_REF] Baddeley | Spatstat: an R package for analyzing spatial point patterns[END_REF].

Number of seeds = number of fields in the original pattern minus 1.

3 original patterns ×

simulations × 2 tessellation methods.

LS mean± sd log10(impurity rate): original = -3.55 ± 0.14; Voronoi = -3.51 ± 0.06; rectangular = -3.18 ± 0.06

LS mean± sd of arc sinus(proportion above 0.1%) 1/2 : original = 0.88 ± 0.04; Voronoi = 0.90 ± 0.017; rectangular = 0.98 ± 0.017

LSmean± sd arc sinus (proportion above 1%) 1/2 : original = 0.52 ± 0.02; Voronoi = 0.49 ± 0.01; rectangular = 0.51 ± 0.01

Acknowledgments

Original maps were provided by the Institute for the Protection and Security of the Citizen (Joint Research Centre of the European Union) and AUP (Agence Unique de Paiement / French Payment Agency CAP Support). We acknowledge support from the French national program "ACI OGM et environnement". We wish to thank Ms S. Tanis-Plant for fruitful discussion and advice in English.