Stochastic order and attractiveness for particle systems with multiple births, deaths and jumps

An approach to analyse the properties of a particle system is to compare it with different processes to understand when one of them is larger than other ones. The main technique for that is coupling, which may not be easy to construct. We give a characterization of stochastic order between different interacting particle systems in a large class of processes with births, deaths and jumps of many particles per time depending on the configuration in a general way: it consists in checking inequalities involving the transition rates. We construct explicitly the coupling that characterizes the stochastic order. As a corollary we get necessary and sufficient conditions for attractiveness. As an application, we first give the conditions on examples including reaction-diffusion processes, multitype contact process and conservative dynamics and then we improve an ergodicity result for an epidemic model.


Introduction
The use of interacting particle systems to study biological models is becoming more and more fruitful. In many biological applications a particle represents an individual from a species which interacts with others in many different ways; the empty configuration is often an absorbing state and corresponds to the extinction of that species. An important problem is to find conditions which give either the survival of the species, or the almost sure extinction. When the population in a system is always larger (or smaller) than the number of individuals of another one there is a stochastic order between the two processes and one can get information on the larger population starting from the smaller one and vice-versa.
Attractiveness is a property concerning the distribution at time t of two processes with the same generator : if a process is attractive the stochastic order between two processes starting from different configurations is preserved in the time evolution (see Section 2.1).
The main technique to check if there is stochastic order between two systems is coupling: if the transitions are intricate an increasing coupling may be hard to construct. The main result of the paper (Theorem 2.4, Section 2.1) gives a characterization of the stochastic order (resp. attractiveness) in a large class of interacting particle systems: in order to verify if two particle systems are stochastically ordered (resp. one particle system is attractive), we are reduced to check inequalities involving the transition rates.
A first motivation is a general understanding of the ordering conditions between two processes. The analysis of interacting particle systems began with Spin Systems, that are processes with state space {0, 1} Z d . We refer to [11] and [12] for construction and main results. The most famous examples are Ising model, contact process and voter model. These processes have been largely investigated, in particular their attractiveness (see [11, Chapter III, Section 2]). Many other models taking place on X Z d , where X = {0, 1, . . . M } ⊆ N, that is with more than one particle per site, have been studied. Reaction-diffusion processes, for example, are processes with state space N Z d (hence non compact), used to model chemical reactions. We refer to [4,Chapter 13] for a general introduction and construction. In such particle systems a birth, death or jump of at most one particle per time is allowed. But sometimes the model requires births or deaths of more than one particle per time. This is the case of biological systems with mass extinction ( [17], [18]), or multitype contact process ( [6], [7], [14], [16]). A partial understanding of attractiveness properties can be found in [20] for the multitype contact process. A system with jumps of many particles per time has been investigated in [8,Theorem 2.21], where the authors found necessary and sufficient conditions for attractiveness for a conservative particle system with multiple jumps on N Z d with misanthrope type rates.
Those examples and the need for more realistic models for metapopulation dynamics systems ( [9]) has led us to consider systems ruled by births, deaths and migrations of more than one individual per time with general transition rates to get an exhaustive analysis of the stochastic order behaviour and attractiveness. Our method relies on [8], that it generalizes.
The main applications allow to investigate the ergodic properties of a process. A process is ergodic if there exists a unique invariant measure to which it converges starting from any initial configuration: if the process is attractive, it is enough to check the convergence of the largest and the smallest initial configurations. This is a first application of Theorem 2.4. In Section 2.2.5 we combine attractiveness and a technique called u−criterion (see [4]) to get ergodicity conditions on a model for spread of epidemics, either if there is a trivial invariant measure or not.
For many biological models the empty configuration 0 is an absorbing state and the main question is if the particle system may survive, that is if there is a positive probability that the process does not converge to the Dirac measure δ 0 concentrated on 0, which is a trivial invariant measure. In order to prove that a metapopulation dynamics model (see [3]) survives, we largely use comparison (therefore stochastic order) with auxiliary processes: this is a second application of the result. Instead of constructing a different coupling for each comparison, we just check that inequalities of Theorem 2.4 are satisfied on the transition rates. Moreover the main technique we use to get survival is a comparison with oriented percolation (see [6]), and attractiveness is a key tool in many steps of the proofs.
The survival of a process does not imply the existence of a non trivial invariant measure: one can have the presence of particles in the system for all times but no invariant measures. If the process is attractive and the state space is compact, a standard approach allows to construct such a measure starting from the largest initial configuration: this is the third application. Once we get survival, we use this argument to construct non trivial invariant measures for metapopulation dynamics models. In Section 2.2.4 we introduce a metapopulation dynamics model with mass migration and Allee effect investigated in [3].
The transition rates of the particle systems we analyse in this paper depend on two sites x, y, on the number of particles at x and y and on the number of particles k involved in a transition: they are of the form b(k, α, β)p(x, y), where α and β are respectively the number of particles on x and y and p(x, y) is a probability distribution on Z d given by a bistochastic matrix (we require neither symmetry nor translation invariance). Moreover we allow birth and death rates on a site x depending only on the configuration state on x.
In other words we work with three different types of transition rates: given a configuration η, on each site y we can have a birth (death) of k individuals depending on the configuration state η(y) on the same site y with rate P k η(y) (P −k η(y) ) and depending also on the number of particles on the other sites x = y with rate x R 0,k η(x),η(y) p(x, y) ( x R −k,0 η(y),η(x) p(y, x)). We consider a death rate R −k,0 η(y),η(x) p(y, x) instead of the more natural R −k,0 η(y),η(x) p(x, y) to simplify the proofs and since we are interested in applications given by a symmetric probability distribution p(·, ·). This represents a possible different interaction rule between individuals of the same population and individuals from different populations. We can have a jump of k particles from x to y with rate Γ k η(x),η(y) p(x, y), which represents a migration of a flock of individuals (see Section 2.1). We require that the birth/death and jump rates differ only on the term b(k, α, β), that is the conservative and non-conservative rates depend on the same probability distribution p(x, y).
In Section 2.1 we recall some classical definitions and propositions needed in the sequel, we introduce the particle system with more details and we state the main result, Theorem 2.4. In Section 2.2 we derive the conditions on several examples (multitype contact processes, conservative dynamics and reaction-diffusion processes); we also detail the conditions on models with transitions of at most one particle per time. In Section 2.2.5 we apply the attractiveness conditions and the so-called u-criterion technique to improve an ergodicity result for a model of spread of epidemics. Others applications to the construction of non-trivial invariant measures in metapopulation dynamics models will be presented in a subsequent paper (see [3]). In Section 3 we prove Theorem 2.4: the coupling is constructed explicitly through a downwards recursive formula in Section 3.2, where a detailed analysis of the coupling mechanisms is presented. We have to mix births, jumps and deaths in a non trivial way by following a preferential direction. Section 4 is devoted to the proofs needed for the application to the epidemic model. Finally we propose some possible extensions to more general systems.

Main result and applications 2.1 Stochastic order an attractiveness
Denote by S = Z d the set of sites and let X ⊆ N be the set of possible states on each site of an interacting particle system (η t ) t≥0 on the state space Ω = X S , with semi-group T (t) and infinitesimal generator L given, for a local function f , by where χ x,y α,β is the indicator of configurations with values (α, β) on (x, y), that is and S −k,k x,y , S k y , S −k y , where k > 0, are local operators performing the transformations whenever possible The transition rates have the following meaning: -p(x, y) is a bistochastic probability distribution on Z d ; -Γ k α,β p(x, y) is the jump rate of k particles from x, where η(x) = α, to y, where η(y) = β; -R 0,k α,β p(x, y) is the part of the birth rate of k particles in y such that η(y) = β which depends on the value of η in x (that is α); -R −k,0 α,β p(x, y) is the part of the death rate of k particles in x such that η(x) = α which depends on the value of η in y (that is β); -P ±k β is the birth/death rate of k particles in η(y) = β which depends only on the value of η in y: we call it an independent birth/death rate.
We call addition on y (subtraction from x) of k particles the birth on y (death on x) or jump from x to y of k particles. By convention we take births on the right subscript and deaths on the left one: formula (2.1) involves births upon β, deaths from α and a fixed direction, from α to β, for jumps of particles. We define, for notational convenience We refer to [11] for the classical construction in a compact state space. Since we are interested also in non compact cases, we assume that (η t ) t≥0 is a well defined Markov process on a subset Ω 0 ⊂ Ω, and for any bounded local function f on Ω 0 , We will be more precise on the induced conditions on transition rates in the examples. We state here only a common necessary condition on the rates.
In other words, for each α, β there exists a maximal number of particles involved in birth, death and jump rates. Notice that N (α, β) is not necessarily equal to N (β, α), which involves deaths from β, births upon α and jumps from β to α.
The particle system admits an invariant measure µ if µ is such that P µ (η t ∈ A) = µ(A) for each t ≥ 0, A ⊆ Ω, where P µ is the law of the process starting from the initial distribution µ. An invariant measure is trivial if it is concentrated on an absorbing state when there exists any. The process is ergodic if there is a unique invariant measure to which the process converges starting from any initial distribution (see [11,Definition 1.9]). Given two processes (ξ t ) t≥0 and (ζ t ) t≥0 , a coupled process (ξ t , ζ t ) t≥0 is Markovian with state space Ω 0 × Ω 0 , and such that each marginal is a copy of the original process. We define a partial order on the state space: We denote by M the set of all bounded, monotone continuous functions on Ω. The partial order on Ω induces a stochastic order on the set P of probability measures on Ω endowed with the weak topology: The following theorem is a key result to compare distributions of processes with different generators starting with different initial distributions.
Theorem 2.2 Let (ξ t ) t≥0 and (ζ t ) t≥0 be two processes with generators L and L and semigroups T (t) and T (t) respectively. The following two statements are equivalent: The proof is a slight modification of [11, proof of Theorem II.2.2].
Definition 2.3 A process (ζ t ) t≥0 is stochastically larger than a process (ξ t ) t≥0 if the equivalent conditions of Theorem 2.2 are satisfied. In this case the process (ξ t ) t≥0 is stochastically smaller than (ζ t ) t≥0 and the pair (ξ t , ζ t ) t≥0 is stochastically ordered.
Attractiveness is a property concerning the distribution at time t of two processes with the same generator which start with different initial distributions. By taking T = T , Theorem 2.2 reduces to [11,Theorem II.2.2] and Definition 2.3 is equivalent to the definition of an attractive process (see [11,Definition II.2.3]). If an attractive process in a compact state space starts from the larger initial configuration, it converges to an invariant measure.
Remark 2.5 The restriction K ≤ N (α, β) ∨ N (γ, δ) avoids that an infinite number of K, I a , I b , I c , I d result in the same rate inequality. Since Γ k α,β = 0 for each k > N (α, β), if K > N (α, β) no terms are being added to the left hand side of (2.13), and adding more terms on the right hand side does not give any new restrictions. A similar statement holds for (2.14), with the corresponding condition K ≤ N (γ, δ). We follow the approach in [8]: in order to characterize the stochastic ordering of two processes, first of all we find necessary conditions on the transition rates. Then we construct a Markovian increasing coupling, that is a coupled process (ξ t , ζ t ) t≥0 which has the property that ξ 0 ≤ ζ 0 implies for all t ≥ 0. Here P (ξ 0 ,ζ 0 ) denotes the distribution of (ξ t , ζ t ) t≥0 with initial state (ξ 0 , ζ 0 ).

No multiple births, deaths or jumps
Proposition 2.7 If N = 1 then a change of at most one particle per time is allowed and Conditions (2.13) and (2.14) become (2.9). Since N = 1 the left hand side of (2.13) is null; if β = δ the only case for which the left hand side of (2.13) is not null is j 1 = 0, which gives Since N = 1, the value K = 1 covers all possible sets I a and I b , namely I a = {k : m 1 ≥ k > 0} and I b = {γ − α + m 1 ≥ l > 0}. If m 1 > 0, we get (2.15). If γ = α and m 1 = 0 we get (2.16). One can prove (2.17) in a similar way.
If β = δ and γ ≥ α, Formula (2.15) expresses that the sum of the addition rates of the smaller process on y in state β must be smaller than the corresponding addition rates on y of the larger process on y in the same state. If β = δ and γ = α we also need that the birth rate of the smaller process on y is smaller than the one of the larger process, that is (2.16). Conditions (2.17)-(2.18) have a symmetric meaning with respect to subtraction of particles from x. Proposition 2.7 will be used in a companion paper for metapopulation models, see [3]. If R 0,k α,β = 0 for all α, β, k, the model is the reaction diffusion process studied by Chen (see [4]) and the attractiveness Conditions (2.15), (2.17) (the only ones by Remark 2.8) reduce to In other words we need Γ 1 α,β to be non decreasing with respect to α for each fixed β, and non increasing with respect to β for each fixed α. In [4], the author introduces several couplings in order to find ergodicity conditions of reaction diffusion processes. All these couplings are identical to the coupling H introduced in Section 3.2 (and detailed in Appendix A if N = 1), on configurations where an addition or a subtraction of particles may break the partial order, but differ from H on configurations where it cannot happen.
Spread of tubercolosis model ( [17]). Here M represents the number of individuals in a population at a site x ∈ Z d . The transitions are: where y ∼ x is one of the 2d nearest neighbours of site x. Given two systems with parameters (λ, φ, M ) and (λ, φ, M ), the proof of [ In the following examples we suppose S = S, that is we consider necessary and sufficient conditions for attractiveness.
2-type contact process ( [14]). In this model M = 2. Since a value on a given site does not represent the number of particles on that site, we write the state space {A, B, C} Z d . The value B represents the presence of a type-B species, C the presence of a type-C species and A an empty site. If A = 0, B = 1, C = 2 then the transitions are By taking β = 0, δ = 1, α = γ = 2 we get 2dλ 2 ≤ 0, which is not satisfied since λ 2 > 0. As already observed, see [20, Section 5.1], one can get an attractive process by changing the order between species: namely by taking A = 1, B = 0 and C = 2 the process is attractive.

Conservative dynamics
If Π 0,k α,β = 0, for all (α, β) ∈ X 2 , k ∈ N, we get a particular case of the model introduced in [8], for which neither particles births nor deaths are allowed, and the particle system is conservative. Suppose that in this model the rate Γ k α,β (y − x) has the form Γ k α,β p(y − x) for each k, α, β. Necessary and sufficient conditions for attractiveness are given by [8, Theorem 2.21]: Condition (2.23) follows in a similar way.

Metapopulation model with Allee effect and mass migration
The third model investigated in [3] is a metapopulation dynamics model where migrations of many individuals per time are allowed to avoid the biological phenomenon of the Allee effect (see [1], [19]). The state space is compact and on each site for each α, β ∈ X, and p(x, y) = 1 2d 1l {y∼x} . In other words each individual reproduces with rate 1, but dies with different rates: either φ A if the local population size is smaller than M A (Allee effect) or φ if it is larger. When a local population has more than M − N individuals a migration of more than one individual per time is allowed. Such a process is attractive by [3, Proposition 5.1, where N and M play opposite roles], which is an application of Theorem 2.4.

Individual recovery epidemic model
We apply Theorem 2.4 to get new ergodicity conditions for a model of spread of epidemics. The most investigated interacting particle system that models the spread of epidemics is the contact process, introduced by Harris [10]. It is a spin system (η t ) t≥0 on {0, 1} Z d ruled by the transitions See [11] and [12] for an exhaustive analysis of this model. In order to understand the role of social clusters in the spread of epidemics, Schinazi [17] introduced a generalization of the contact process. Then, Belhadji [2] investigated some generalizations of this model: on each site in Z d there is a cluster of M ≤ ∞ individuals, where each individual can be healthy or infected. A cluster is infected if there is at least one infected individual, otherwise it is healthy. The illness moves from an infected individual to a healthy one with rate φ if they are in the same cluster. The infection rate between different clusters is different: the epidemics moves from an infected individual in a cluster y to an individual in a neighboring cluster x with rate λ if x is healthy, and with rate β if x is infected.
We focus on one of those models, the individual recovery epidemic model in a compact state space: each sick individual recovers after an exponential time and each cluster contains at most M individuals. The non-null transition rates are and Γ k η(x),η(y) = 0 for all k ∈ N. The rate γ represents a positive "pure birth" of the illness: by setting γ = 0, we get the epidemic model in [2], where the author analyses the system with M < ∞ and M = ∞ and shows [2,Theorem 14] that different phase transitions occur with respect to λ and φ. Moreover ([2, Theorem 15] the disease dies out for each cluster size M . By using the attractiveness of the model we improve this ergodicity condition. Notice that a dependence on the cluster size M appears. with φ < 1 and either i) γ = 0, or ii) γ > 0 and β − λ ≤ γ/(2d).
Then the system is ergodic.
Notice that if γ > 0 and β ≤ λ hypothesis ii) is trivially satisfied. If M = 1 and γ = 0 the process reduces to the contact process and the result is a well known (and already improved) ergodic result (see for instance [11,Corollary 4.4, Chapter VI]); as a corollary we get the ergodicity result in the non compact case as M goes to infinity. In order to prove Theorem 2.10, we use a technique called u-criterion. It gives sufficient conditions on transition rates which yield ergodicity of an attractive translation invariant process. It has been used by several authors (see [4], [5], [13]) for reaction-diffusion processes.
First of all we observe that The process is attractive for all λ, β, γ, φ, M .
The key point for attractiveness is that R 0,1 η(x),η(y) is increasing in η(x).
Given ǫ > 0 and {u l (ǫ)} l∈X such that u l (ǫ) > 0 for all l ∈ X, let F ǫ : X × X → R + be defined by for all x, y ∈ X. When not necessary we omit the dependence on ǫ and we simply write this is a metric on X and it induces in a natural way a metric on Ω. Namely, for each η and ξ in Ω we define The key idea consists in taking a "good sequence" {u l } l∈X and in looking for conditions on the rates under which the expected value E(·) (with respect to a coupled measure P) of the distance between η M t and η 0 t converges to zero as t goes to infinity, uniformly with respect to x ∈ S. We will use the generator properties and Gronwall's Lemma to prove that if there exists ǫ > 0 and a sequence {u l (ǫ)} l∈X , u l (ǫ) > 0 for all l ∈ X such that the metric F satisfies uniformly with respect to x ∈ S so that the distance ρ α (·, ·) between the larger and the smaller process converges to zero, and ergodicity follows.
Hence we are left with checking the existence of ǫ > 0 and positive {u l (ǫ)} l∈X which satisfy (2.35). Such a choice is not unique.

Remark 2.13 Given
Another possible choice is Definition 2.14 Given ǫ > 0 and U > 0, we set u 0 (ǫ) = U and we define (u l (ǫ)) l∈X recursively through Definition 2.14 gives a better choice of {u l (ǫ)} l∈X , indeed the u-criterion is satisfied under the more general assumption (2.29). Proofs of Theorems 2.10 and 2.12 are detailed in Section 4.

Coupling construction and proof of Theorem 2.4
In this section we prove the main result: we begin with the necessary condition, based on [8, Proposition 2.24].
which gives, by Assumption (2.6), We have, by using (2.1), We write (L1l V )(η) by using the corresponding rates of S.
We define the set C + y of sites which interact with y with an increase of the configuration on y, Denote by x = (x 1 , . . . , x d ) the coordinates of each x ∈ S. We define, for each n ∈ N, We may suppose C + y = ∅, since otherwise (2.13)-(2.14) would be trivially satisfied. Given K ∈ N, we fix {p i z } i≤K,z∈C + y such that for each i, z, p i z ∈ X and p i z ≤ ξ(z). Moreover we fix {p i y } i≤K such that p i y > δ, for each i. For n ∈ N, let The union of increasing cylinder sets I y (n) is an increasing set, to which neither ξ nor η belong. We compute, using (3.3), So, by setting Taking the monotone limit n → ∞ gives z∈S k∈J(pz,αz,β) to the complement of which ξ and η belong, and the application of inequality (3.2) to ξ, η, Ω\D x (n) (which is an increasing set since it is the complement of an increasing one) leads to Finally, taking p i gives (2.13)-(2.14).

Coupling construction
The (harder) sufficient condition of Theorem 2.4 is obtained by showing (in this subsection) the existence of a Markovian coupling, which appears to be increasing under Conditions (2.13)-(2.14) (see Subsection 3.3). Our method is inspired by [8, Propositions 2.25, 2.39, 2.44], but it is much more intricate since we are dealing with jumps, births and deaths. Let ξ t ∼ S and η t ∼ S such that ξ t ≤ η t . The first step consists in proving that instead of taking all possible sites, it is enough to consider an ordered pair of sites (x, y) and to construct an increasing coupling concerning some of the rates depending on η t (x), η t (y) and ξ t (x), ξ t (y) (remember that we choose to take births on y, deaths on x and jumps from x to y) and a small part of the independent rates (by this, we mean deaths from x with a rate depending only on η t (x) and ξ t (x) and births upon y with a rate depending only on η t (y) and ξ t (y)). We do not have to combine any "dependent reaction" R ·,· ξt(x),ξt(y) or jumps rate Γ · ξt(x),ξt(y) on y with any rate R ·,· ηt(z),ηt(y) or Γ · ηt(z),ηt(y) if z is different from x.
Definition 3.2 For fixed (x, y) ∈ S 2 , for all η ∈ Ω and k ∈ N let otherwise. and q(z, w) = p(x, y) if z = x and w = y, 0 otherwise.
An ordered pair of sites (x, y) is an attractive pair for ( S, S) if there exists an increasing coupling for (ξ t , η t ) t≥0 where ξ t ∼ S( Γ, R, h( P ), q), η t ∼ S(Γ, R, h(P ), q). For notational convenience we call these new systems S (x,y) and S (x,y) .
Notice that S (x,y) = S (y,x) , because we take into account births on the second site, deaths on the first one, and only particles' jumps from the first site to the second one. The same remark holds for S (x,y) . In other words in order to see if a pair is attractive, we reduce ourselves to a system with only part of the rates depending on the pair, and a part of the independent rates depending on p(x, y) (P k ηt(y) p(x, y) and P −k ηt(x) p(x, y)).

Proposition 3.3
The process η t ∼ S is stochastically larger than ξ t ∼ S if all its pairs are attractive pairs for ( S, S).
Proof . If for each pair (x, y) we are able to construct an increasing coupling for ( S (x,y) , S (x,y) ), we define an increasing coupling for ( S, S) by superposition of all these couplings for pairs. Indeed: it is a coupling since by Definition 3.2 the sum of all marginals gives the original rates; it is increasing since each coupling for a pair is increasing.
The rest of this section is devoted to the construction of an increasing coupling for (x, y) and (α, β) ≤ (γ, δ).

Definition 3.5
There is a lower attractiveness problem on β if there exists k such that β + k > δ and Π 0,k α,β + Γ k α,β > 0; β is k-bad and k is a bad value (with respect to β). There is a higher attractiveness problem on γ if there exists l such that γ − l < α and Π −l,0 γ,δ + Γ l γ,δ > 0; γ is l-bad and l is a bad value (with respect to γ). Otherwise β is k-good (resp. γ is l-good). There is an attractiveness problem on (α, β), (γ, δ) if there exists at least one bad value.
In other words we distinguish bad situations, where an addition of particles allows lower states to go over upper ones (or upper ones to go under lower ones) from good ones, where it cannot happen. Notice that Definition 3.5 involves addition of particles upon β and subtraction of particles from γ. If we are interested in attractiveness problems coming from addition upon α and subtraction from δ we refer to (β, α), (δ, γ).
We choose to define a coupling rate that moves both processes only if we are dealing with an attractiveness problem, otherwise we let the two processes evolve independently through uncoupled rates. Conditions (2.13) − (2.14) do not involve configurations without an attractiveness problem, so a different construction for them does not change the result. Since N is finite, we can construct the coupling by a downwards recursion on the number of particles involved in a transition.
Our purpose now is to describe a coupling for ( S (x,y) , S (x,y) ), that we denote by H(x, y) (or simply H), which will be increasing under Conditions (2.13) − (2.14).
First of all we detail the construction on terms involving the larger number N of particles and we prove that under Conditions (2.13)-(2.14) none of these coupling terms breaks the partial order: this is the claim of Proposition 3.18. Remark 3.6 By Hypothesis 2.1, at least one of the terms Π 0,N α,β , Π 0,N γ,δ , Π −N,0 α,β , Π −N,0 γ,δ , Γ N α,β and Γ N γ,δ is not null. We assume all these terms (and the smaller ones) positive. Otherwise the construction works in a similar way with some null terms.

Definition 3.7 Let
If there is a lower attractiveness problem on β, then β + N > δ ( N + > 0) and an addition of N particles upon β breaks the partial order. Such a problem comes both from birth (Π 0,N α,β p) and from jump (Γ N α,β p) rates.
Notice that if l ≥ N + , then β + N ≤ δ + l, and additions of N particles upon β and of l particles upon δ do not break the partial order on y. The construction consists in coupling the terms on the left hand side (which involve N particles and break the partial order) to the ones on the right hand side (in such a way that the final configuration preserves the partial order on y and on x) by following a basic coupling idea. We couple jumps on the lower configuration with jumps on the upper one and births on the lower configuration with births on the upper one. Only if this is not enough to solve the attractiveness problem, we mix births with jumps. If there is a higher attractiveness problem on γ, then γ−N < α ( N − > 0) and a subtraction of N particles from γ breaks the partial order. In this case the problem comes from Π −N,0 γ,δ p and Γ N γ,δ p; we use a symmetric construction starting from Condition (2.14) with K = 1, We denote by H k,k,·,· α,β,γ,δ (resp. H ·,·,l,l α,β,γ,δ ) the coupling terms which involve jumps of k (resp. l) particles from x to y on the lower (resp. upper) configuration; H 0,k,·,· α,β,γ,δ (resp. H ·,·,0,l α,β,γ,δ ) are the coupling terms concerning births of k (resp. l) particles on y on the lower (resp. upper) configuration and H −k,0,·,· α,β,γ,δ (resp. H ·,·,−l,0 α,β,γ,δ ) are the symmetric ones for death rates. For instance H k,k,0,l α,β,γ,δ combines the jump of k particles from x to y on the lower configuration and the birth of l particles on y on the upper one.
The coupling construction takes place in three main steps.
Step 1) Suppose both β and γ are N −bad; if this is not the case one of them is good and the construction works in an easier way. We begin with jump rates. We couple the lower configuration N −jump rate Γ N α,β p with jumps on the upper configuration. We first couple it with Γ N γ,δ p, because α ≤ γ, α − N ≤ Then if the lower attractiveness problem is not solved, that is H N,N,N,N α,β,γ,δ = Γ N γ,δ p, we have a remainder of the lower configuration jump rate that we couple with the upper configuration jump rate with the largest change of particles left, N − 1. We go on by coupling the new remainder of Γ N α,β p, if positive, with Γ l γ,δ p at l th step. The final pairs of values we reach are (α − N, β + N ), (γ − l, δ + l), which always preserve the partial order on x since α − N ≤ γ − l when l ≤ N . The partial order on y is preserved only if β + N ≤ δ + l, that is if l ≥ N + . For this reason we stop the coupling between jumps at step N + : this is the meaning of formula (3.9). More precisely, when Γ N α,β > Γ N γ,δ and N − 1 ≥ N + , we get the second coupling rate where J N,l α,β is the remainder of the jump (hence notation J) rate Γ N α,β p left over after the l th step of the coupling construction; (3.12) means that we couple the remainder from the (l + 1) th step, J N,l+1 α,β , with Γ l γ,δ p at the l th step.
We proceed this way until either we have no remainder of Γ N α,β p, or we have reached N + with the remainder Note that, since β + N ≤ δ + l for l ≥ N + , none of the coupled transitions until this point have broken the partial ordering. Proceeding until N + − 1 would break the partial ordering. Therefore we stop at N + the construction in Step 1 and we will couple the remainder J N, N + α,β with upper birth rates at Step 3a.
Step 1 is detailed in Tables 1 and 2, where N d+ corresponds to the first l (going downwards from N ) such that the minimum in (3.12) is J N,l+1 α,β . Table 1: We have to distinguish between two situations: • if the minimum given by H N,N,l,l α,β,γ,δ is always the second term, we do not reach N d+ at N + yet. Since we have decided to stop the coupling construction at step N + , we need to couple with birth rates (Step 3a) in order to solve the attractiveness problem (Table 2 when J N, N + α,β > 0) and we put N d+ = N + − 1.
If there is a higher attractiveness problem, we repeat the same construction for the coupling terms involving the jump rate Γ N γ,δ p starting from (3.11) and we define a value N d− analogous to the previous N d+ . The recursive formula symmetric to (3.12) where J k,N γ,δ represents the remainder of the jump rate Γ N γ,δ p left over after the k th step of the coupling construction. We need to couple the remainder of Γ N γ,δ p with the lower death rates in Step 3a if  If there is no lower (higher) attractiveness problem we put N d+ = N + 1 (N d− = N + 1). If either β or γ is N −good, we use only one of the previous constructions. Suppose for instance that γ is N −good: then the construction involving Γ N α,β p works in the same way, but the symmetric one is not required and we use the coupling terms H ·,·,N,N α,β,γ,δ only to solve the lower attractiveness problems induced by either Γ N α,β p (at Step 1) or by Π 0,N α,β p (at Step 3b). Therefore H 0,N,N,N α,β,γ,δ (defined at Step 3a) might be non null, but we define H −l,0,N,N α,β,γ,δ := 0 for each l > 0. (3.19) If β is N −good a symmetric remark holds.
Step 2) Suppose β is N −bad. The birth rate Π 0,N α,β p could break the partial order on β. We work as in Step 1 and we begin with the coupling term H 0,N,0,N α,β,γ,δ :=(Π 0,N α,β ∧ Π 0,N γ,δ )p. (3.20) If the attractiveness problem is not solved we couple the remainder of Π 0,N α,β p with the birth rate of the upper configuration with the largest change of particles, Π 0,N −1 γ,δ p and going down we couple it with Π 0,l γ,δ p at l th -step, until l = N + . We define recursively starting from (3.20) the terms where B N,l α,β is the remainder of the birth (hence notation B) rate Π 0,N α,β p left over after the l th step of the coupling construction. We proceed as in Step 1 for transitions involving births: while the minimum in (3.21) is the second term we go on downwards in l. As soon as the minimum is the first term, we have no remainder of Π 0,N α,β p, so the lower attractiveness problem is solved and we define N B to be this first such l. Otherwise we have reached N + with the remainder We have to distinguish between two situations: • if the minimum given by H 0,N,0,l α,β,γ,δ is always the second term, we do not reach N B at N + yet. Since we stop the coupling construction at step N + , we need to couple the remainder with jump rates (Step 3b) in order to solve the attractiveness problem (Table 4  • if N B ≥ N + then B N, N + α,β = 0 (see Table 3). Therefore there is no need to continue a coupling involving Π 0,N α,β p, since the attractiveness problem is solved. In this case H 0,N,0,l α,β,γ,δ = 0 for N + ≤ l < N B by definition, we do not need Step 3b and we define H 0,N,l,l α,β,γ,δ := 0 for each l > 0, (3.23) In both cases, we define H 0,N,0,l α,β,γ,δ := 0 for 0 < l < N + . (3.24) The attractiveness problem coming from either birth or death rates is solved. Indeed Lemma 3.9 Under Condition (2.13), Proof . If this is not the case, by (3.17), Table 1 and definition of N d+  Table 3 and definition of N B The sum of (3.25) and (3.26) contradicts (3.9).
Step 3a) Suppose β is N −bad and N d+ = N + − 1. We come back to Step 1, where even if the remaining part of Γ N α,β p was still positive at step N + , we decided to stop. We refer to Table 5. By Lemma 3.9, N B ≥ N + . We use the upper configuration birth rate remaining from Step 2 in order to solve the attractiveness problem: we couple the remainder from Step 1 of Γ N α,β p with the remainder from Step 2 of Π 0,N γ,δ p, where J N,l α,β is the remainder of the jump (J ) rate Γ N α,β p left over after l th step. Notice that J N,N +1 We proceed with (3.28) until the minimum is the first term, in which case there is no remainder of Γ N α,β p, so the attractiveness problem is solved, and we define N dB to be the first such l. In other words when l > N dB , the coupling term is the remainder of the upper configuration l-birth rate, when l = N dB it is the remainder of the lower configuration N -jump rate, and when l < N dB the coupling terms are null.

Remark 3.11
If l > N dB then the minimum of (3.28) is the second term, which depends on Step 1: -if l > N B , by Table 3 α,β = 0 and H N,N,0,l α,β,γ,δ = Π 0,l γ,δ p. In other words even if the minimum is the second term in H N,N,0,l α,β,γ,δ , it could be null, when solving the lower attractiveness problem left no remainder of Π 0,l γ,δ p. It means that positive coupling terms begin below N B .
is N −good: in these cases H N,N,0,l α,β,γ,δ = 0 for each l > 0.  Step 3b) Suppose β is N −bad and N B = N + − 1. We come back to Step 2, where even if the remainder of Π 0,N α,β p was still positive at step N + , we decided to stop. We refer to Table  6. By Lemma 3.9, N d+ ≥ N + . We cannot couple B N, N + α,β with the upper configuration jump rates remaining from Step 2 with l > γ − α, because the final states we would reach are (α, β + N ), (γ − l, δ + l), which break the partial order. Therefore we put H 0,N,l,l α,β,γ,δ := 0 for each l > γ − α (3.29) Then we couple the remainder B N, N + α,β from Step 2 of Π 0,N α,β p with the upper configuration jump rates remaining from Step 1 that do not break the partial order in α to solve the attractiveness problem, that is Γ l γ,δ p with l ≤ γ − α, if the minimum is the second term we proceed downwards in l with Γ l γ,δ p with l < γ − α H 0,N,l,l α,β,γ,δ : where B N,l α,β is the remainder of the birth (B) rate Π 0,N α,β p left over after step l th in Step 3b. Notice that B N,N +1 α,β = B N, N + α,β . We proceed with (3.31) until the minimum is the first term, in which case there is no remainder of Π 0,N α,β p, so the attractiveness problem is solved, and we define N Bd to be the first such l. In other words when l > N Bd , the coupling term is the remainder of the upper configuration l-jump rate, when l = N Bd it is the remainder of the lower configuration N -birth rate, and when l < N Bd the coupling terms are null.  Table 1 and H 0,N,l,l α,β,γ,δ = 0; (3.32) • if l = N d+ ≥ N + (see Table 1), then J N,l+1 α,β ∧ Γ l γ,δ p = J N,l+1 α,β and H 0,N,l,l α,β,γ,δ = Γ l γ,δ p − J N,l+1 α,β = 0 by Table 1; hence In other words if there is no remainder of Γ l γ,δ p after Step 2, the coupling term H 0,N,l,l α,β,γ,δ is null even if the attractiveness problem is not solved yet. It means that positive coupling terms begin below N d+ .
We give it the same value if β is N −good: in these cases H 0,N,l,l α,β,γ,δ = 0 for each l > 0. We refer to Table 6.
We couple an upper configuration jump rate Γ l γ,δ p with a lower configuration birth rate Π 0,k α,β p if β + k > l and α ≤ γ − l in order to solve a lower attractiveness problem, and with a lower configuration death rate Π −k,0 α,β p if α > γ − l in order to solve a higher attractiveness problem: we cannot couple the same higher jump rate both with lower birth and death rates. A symmetric remark holds for lower jump rates.
We explicitly constructed the coupling rates such that the largest change of N particles breaks the partial order: the following proposition proves that coupling H is increasing for such rates, under Conditions (2.13) and (2.14) if β or/and γ are bad values. Indeed it states that all coupling terms that would break the order of configurations are equal to 0. Proof i) Suppose l < N + . Then H 0,N,0,l α,β,γ,δ = H N,N,l,l α,β,γ,δ = 0 by (3.15) and (3.24). • Suppose N d+ = N + − 1, then N B ≥ N + by Lemma 3.9 and H 0,N,l,l α,β,γ,δ = 0 for each l by (3.23). Since H N,N,0,l α,β,γ,δ is null for each l < N dB by Table 5, we prove that N dB ≥ N + . Assume by contradiction that N dB < N + . By Definition of N dB , (3.28), Table 3 and Remark 3.11, (3.39) Hence by (3.43), Claim (ii) is proved by symmetric arguments.

Remark 3.19
As a consequence, Tables 5 and 6 (and the symmetric ones) do not contain any coupling term breaking the partial order between configurations.

Sufficient condition
We complete the proof of Theorem 2.4 by Proposition 3.20 Under Conditions (2.13)-(2.14), H is increasing.
In order to prove Proposition 3.20, we define a new system S (in fact a new pair of systems S := ( S, S) by Remark 3.4), depending on S and H, whose rates are those of S to whom we subtract the coupled rates of H involving changes of N particles.
Our plan consists in working by induction on the largest change of particles n(S) which causes either a lower or a higher attractiveness problem: given the particle system S, it is defined by (3.45)
We prove that if S satisfies Conditions (2.13)-(2.14) and n(S) = N , then H = H(S) is increasing. The induction hypothesis is: if a particle system S * satisfies Conditions (2.13)-(2.14) and n(S * ) ≤ N − 1, then H(S * ) is increasing. We give an outline of the proof: suppose that the induction hypothesis is satisfied. We defined a new system S. By Proposition 3.23, n ≤ N − 1 and by Proposition 3.24, it satisfies Conditions (2.13)-(2.14). Therefore we can use the induction hypothesis and H = H(S) is increasing. This implies, by Proposition 3.27, that H(S) is increasing. Proof . If β is N -bad, we prove that Γ Proof . We prove that for all K, h, j, m, I a , I b , I c , I d in Theorem 2.4, (3.47) We prove (3.46). Since, by symmetry, the proof of (3.47) is similar, we skip it.
Remark 3.25 Let A = {a ∈ X : a ≤ K, j a ≥ N + }, then for each k > δ − β + j i > N such that j i ≥ j a we have Γ k α,β = 0 by Definitions of N and 3.21. Therefore implies (3.46) and we can suppose without loss of generality If γ is N -bad a similar remark involving N − and variables N d− , N D , N Dd and N dD holds by symmetry.
If β is N -good then Condition (3.46) is trivially satisfied. We suppose that both β and γ are N -bad. If γ is N -good the proof is similar but easier, then we skip it. Hence we suppose N + > 0 and N − > 0. By Proposition 3.23, Π 0,N α,β = 0 and by Definition 3.21, Π 0,l α,β = Π 0,l α,β for each l < N . Therefore By Remark 3.9 either N B ≥ N + or N d+ ≥ N + . We detail the case N d+ ≥ N + , which contains all the technical difficulties of the proof. The other proof (for N d+ = N + − 1) is similar.
• If N / ∈ I a , (3.61) writes   Case B: Suppose N = N d+ . In this case H N,N,l,l α,β,γ,δ = 0 for each l < N , the left hand side of (3.57) is equal to l∈I b \{N } Γ l γ,δ p, Formula (3.51) still holds but terms H k,k,N,N α,β,γ,δ are not null. One works as in Case A, when N d− = N , by using the symmetric construction of H k,k,N,N α,β,γ,δ and by checking the condition in different cases.
The equality of recursive formulas corresponding to Step 3b is proved in a similar way.
Since all terms are defined by the same downwards induction formula, we just need to check that they coincide for the initial coupling rates. The one involving Γ N −1 α,β p and Γ N −1 γ,δ p is given by (3.64) with l = N − 1, that is for death rates by symmetric arguments. Therefore H and H are identical for all coupling rates involving a change of less than N particles which cause an attractiveness problem: the conclusion includes uncoupled rates. The claim follows since H is increasing and by Proposition 3.18.

Proof of Proposition 3.20
We do an induction in two steps. First of all we suppose that β is N -bad and γ is N -good, that is β + N > δ and γ − N ≥ α. We skip the similar symmetric case. Let S * N := S and suppose that Proposition 3.20 holds for each system S * such that n(S * ) ≤ N − 1. Notice that since γ is N -good, then γ is n(S * )-good for each S * and, by Remark 3.22, definition of n(S * ) involves for each S * only rates that cause a lower attractiveness problem. Remember that Definition 3.21 of S depends on the original system, that is S = S(S). We define S * j = S(S * j+1 ) and n j = n(S * j ) for j ∈ N, j ≤ N − 1. If β is (N − 1)-bad, by Proposition 3.23, n N −1 ≤ N − 1, and S * N −1 satisfies Conditions (2.13)-(2.14) by Proposition 3.18. We define a coupling H N −1 for S * N −1 as in Definition 3.15, and by induction hypothesis it is increasing. By Proposition 3.27, if H N −1 is increasing then so is H N . We have to check the induction basis: we proceed downwards with definitions of the new systems S * j until β + j ≤ δ. If β < δ, then S * i with j = δ − β is attractive, since β + j ≤ δ and there are no attractiveness problems. If β = δ, then the attractiveness of a system with n 1 = 1 under Condition (2.13) is proved in Appendix A.
We proved that all systems such that there is only a lower (or only a higher) attractiveness problem are attractive under Conditions (2.13)-(2.14). Now we prove that this holds also for a system where both β and γ are N -bad.
If β + N > δ and γ − N < α, the definition of n(S) is given by (3.45) but everything works in a similar way: we define S * N = S, S * j = S(S * j+1 ) and n j = n(S * j ) for j ≤ N − 1. If both β and γ are (N − 1)-bad, by Proposition 3.23, n N −1 ≤ N − 1 and S * N −1 satisfies Conditions (2.13)-(2.14) by Proposition 3.18. We define a coupling H N −1 for S * N −1 as in Definition 3.15; by induction hypothesis it is increasing. By Proposition 3.27, if H N −1 is increasing then so is H N . We have to check the induction basis: we proceed downwards with the definition of a new system S * j until either β or γ are j-good, that is until j = (δ − β) ∨ (γ − α). If δ − β = γ − α, then only one attractiveness problem is present and the claim follows from the first part of the proof. If j = δ − β = γ − α > 0, then S * j is attractive since there are no attractiveness problems; if j = δ − β = γ − α = 0, then the attractiveness of a system with n 1 = 1, β = δ and γ = α under Conditions (2.13)-(2.14) is proved in Appendix A.

Sufficient conditions on more general systems
In order to show the sufficient conditions of Theorem 2.4, we restricted ourselves to transition rates on a given pair of sites (x, y): if Conditions (2.13)-(2.14) are satisfied we can construct an increasing coupling for the system S (x,y) and the final increasing coupling is given by superposition of couplings for all pairs of sites (see Section 3.2 and Proposition 3.3). We use neither the translation invariance of P ±k η(x) , nor the fact that the smaller and the larger systems share the same p(x, y). Therefore we can state the result for a more general pair of systems G = {R 0,±k η(x),η(y) (x, y), Γ k η(x),η(y) (x, y), P ±k η(x) (x), p(x, y))} and G = { R 0,±k η(x),η(y) , Γ k η(x),η(y) , P ±k η(x) (x), p(x, y)}: Corollary 3.28 A particle system η t ∼ G is stochastically larger than ξ t ∼ G if for each (x, y) ∈ S 2 , (α, β) ≤ (γ, δ), (α, β) ∈ X 2 , (γ, δ) ∈ X 2 , with p = p(x, y) and p = p(x, y) Such conditions are not necessary. We use Corollary 3.28 with the comparison technique with oriented percolation in [3] to prove survival of species in metapopulation models.

Proof of Theorem 2.10
Since F (x, y) = F (0, y − x) = F (y − x, 0), with a slight abuse of notation we write F (y − x) instead of F (x, y).
Proof of Proposition 2.12.
We treat separately the cases γ = 0 and γ > 0. i) If γ = 0 then the Dirac measure δ 0 is invariant. We denote by P(·) the independent coupling measure and by E(·) its expected value. We fix x ∈ S and we compute the generator on F (η M t (x)) in Definition ( Notice that if k = 0 then η t (x) = k + l = l. If k ≥ 0 and k + l = M , then the last term in the right hand side of (4.3) is smaller or equal to β y∼x ζ t (y)u l + lφu l − lu l−1 .
The same inequality holds when l = 0 and k = M , that is ξ t (x) = η t (x) = M , since the last term in the right hand side of (4.3) is null and lφu l − lu l−1 = 0. If k = 0 and l = 0, Now we prove that Definition 2.14 satisfies the hypothesis of the u-criterion under (2.29). We begin with a technical proposition.
Proposition 4.2 Let (u l (ǫ)) l∈X be given by Definition 2.14. If φ < 1 and then there existsǭ > 0 such that u l (ǫ) is positive, decreasing in l for each l ∈ X and in ǫ for each 0 < ǫ ≤ǭ.