S Nundloll
email: sapna.nundloll@inria.fr

L Mailleret
email: ludovic.mailleret@sophia.inra.fr

F Grognard

Global Stability in a General Impulsive Biological Control Model with Harvest

In this paper, the effects of periodic partial harvesting of a continuously grown crop on augmentative biological control are analysed. Partial harvesting can remove a proportion of both pests and biological control agents, so its influence on the control efficiency cannot be a priori neglected. An impulsive model consisting of a general predator-prey model in ode is used. It is augmented by a discrete component to depict releases of biological control agents and the periodic partial harvesting. A stability condition for pest eradication is expressed as the minimal value of the budget per unit time to spend on predators. We consider the partial harvesting period to be fixed so that the only manipulated variable is the release period. One period is taken as the integer multiple of the other. We show that when the releases are carried out more often than the harvests, the release period influences the minimal budget. Conversely, there is no effect on this budget when releases take place as often as or less frequently than harvests.

INTRODUCTION

Biological control is the reduction of pest populations to harmless levels through the release of their natural enemies. The target pest species and the setting usually determine the type of control required, namely whether pest eradication is necessary or not. For an exhaustive list of definitions and applications, we refer the reader to [START_REF] Eilenberg | Suggestions for unifying the terminology in biological control[END_REF]. In this paper, we consider the implementation of biological control through natural predators to protect continuously grown crops which have zero tolerance to pest invasions. Examples are flowers of high cosmetic value, as well as fruits and vegetables, which, upon being attacked by a pest, become inedible. There are three aspects in the formulation of this problem:

• Firstly, the most appropriate method of protection in such a case is inundative or augmentative control. This involves a calculated number of predators to be periodically injected into the ecosystem, independently of the detection of pest insects. Such a preemptive approach is becoming increasingly popular as it appears to achieve more acceptable pest control (see [START_REF] De | Biological control of thrips on ornamental crops: interactions between the predatory mite Neoseiulus Cucumeris (acari: phytoseiidae) and western flower thrips Frankliniella Occidentalis (thysanoptera: thripidae) on cyclamen[END_REF], [START_REF] Jacobson | Suppressing establishment of frankliniella occidentalis pergande (thysanoptera : Thripidae) in cucumber crops by prophylactic release of amblyseius cucumeris oudemans (acarina : Phytoseiidae)[END_REF] for real life experiments). • Secondly, pest eradication is desirable.

• Thirdly, over their growing period, these crops are partially harvested on a regular basis. This practice is likely to remove a proportion of the insects -including the predatory ones -present throughout the plantation. Harvest has already been shown in [START_REF] Volterra | Variations and fluctuations in the numbers of coexisting animal species, volume 22 of The Golden Age of Theoretical Ecology[END_REF] among others to influence, even counterintuitively, predator-prey dynamics. Hence, it is incorporated into the model under study.

⋆ This work was partially funded by the INRA ECOGER program.

The model consists of ODEs augmented by a discrete component to integrate the effect of partial harvest and releases that by their very nature are discrete phenomena.

Few papers in the literature on impulsive crop protection seem to focus on stability of the pest-free state: [START_REF] Liu | The dynamical behaviors of a Lotka-Volterra predator-prey model concerning integrated pest management[END_REF], [START_REF] Negi | Dynamics in a Beddington-DeAngelis prey-predator system with impulsive harvesting[END_REF] for instance study the chaotic properties of similar models. Yet this is of practical importance, especially for high valued crop cultures.

We are able to express a stability condition as the minimal number of predators per unit time required to drive the pests to zero. This time unit can be chosen to be the budget year over which the number of predators bought by the grower corresponds to the budget allocated to biological control. In line with the work of [START_REF] Mailleret | Positive Systems[END_REF], we investigate how the frequency of releases is to be varied with respect to the (fixed) harvesting frequency to minimise the minimal budget value. We consider the harvest period as the reference since it is set by market constraints and crop dynamics. One period is taken as the integer multiple of the other. This feature is key in solving for the stability condition.

We show that for a given harvest period, when releases take place less often or as often as harvests, the minimal budget is at a calculated value which is independent of release period. However, when releases take place more often than harvests, the minimal budget required always exceeds this value.

In section 2, the system model is described. Next, in section 3, the mathematical analysis of the system's stability and the formulation of the stability condition in terms of the minimal budget are presented. A brief interpretation of the mathematical results follows in section 4. Finally, we conclude with a discussion on their implications.

MODEL DESCRIPTION

We consider the case at the onset of pest invasion where the crop -the pest food supply -is in abundance. At this stage, it is sufficient to model a two-dimensional system consisting only of the pest x and predator y species.

         ẋ = f (x) -g(x)y ẏ = h(x)y -dy x(nT + h ) = (1 -α x )x(nT h ) y(nT + h ) = (1 -α y )y(nT h ) + δ (nT h mod T r ) µT r y(mT + r ) = (1 -δ (mT r mod T h ) α y )y(mT r ) + µT r
(1) ∀n , m ∈ N and where 0 ≤ α x , α y ≤ 1.

Continuous part f (x) is the growth velocity or feeding input of the pests. It represents the growth function of the pest species and in our model, it also encompasses any non-predatory losses of the pest population (e.g. logistic growth). We assume that the predator population is never large enough for intra-predator interaction to take place so the functional and numerical responses can be expressed solely in terms of the prey numbers, i.e. as g(x) and h(x) respectively. The functions discussed are not specified so they are representative of most functions reported in the literature. Only the following hypotheses are made. Hypothesis 1. Let f (x), g(x) and h(x) be locally Lipschitz continuous in R + such that

• f (0) = 0 • g(0) = 0, g ′ (0) > 0 and g(x) > 0 ∀x > 0 • h(0) = 0 and h(x) > 0 ∀x > 0 • f (x)
g(x) and g(x)

x are upper bounded for x ≥ 0 g(x) is increasing for small pest population levels. We also consider that, in the presence of pests, predation always takes place with a negative impact on x (g(x) > 0) and a positive impact on y (h(x) > 0). In classical density dependent models, g(x) is bounded or linear, so that g(x)

x is always bounded. The boundedness of f (x)

g(x)
means that there is no value of x where the pest growth f (x) overwhelmingly dominates the predation g(x), which would render the biological control impossible.

Discrete part Partial crop harvests and predator releases occur respectively every T h and T r . α x and α y represent the respective proportions of the prey and predator populations affected at each harvest. These parameters are allowed be different since in reality, it is likely that each species occupies different parts of the plant. We also assume that the insects are uniformly distributed throughout our plantation. The + -superscript denotes the instant right after the relevant impulse is applied. The δ-function is defined thus to identify instants of simultaneous partial harvest and predator release.

δ(θ) = 1 if θ = 0 0 otherwise
µ refers to the total number of predators purchased per time unit. So as to give an economic dimension to the solution, the time unit will be the budget time period for which a crop grower will invest that number µ. In this same unit, µT r is the number of predators released every T r .

MATHEMATICAL ANALYSIS

While the formalism in model ( 1) is more general, we restrict ourselves to the case where either one of the periods (release or partial harvests) is the integer multiple of the other. We study the stability of the system around x = 0. In addition of being invariant, it is the target state. Our analysis takes place separately for the case when releases are more frequent than harvests, and when they are less frequent.

Pest-free stability analysis

Releases more frequent than harvests Proposition 2. Let T h = kT r where k ∈ N * and Hypotheses 1 be satisfied. Then, in the absence of pests, model (1) possesses a globally stable periodic solution (x ph (t) , y ph (t)) such that

x ph (t) = 0 y ph (t) = y * e -d(t mod T h ) + µT r e -d(t mod Tr) Λ (2) with Λ = ⌊ t mod T h Tr ⌋-1 j=0
e -jdTr , and where

y * = 1-e -dT h 1-e -dTr (1 -α y ) + α y 1 -(1 -α y )e -dT h µT r (3)
Proof. When T h = kT r , in the absence of pests and using Hypotheses 1, the system is simplified to

         ẋ = 0 ẏ = -dy x(mT + r ) = (1 -δ (m mod k) α x )x(mT r ) y(mT + r ) = (1 -δ (m mod k) α y )y(mT r ) + µT r ∀m ∈ N
The pest population stays nil since in the absence of pests, their population does not change either. The solution

x ph (t) = 0 is trivial. On the other hand, the predator population will vary according to the number of predators manually injected into the system and, since the population is nonzero, according to the partial harvest effect. The absence of their source of food will cause an exponential decay of the population. We demonstrate that these forces will provoke the predator population to reach a periodic pattern of period equal to T h , which we shall refer to as the period of reference. The instant following a coinciding partial harvest and release is taken as the point of reference.

To prove Proposition 2, we first express the predator population right after a release in terms of the point of reference as

y(nT h + iT + r ) = y(nT + h )e -idTr + µT r i-1 j=0 e -jdTr (4)
where i ∈ [0, 1, . . . , (k -1)]. This can be verified by induction [START_REF] Nundloll | The effects of partial crop harvest on biological pest control[END_REF].

To evaluate the evolution of y according to the period of reference T h , we need to calculate the value of y((n + 1)T + h ), which is equivalent to y(nT h + kT + r ), in terms of y(nT + h ). At this point however, we suppose that partial harvesting takes place before predator release so as not to directly waste predators; so we first express it in terms of y(nT h + (k -1)T + r ) then expand the expression using (4) as follows

y (n + 1)T + h = y nT h + (k -1)T + r e -dTr (1 -α y ) + µT r =   y(nT + h )e -d(k-1)Tr + µT r k-2 j=0 e -jdTr   e -dTr (1 -α y ) + µT r = y(nT + h )e -dT h + µT r (1 -α y ) k-1 j=1 e -jdTr + µT r = y(nT + h )e -dT h + µT r   (1 -α y ) k-1 j=0 e -jdTr + α y   (5)
In this linear dynamical system, the coefficient of y(nT + h ), e -dT h is less than one in magnitude, so the sequence will converge to a limit, the equilibrium of ( 5). This equilibrium yields (3) and the convergence of y(t) to a periodic solution y ph (t) based on y * . Now that we have established the existence of the periodic solution y ph (t), we seek to formulate it. We focus on a reference period over nT h < t ≤ (n + 1)T h during which y ph (t) is piecewise continuous, with the continuous components separated by predator releases. The continuous intervals are defined over nT h + iT r < t ≤ nT h + (i + 1)T r where i ∈ [0, 1, . . . , k -1]. For a given value of t, the value of i is easily identified as being i = ⌊ t mod T h Tr ⌋. The value of y ph (t) is then of the form y ph (t) = y ph (nT h + iT + r )e -d(t mod Tr) and, from (4) with y(nT + h ) = y * , we have that

y ph (nT h + iT + r ) = y * e -idTr + µT r i-1 j=0 e -jdTr
so that

y ph (t) =   y * e -idTr + µT r i-1 j=0 e -jdTr   e -d(t mod Tr) = y * e -d(t mod T h ) + µT r e -d(t mod Tr) i-1 j=0 e -jdTr
Using the original notation where i = ⌊ t mod T h Tr ⌋ leads to the same form as proposed in (2), thereby completing our proof.

Releases less frequent than harvests When harvesting is more frequent than the release of predators, we have a similar result about the existence of a periodic solution. Proposition 3. Let T r = kT h where k ∈ N * and Hypotheses 1 be satisfied. Then, in the absence of pests, model (1) possesses a globally stable periodic solution (x pr (t) , y pr (t)) = 0, y * e -d(t mod Tr) (1α y )

⌊ t mod Tr T h ⌋ where y * = µT r 1 -(1 -α y ) k e -dTr
The proof is constructed in the same manner as for Proposition 2, and is detailed in [START_REF] Nundloll | The effects of partial crop harvest on biological pest control[END_REF].

The form of the y ph and y pr functions is illustrated on Figure 1. Harvests are apparent at every nT h instant, while the release of predators dominates the harvest at every mT r instant. Between impulses, the population decays exponentially since it has no prey to feed on.

Stability analysis in the presence of pests

We will now show that, when preys are present at the initial time, convergence of the predator population also takes place to that same periodic solution, while the preys go extinct provided some condition on the parameters is verified.

Since we will study the convergence of the solutions to (0, y p (t)) (where the p subscript stands as well for ph or pr), it will be convenient to describe the system in terms of the deviation coordinates with respect to the reference periodic solution, (x(t), ỹ(t)) = (x(t)x p (t), y(t)y p (t)). This yields ẋ = f (x)g(x)y = f (x)g(x)(ỹ + y p (t)) ẏ = h(x)ydyh(x p )y p + dy p = h(x)(ỹ + y p (t))dỹ

The harvest effects on x and for ỹ are obviously unchanged compared to those on x and y. On the other hand, the release effects on y disappear in ỹ; indeed, we have

ỹ(mT + r ) = y(mT + r ) -y p (mT + r ) = y(mT r ) + µT r -(y p (mT r ) + µT r ) = ỹ(mT r )
To perform a global and a local stability analysis, we will also need the computation of the linear approximation of the deviation system around the periodic solution (0, y p (t)):

ẋ = (f ′ (0) -g ′ (0)y p (t))x ẏ = h ′ (0)y p (t)x -dỹ (6)
We obtain two different constraints for the Local Asymptotic Stability (LAS) and Global Asymptotic Stability (GAS) of the periodic solution in system (1). The latter is obviously stronger than the former, but is sufficient in the case where pests outbreaks do not immediately take large proportions.

Releases more frequent than harvests In order to state the following theorem, we first need to define the function

µ h (S, r) = d S + ln (1 -α x ) rT h 1 1 - αy(1-e -dT h ) 1-(1-αy)e -dT h Γ where Γ = e -dT h /k
k(1-e -dT h /k ) , and S and r are constants. This function is increasing in S and r because the sign of the partial derivatives is determined by the sign of factor Γ, which is positive [START_REF] Nundloll | The effects of partial crop harvest on biological pest control[END_REF]. Theorem 4. When T h = kT r with k ∈ N * , the solution (x(t), y(t)) = (0, y ph (t)) of ( 1) is LAS iff

µ > µ h f ′ (0) g ′ (0) , g ′ (0) (7)
and is GAS if µ > µ h (S g , r g ) ( 8)

where S g = sup x≥0 f (x) g(x) and r g = sup x≥0 g(x)

x

Proof. We start with the proof of global convergence under condition (8). In this proof, we will first show that x goes to zero, from which we will derive that ỹ goes to 0 also (so that y(t) converges to y ph (t)).

Let the initial condition for system (6) be (x 0 , ỹ0 ) at time t 0 = 0 + , that is after the harvest and the predator release that take place at the initial time. Analysing (6) and noting that y ph (t) + ỹ = y(t) ≥ 0, we have ẏ ≥ -dỹ so that ỹ(t) ≥ min(0, ỹ0 )e -dt .

In order to analyse the ẋ equation, we define the function

G(x) = x x0 1 g(s) ds (9)
which can easily be seen to be an increasing function of x since g(s) > 0.

We also have that

g(s) < sup x≥0 g(x) x s, therefore lim x > →0 G(x) = -∞.
In order to show the extinction of the pests we will then prove that G(x) goes to -∞ as t goes to infinity. We start by writing the G-dynamics:

dG(x) dt = 1 g(x) ẋ = f (x) g(x) -ỹ -y ph (t) ≤ f (x) g(x)
-min(0, ỹ0 )e -dty ph (t)

We then consider the evolution of G between two successive harvests, that is between the times nT + h and (n+1)T h for a given n:

G(x((n + 1)T h )) ≤ G(x(nT + h )) + (n+1)T h nT + h f (x(s)) g(x(s))
-min(0, ỹ0 )e -dsy ph (s) ds

Since no impulse is present inside the integral, we can drop the + superscript in its lower extremity.

We now analyse how the harvest that takes place at time (n + 1)T h impacts G. We have

G(x((n + 1)T + h )) = G(x((n + 1)T h )) + x((n+1)T + h ) x((n+1)T h ) 1 g(s) ds ≤ x((n+1)T + h ) x((n+1)T h ) 1 g(s) ds + G(x(nT + h )) + (n+1)T h nT h f (x(s)) g(x(s)) -min(0, ỹ0 )e -ds -y ph (s) ds (10)
The first term represents the influence of harvest on G and can easily be approximated because x((n + 1)T h ) > x((n + 1)T + h ) = (1α x )x((n + 1)T h ). Substituting in S g and r g , we have 11) into (10) then yields a bound on the application between successive moments after harvest.

(1-αx)x((n+1)T h ) x((n+1)T h ) 1 g(s) ds ≤ (1-αx)x((n+1)T h ) x((n+1)T h ) 1 r g s ds = ln(1 -α x ) r g (11) Introducing (

G(x((n + 1)T

+ h )) ≤ ln(1 -α x ) r g + G(x(nT + h )) + (n+1)T h nT h f (x(s)) g(x(s))
-min(0, ỹ0 )e -dsy ph (s) ds

We can now evaluate an upper-bound for G at any time t ≥ 0. Defining l as the integer part of t T h , we have:

G(x(t)) -G(x 0 ) ≤ t 0 f (x(s)) g(x(s)) -min(0, ỹ0 )e -ds -y ph (s) ds +l ln(1 -α x ) r g ≤ t 0 S g -min(0, ỹ0 )e -ds -y ph (s) ds + l ln(1 -α x ) r g = - t 0 min(0, ỹ0 )e -ds ds + t lT h [S g -y ph (s)] ds +l T h 0 [S g -y ph (s)] ds + l ln(1 -α x ) r g = min(0, ỹ0 ) d (e -dt -1) + t lT h [S g -y ph (s)] ds +l T h 0 [S g -y ph (s)] ds + l ln(1 -α x )
r g The first two terms are bounded (the first one is obvious and the second one is upper-bounded by S g T h ). We then have to analyse the third one, which has been obtained through the periodicity of y ph (t) and the fourth in order to know if G(x(t)) goes to -∞ when t goes to infinity. In fact, it suffices to have

T h 0 [S g -y ph (s)] ds + ln(1 -α x ) r g < 0
to achieve this. It is more clearly rewritten in the form

T h 0 y ph (t)dt > S g T h + ln(1 -α x ) r g (12)
In order to obtain (8), we are now left with the computation of T h 0 y ph (t)dt, which is detailed in [START_REF] Nundloll | The effects of partial crop harvest on biological pest control[END_REF].

T h 0 y ph (t)dt = µT h d 1 - α y 1 -e -dT h 1 -(1 -α y ) e -dT h Γ (13)
Introducing ( 13) into (12) then yields (8), which shows that this last condition is sufficient for having x going to 0 as t goes to ∞.

Since x goes to zero, there exists a finite time t f after which h(x) ≤ d 2 for all times. Therefore, after this time,

ẏ = h(x)(y ph (t) + ỹ) -dỹ ≤ h(x)y ph (t) - d 2 ỹ
We have seen that h(x)y ph (t) goes to zero as t goes to infinity; so does ỹ also.

In order to have the global asymptotic stability, we are only left with the local asymptotic stability to prove. In order to do that, we only have to consider the discrete system that maps the state at time nT + h onto the state at time (n + 1)T + h with respect to the linear equation ( 6) and the discrete part [START_REF] Shulgin | Pulse vaccination strategy in the sir epidemic model[END_REF]]. After some computations, we obtain:

x ỹ (n + 1)T

+ h = B x ỹ nT + h where B =   (1 -αx)e (n+1)T h nT h f ′ (0)-g ′ (0)y ph dτ 0 ‡ (1 -αy )e -d (n+1)T h nT h dτ   Since the matrix is triangular, it is stable if |B 11 | < 1, i.e. (n+1)T h nT h y ph dτ > f ′ (0)T h + ln(1 -α x ) g ′ (0) (14)
Similarly to what was done earlier, it can be shown that ( 14) is equivalent to (7), so that the necessary and sufficient condition for local stability is proven.

It is directly seen that ( 14) is satisfied when ( 12) is because µ h (S, r) is increasing in S and r and we have

f ′ (0) g ′ (0) = lim x ≥ →0 f (x) g(x) ≤ sup x≥0 f (x) g(x) g ′ (0) = lim x ≥ →0 g(x) x ≤ sup x≥0 g(x) x
This completes the proof of global stability, since we have shown global convergence and local stability when (8) is satisfied.

Releases less frequent than harvests If we now consider the case where predators releases take place less often than harvests, we also obtain global and local stability results based on the following function

µ r (S, r) = d S + ln(1 -α x ) rT h 1 -(1 -α y )e -dT h
1e -dT h which is increasing in S and r since the last fraction is positive and α x ≤ 1. Theorem 5. When T r = kT h with k ∈ N * , the solution (x(t), y(t)) = (0, y pr (t)) of ( 1) is LAS iff

µ > µ r f ′ (0) g ′ (0) , g ′ (0) (15)
and is GAS if

µ > µ r sup x≥0 f (x) g(x) , sup x≥0 g(x)
x

This proof does not depart very much from the one of Theorem 4, and is given in detail in [START_REF] Nundloll | The effects of partial crop harvest on biological pest control[END_REF].

Comment As we have seen, when the condition (8) or ( 5) is satisfied, the extinction of the pests is GAS. When the local condition ( 7) or ( 15) is not verified, the extinction of the pests is not stable and a bifurcation analysis similar to what is done in [START_REF] Lakmeche | Bifurcation of non-trivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment[END_REF], [START_REF] Liu | The dynamical behaviors of a Lotka-Volterra predator-prey model concerning integrated pest management[END_REF] would show the presence of a limit cycle when µ is close to the limit. When µ satisfies condition ( 7) or ( 15) only, the pests extinction is locally stable and we cannot rule out that it is globally stable (since our global condition is only sufficient). Such a budget has the advantage of being smaller than the one that guarantees global stability. It allows for good control of limited pest invasions; however the culture is at risk of being destroyed by a large pest outbreak. We also need to note that when S + ln(1-αx) rT h < 0, for any of the local or global condition, the condition is trivially verified. Indeed, it implies simply that no biological control is needed for exterminating the pests; in fact, the partial harvesting is effective enough for this purpose (as α x is large enough).

INTERPRETATION OF RESULTS

It is clear from (15) that µ r is independent of T r . The latter's influence for non-trivial values of µ h is formulated in the following theorem, Theorem 6. Let T h = kT r where k ∈ N * .

The minimal budget is monotonically decreasing with respect to the release period T r for non-negative values of µ h , i.e.

∂µ h ∂T r < 0 (16)
the proof of which can be found in [START_REF] Nundloll | The effects of partial crop harvest on biological pest control[END_REF].

We can deduce that we hit the smallest minimal value for the budget for the largest possible T r in this case that corresponds to when k = 1. This happens when the release frequency equals the partial harvest frequency. Figure 2 represents the analytical results formulated in the Theorem 6, and obtained from the expressions ( 7) and (15). , and the rate of growth f ′ (0), functional response g ′ (0) and numerical response h ′ (0) when the ecosystem is pestfree, i.e. x p (t) = 0, are all equal to 1.

These results imply that it is clearly less costly to protect a greenhouse culture for lower frequencies of release. Of additional economic interest, in this case, the biological treatment is always combined with partial harvesting, so that there is little or no extra cost linked to the presence of workers on-site. However, this result runs counter with that obtained by [START_REF] Mailleret | Positive Systems[END_REF] who demonstrated that the higher the release frequency, the smaller the worst-case damages. Merging the two, for the set of possibilities considered, seems to indicate that the harvest frequency is a threshold that should not be exceeded when releasing predators for efficient biological control.

CONCLUSION

Our study was based on the biological control of continuously grown crops which are partially harvested on a regular basis. We formulated the stability condition for the system in terms of the minimal number of predators to inject over a period of time, formalising mathematically the concept of a minimal budget.

We investigated the combined effects of releases and partial harvests on this minimal budget in terms of the relative frequencies of their implementation. We considered the case where these two events occurred at periods such that one was the integer multiple of the other, and with the two events coinciding over the longer period. We found that the harvest frequency provided a threshold for the release frequency, below which biological control is less costly.

Combined with the findings of [START_REF] Mailleret | Positive Systems[END_REF], we concluded that for the set of possibilities that was studied, the current 'best' strategy is when release and harvest frequencies are equal.

This approach has, however, its shortcomings. One is that it is not yet generalised to other scenarios where neither period is the integer multiple of the other. It is possible that these intermediate ratios induce other dynamics in the system. Whether they might stabilise it giving even lower minimal budget values or favour chaos remains to be seen.

Nevertheless, we consider that our results already have a practical economic advantage. Indeed, coinciding periods imply little or no additional costs incurred in terms of labour: the task of predator release can be assigned to workers in charge of partial harvesting.

Fig. 1 .

 1 Fig. 1. Top Periodic solution y ph (t) in the case where k = 3. Releases of predators are apparent at every mT r instant. Bottom Periodic solution y pr (t) in the case where k = 3.Harvests are apparent at every nT h instant, while the release of predators dominates the harvest at every mT r instant. Between impulses, the population decays exponentially since it has no prey to feed on.

Fig. 2 .

 2 Fig. 2. Variation of the minimal number of predatorsrequired per budget year µ as a function of release to harvest period ratio. Parameters are given the values (in arbitrary units): α x = α y = 0.5, d = 1, and the rate of growth f ′ (0), functional response g ′ (0) and numerical response h ′ (0) when the ecosystem is pestfree, i.e. x p (t) = 0, are all equal to 1.


E Tassi
email: tassi@cpt.univ-mrs.fr

Hamiltonian derivation of a gyrofluid model for collisionless magnetic reconnection

We consider a simple electromagnetic gyrokinetic model for collisionless plasmas and show that it possesses a Hamiltonian structure. Subsequently, from this model we derive a two-moment gyrofluid model by means of a procedure which guarantees that the resulting gyrofluid model is also Hamiltonian. The first step in the derivation consists of imposing a generic fluid closure in the Poisson bracket of the gyrokinetic model, after expressing such bracket in terms of the gyrofluid moments. The constraint of the Jacobi identity, which every Poisson bracket has to satisfy, selects then what closures can lead to a Hamiltonian gyrofluid system. For the case at hand, it turns out that the only closures (not involving integro/differential operators or an explicit dependence on the spatial coordinates) that lead to a valid Poisson bracket are those for which the second order parallel moment, independently for each species, is proportional to the zero order moment. In particular, if one chooses an isothermal closure based on the equilibrium temperatures and derives accordingly the Hamiltonian of the system from the Hamiltonian of the parent gyrokinetic model, one recovers a known Hamiltonian gyrofluid model for collisionless reconnection. The proposed procedure, in addition to yield a gyrofluid model which automatically conserves the total energy, provides also, through the resulting Poisson bracket, a way to derive further conservation laws of the gyrofluid model, associated with the so called Casimir invariants. We show that a relation exists between Casimir invariants of the gyrofluid model and those of the gyrokinetic parent model. The application of such Hamiltonian derivation procedure to this two-moment gyrofluid model is a first step toward its application to more realistic, higher-order fluid or gyrofluid models for tokamaks. It also extends to the electromagnetic gyrokinetic case, recent applications of the same procedure to Vlasov and driftkinetic systems.

Introduction

Fluid models represent a widespread and effective tool for investigating important phenomena in fusion plasmas such as instabilities, turbulence and reconnection events. Indeed, fluid models offer a considerable advantage, in terms of required computational resources, with respect to kinetic models. On the other hand, compared to these, they obviously suffer from limitations in the range of scales and frequencies of the phenomena that they can describe. The derivation of sophisticated fluid models aiming at partially remedying such limitations is an active line of research since a long time. An essential part of the problem is related to the closure adopted to truncate the fluid hierarchy of equations obtained by taking moments of a parent kinetic model. A considerable effort, directed also toward applications to space plasma turbulence, has been carried out to derive fluid and gyrofluid models (see, e.g. Refs. [1][2][3][4][5][6][7][8]) by taking moments of kinetic or gyrokinetic models and imposing closures designed to satisfy desired properties, such as for instance consistency with kinetic linear theory or energy conservation. Relatively little attention has been paid, however, to ensure that the adopted closures satisfy a further criterion, which is the preservation of a Hamiltonian structure in the derivation of the non-dissipative part of the model. In other words, the parent kinetic model, in its non-dissipative limit, is supposed to possess a Hamiltonian structure, and, unless forcing or dissipative terms are voluntarily added when applying the fluid closure, the resulting fluid model is also supposed to possess a Hamiltonian structure. Apart from being an important feature from a fundamental point of view, the preservation of a Hamiltonian structure provides also some practical advantages. On one hand, it provides an unambiguously defined conserved total energy (the Hamiltonian functional). Moreover, given that fluid models for plasmas are often expressed in Eulerian form, their Hamiltonian structure is typically noncanonical (see, e.g. Ref. [9]), which implies the existence of additional conserved quantities, denoted as Casimirs. Their knowledge can provide useful information on the nonlinear dynamics of the fluid model. In the presence of dissipation, the decay rate and the cascade of such invariants can provide a further way to characterize plasma turbulence. Moreover, they can be used as additional invariants to test the conservation properties of numerical codes. Namely because with Hamiltonian plasma models one has to deal with noncanonical Poisson brackets, a delicate point in preserving the Hamiltonian structure throughout the derivation, is that of not violating the Jacobi identity. The latter is one of the properties defining a Poisson bracket. Although it is always satisfied for canonical Poisson brackets, showing its validity for noncanonical Poisson brackets is often far from obvious.

Recently, derivation procedures that preserve the Hamiltonian structure of the parent model have been applied to kinetic and drift-kinetic models [10][11][12]. A key element in these procedures consists of inserting a generic fluid closure relation in the bracket obtained from the parent kinetic model, and then imposing the constraint that the Jacobi identity be satisfied. Such constraint selects fluid closures that automatically lead to Hamiltonian fluid models. The purpose of the present paper is two-fold. In the first place, we show in Sec. 2 the existence of a Hamiltonian structure for a simple gyrokinetic electromagnetic model in slab geometry. In the second place, in Sec. 3 we apply the procedure of Ref. [11] to this model in order to derive a two-moment Hamiltonian gyrofluid model for magnetic reconnection driven by electron inertia. This type of reconnection is believed to be relevant for sawtooth oscillations in tokamaks (see, e.g. Ref. [START_REF] Biskamp | Magnetic Reconnection in Plasmas[END_REF]). Because of the relative simplicity of the models under consideration, and of the limited number of fluid moments involved, the present results cannot be directly relevant for a realistic modelling of tokamak plasmas. They provide, however, an explicit example of a structure-preserving derivation which is a first step toward applications to more realistic higher-order fluid models. They also provide new knowledge about the Hamiltonian structures of the more basic gyrokinetic and gyrofluid models involved in the paper. The results are summarized in Sec. 4, where also limitations of the present analysis and possible future directions of investigation are discussed.

Parent gyrokinetic model and its Hamiltonian structure

In a Cartesian coordinate system (x, y, z) we consider a gyrokinetic model in slab geometry in the δf approximation. The model consists of the following evolution equations

∂g i ∂t + c B J 0 φ - v c A , g i + v ∂ ∂z g i + F i T i eJ 0 φ - v c A = 0, (1)
∂g e ∂t + c B φ - v c A, g e + v ∂ ∂z g e - F e T e e φ - v c A = 0, (2)
of the quasi-neutrality relation

e 2 n 0 T i (Γ 0 -1)φ + e dW(J 0 g i -g e ) = 0, (3)
and of Ampère's law

∇ 2 ⊥ A -4π e 2 c 2 n 0 Γ 0 M i + 1 M e A + 4π e c dWv(J 0 g i -g e ) = 0. (4)
Eqs. ( 1) and ( 2) govern the evolution of g i = f i + F i (e/T i )(v/c)J 0 A and g e = f e -F e (e/T e )(v/c)A, where f i,e (x, y, z, v, µ, t) indicate the perturbations of the ion and electron equilibrium distribution functions F i,e , respectively, which are defined as

F i (v, µ) = n 0 M i 2πT i 3/2 e - M i v 2 +2µB 2T i , F e (v, µ) = n 0 M e 2πT e 3/2 e -Mev 2 +2(Me/M i )µB 2Te . (5)
We indicate with v the velocity coordinate along the magnetic guide field direction (which we take to be the z direction) and with µ the ion magnetic moment. Ion and electron masses are denoted as M i and M e , respectively, whereas T i,e indicate constant temperatures characterizing the ion and electron equilibrium distribution functions. The constants c and e indicate the speed of light the proton charge, respectively. The flux function A(x, y, z, t) is related to the magnetic field by B = ∇A × ẑ + B ẑ, with B indicating the uniform and constant amplitude of the guide field component. The constant n 0 is the equilibrium density of both electron and ions, whereas φ(x, y, z, t) denotes the electrostatic potential. The symbol ∇ 2 ⊥ in Eq. ( 4) indicates the Laplacian operator restricted to the xy-plane. The integrations in Eqs. ( 3) and ( 4) are carried out over the element dW = 2πBdµ/M i . We assume that all fields are periodic along the x, y and z directions and that f i,e decay to zero sufficiently fast as v → ∞ and µ → +∞. The operators J 0 and Γ 0 are the standard operators appearing in gyrokinetic theory and correspond to the multiplication, in Fourier space, times the zeroth order Bessel function

J 0 (k ⊥ 2µB/M i ω 2
ci ) and times the function

I 0 (k 2 ⊥ ρ 2 i )e -k 2 ⊥ ρ 2 i , respectively, where k 2 ⊥ = k 2
x +k 2 y and where we indicated with ω ci = eB/M i c the ion cyclotron frequency, and with ρ i = T i /M i /ω ci the ion thermal gyroradius. Finally, the canonical bracket appearing in Eqs. ( 1)-( 2) is defined by [f, g] = ∂ x f ∂ y g -∂ y f ∂ x g, for two functions f and g.

In the model ( 1)-( 4) we neglected, as is customary by virtue of the small electron/ion mass ratio, the effects of the gyration radius for the electrons. However, the results present in the paper, can easily be extended to account for them as well.

Showing that the model ( 1)-( 4) possesses a Hamiltonian structure amounts to show that it can be cast in the form

∂χ ∂t = {χ, H}, (6)
where χ indicates a vector of dynamical variables, corresponding, in our case, to g i and g e . H = H(χ) is a functional, denoted as Hamiltonian, and {, } is a Poisson bracket, that is an antisymmetric bilinear operation satisfying the Leibniz identity {F G, H} = F {G, H} + G{F, H}, ∀F, G, H and the Jacobi identity

{{F, G}, H} + {{G, H}, F } + {{H, F }, G} = 0, ∀F, G, H. (7)
Because the Hamiltonian physically represents the conserved total energy of the system, a good candidate for H for our model is the conserved quantity [8]

H(g i , g e ) = 1 2 d 3 xdW T i F i g 2 i +
T e F e g 2 e + e(J 0 g i -g e )L -1

0 dW ′ (J 0 g i -g e ) -e v c (J 0 g i -g e )L -1 1 dW ′ (J 0 g i -g e ) . (8)
The expression (8) assumes the validity of the relations

φ = L -1 0 dW(J 0 g i -g e ), A = L -1 1 dWv(J 0 g i -g e ), (9)
that is the existence of invertible operators L 0,1 that permit, from Eqs. ( 3) and ( 4), to express φ and A in terms of g i and g e . In Fourier space it is not difficult to see that this inversion can be carried out for each k, apart from the case concerning Eq. ( 3) and modes with k x = k y = 0, for which Γ 0 -1 = 0. We restrict then to electrostatic potentials such that φ 0,0,kz = 0. The operators L -1 0 and L -1 1 also turn out to satisfy the symmetry property

d 3 xf (x)L -1 0,1 g(x) = d 3 xg(x)L -1 0,1 f (x)
and the operator J 0 is known (see, e.g. Ref. [5]) to satisfy the relation

d 3 xdWf (x, v, µ)J 0 g(x, v, µ) = d 3 xdWg(x, v, µ)J 0 f (x, v, µ)
. By virtue of these properties, we obtain that the functional derivatives of H with respect to g i,e can be written as

δH δg i = T i F i g i + eJ 0 φ - v c A , δH δg e = T e F e g e -e φ - v c A . (10)
Direct calculations show that, when using (8) as Hamiltonian, with the help of Eqs. (10) and of integration by parts, the Poisson bracket that yields Eqs. ( 1)-( 2), accounting for Eqs. ( 3) and ( 4), is given by

{F, G} = s=i,e d 3 xdW ∓ c eB g s [F gs , G gs ] -v F s T s F gs ∂ ∂z G gs , (11)
where the upper (lower) sign should be taken for ions (electrons). The subscripts on functionals, on the other hand, indicate functional derivatives, so that, for instance F gs = δF/δg s . The operation (11) is the direct sum of two independent Poisson brackets and consequently satisfies the properties of a Poisson bracket. We conclude then that the gyrokinetic model ( 1)-( 4) possesses a Hamiltonian structure. A similar Hamiltonian structure was presented in Ref. [START_REF] De Blank | Proc. of the 28th EPS Conference on Contr. Fusion and Plasmas Phys[END_REF] for the case of a drift-kinetic model.

Derivation of the Hamiltonian gyrofluid model

We define moments of g i,e as P i,e jk = dWv j µ k g i,e , where j and k are non-negative integers. Functionals F (g i , g e ) can be expressed as functionals of the infinite moments P i,e jk by using scalar invariance, so that F (g i , g e ) = F ({P i,e jk } j,k∈N ) and the functional derivatives consequently transform as F g i,e = j,k∈N v j µ k F i,e jk , where we denoted F i,e jk = δ F /δP i,e jk . Inserting this expression into (11) one formally obtains the Poisson bracket in terms of all the moments, which reads

{F, G} =

s=i,e m,n,j,k∈N

d 3 x ∓ c eB P s m+j n+k [F smn , G s jk ] -F smn ∂ ∂z G s jk dWv m+j+1 µ n+k F s T s , (12)
In practice, however, one is interested in deriving a Hamiltonian fluid model, describing the evolution of a finite number of moments, for instance the first M (N ) moments for the parallel (perpendicular) direction, for both the ions and the electrons, where M and N are two positive integers. This means that the functionals of interest will be of the kind F (P i 00 , P i 10 , P i 01 • • • P i M N , P e 00 , P e 10 , P e 01 • • • P e M N ) and will no longer depend on the infinite set of moments. The operation (12) between two such functionals, however, is not closed, in the sense that {F, G}, in general, will not be again a functional of the same kind, but it will depend also on moments of order higher than M (N ) in the parallel (perpendicular) direction. In order to obtain a closed system, the way we follow is that of modifying the Poisson bracket, imposing that such extra higher order moments be functions of the first M (N ) moments in the parallel (perpendicular) direction. Imposing such constraint in the bracket, however, can break the properties defining a Poisson bracket. In particular, the Jacobi identity can be violated by this operation. However, in some cases, which we refer to as to Hamiltonian closures, imposing the closure relation in the Poisson bracket, preserves the properties of the Poisson bracket. Hamiltonian fluid models can then be derived adopting Hamiltonian closures.

In particular, we consider here the case M = 1 and N = 0, that is we restrict to functionals of the first two moments in the parallel direction. The bracket (12) then becomes

{F, G} = s=i,e d 3 x ∓ c eB P s 0 [F s 0 , G s 0 ] ∓ c eB P s 1 ([F s 1 , G s 0 ] + [F s 0 , G s 1 ]) ∓ c eB P s 2 [F s 1 , G s 1 ] - n 0 M s F s 1 ∂ ∂z G s 0 + F s 0 ∂ ∂z G s 1 , (13)
where, in order to simplify the notation, given that the index for the powers of µ has been fixed equal to zero, we have defined P s j = P s j0 , for j = 0, 1, 2. It is evident that, as above anticipated, the result of the operation {F, G} in ( 13) depends also on P s 2 and consequently the resulting equations of motion cannot lead to a closed fluid model, for any Hamiltonian depending on P s 0 and P s 1 only. One remedy to this is to look for two functions P s (P s 0 , P s 1 ), if they exist, such that, when inserting P s 2 = P s into Eq. ( 13), the resulting operation (which will now be closed) still satisfies the Jacobi identity and can consequently still provide a Poisson bracket. It turns out that, for brackets of the form [START_REF] Biskamp | Magnetic Reconnection in Plasmas[END_REF], it has been shown [11] that the constraint of the Jacobi identity ( 7) selects, among all the possible closures (not considering, however, those involving integro/differential operators or an explicit dependence on the spatial coordinates), the following one:

P s 2 = a s P s 0 , for s = i, e, (14)
where a s are two arbitrary constants, which dimensionally must be homogeneous to a velocity squared. We remark that the Hamiltonian closure [START_REF] Waelbroeck | [END_REF] does not depend on first order parallel moments, which is a consequence of having chosen a Maxwellian with zero mean flow as equilibrium distribution function [11].

Once that a valid Poisson bracket, inherited from the original gyrokinetic bracket (11), has been found, in principle, Hamiltonian fluid models can be constructed by choosing the total energy functional (the Hamiltonian) H(P s 0 , P s 1 ) that one would like to be conserved, and then using the definition of Hamiltonian system (6) to get the four dynamical model equations, by replacing χ with the four moments P s 0 and P s 1 .

A relevant example of Hamiltonian closure, which is also consistent with the Hamiltonian of the original gyrokinetic model, is obtained by choosing a s = T s /M s . This corresponds to an isothermal closure, based on the equilibrium temperature, for the linearized pressure. The Poisson bracket then reads

{F, G} = s=i,e d 3 x ∓ c eB P s 0 [F s 0 , G s 0 ] ∓ c eB P s 1 ([F s 1 , G s 0 ] + [F s 0 , G s 1 ]) ∓ c eB T s M s P s 0 [F s 1 , G s 1 ] - n 0 M s F s 1 ∂ ∂z G s 0 + F s 0 ∂ ∂z G s 1 . (15)
Analogously to Ref. [8], one can then make use of the approximate decomposition

g s = F s P s 0 n 0 + M s T s v P s 1 n 0 , for s = i, e, (16)
which respects the moments definition as well as the isothermal closure relation. Inserting the expression (16) into Eq. ( 8) yields

H(P i 0 , P i 1 , P e 0 , P e 1 ) = 1 2 s=i,e d 3 x T s n 0 P 2 s 0 + M s P 2 s 1 n 0 + e 2 d 3 x(Γ 1 P i 0 -P e 0 )L -1 0 (Γ 1 P i 0 -P e 0 ) - e 2c d 3 x(Γ 1 P i 1 -P e 1 )L -1 1 (Γ 1 P i 1 -P e 1 ), (17)
where we made use of the common closure approximations (see, e.g. Ref. [8])

Γ 1 = Γ 1/2 0 , dWJ 0 g i = Γ 1 P i 0 and dWJ 0 vg i = Γ 1 P i 1 .
From such closure approximations it follows also that, using Eqs. ( 3) and ( 4), the following well-known relations hold:

∇ 2 ⊥ A = - 4πe c n 0 (Γ 1 u i -u e ), n 0 e Γ 0 -1 T i φ = -(Γ 1 n i -n e ), (18)
where we introduced the ion and electron density fluctuations n i,e = dWf i,e as well as the parallel fluid velocities u i,e = dWvf i,e /n 0 . Making use of Eqs. ( 18), the fluid Hamiltonian can be rewritten, in a physically more transparent way, as

H = 1 2 s=i,e d 3 x T s n 0 n 2 s + M s n 0 u 2 s - e 2 n 0 2T i d 3 xφ(Γ 0 -1)φ + 1 8π d 3 x|∇ ⊥ A| 2 . ( 19
)
The first and the second term in Eq. ( 19) account for the thermal and kinetic energy, whereas the third and the fourth term represent the electrostatic energy due to the polarization term and the magnetic energy, respectively. Due to the symmetry of the operators L -1 0,1 and Γ 1 involved in Eq. ( 17), one obtains the relations

δH δP i 0 = T i n 0 P i 0 + eΓ 1 φ, δH δP i 1 = M i n 0 P i 1 - e c Γ 1 A, (20)
From Eq. ( 6), using the Poisson bracket [START_REF] De Blank | Proc. of the 28th EPS Conference on Contr. Fusion and Plasmas Phys[END_REF], the Hamiltonian (17) and with the help of the relations (20)-( 21), one can finally obtain the Hamiltonian gyrofluid model

∂n i ∂t + c B [Γ 1 φ, n i ] - n 0 B [Γ 1 A, u i ] + n 0 ∂u i ∂z = 0, (22)
∂D ∂t + c B [Γ 1 φ, D] - T i M i B [Γ 1 A, n i ] + T i M i ∂n i ∂z + en 0 M i ∂Γ 1 φ ∂z = 0, (23)
where we introduced the quantities D = P i 1 = n 0 (u i + eΓ 1 A/M i c) and F = P e 1 = n 0 (u e -eA/M e c), corresponding to the first order fluid momenta and which represent the parallel canonical momenta of the ion and electron fluid, respectively. By multiplying Eq. ( 25) times M e one sees that, because of terms proportional to the electron inertia, the magnetic flux A is not frozen into any fluid. Indeed, in spite of the Hamiltonian character of the model, the presence of electron inertia violates the frozen-in condition, allowing for magnetic reconnection.

The model ( 22)-( 25) and its Hamiltonian structure were presented in Ref. [START_REF] Waelbroeck | [END_REF], but in that work the Hamiltonian structure was "guessed" a posteriori, after deriving the model from the more general gyrofluid model of Ref. [6]. The present derivation, on the other hand, originates from the gyrokinetic description, and guarantees, by construction, the Hamiltonian structure of the resulting model.

The above gyrofluid model is characterized by the presence of four Casimirs, that is functionals C, such that {C, F } = 0, for every functional F . From this definition and from Eq. ( 6) it follows that dC/dt = 0, so that Casimir are invariants of motion and consequently provide information about constraints on the dynamics, also in the nonlinear phase. The four Casimirs for the gyrofluid model can be written as d 3 xC s ± for s = i, e, where

C s ± = M s T s P s 1 ± P s 0 . (26)
The parent gyrokinetic bracket (11), on the other hand, possesses the two Casimirs d 3 xdWg s , for s = i, e. Note that the Casimirs of the gyrofluid model are "inherited" from those of the gyrokinetic model, in the sense that they can be obtained from those by using the decomposition ( 16) and replacing v with plus or minus the equilibrium thermal speed T i,e /M i,e . In the two-dimensional limit ∂/∂ z = 0 the Casimirs of both the gyrokinetic and the gyrofluid model extend to infinite families, given by d 3 xdWG s (g s ) and d 3 xC s ± (C s ± ), respectively, for s = i, e, where G s and C s ± are arbitrary functions. The relevance of the Casimirs of the gyrofluid model for two-dimensional reconnection has been shown in Ref. [START_REF] Comisso | [END_REF].

Conclusions

In this paper we have shown that the electromagnetic gyrokinetic model ( 1)-( 4) possesses a noncanonical Hamiltonian structure. From such model we have then derived a gyrofluid model, which, by construction, also possesses a Hamiltonian structure. The resulting model corresponds to a known model [START_REF] Waelbroeck | [END_REF] for magnetic reconnection. The derivation, however, is new and shows, by taking advantage of the results of Ref. [11], that the constraint of the Jacobi identity selects this model as a Hamiltonian model derived consistently from the Poisson bracket and the Hamiltonian of the parent gyrokinetic model. The Casimir invariants of the gyrokinetic parent model have been shown to be related to those of the final gyrofluid model, which were known to play a role in the cascade toward small scale characterizing collisionless reconnection. Clearly the present results possess evident limitations and they only represent a first step in a long term project aiming at deriving Hamiltonian closures for more realistic models. An obvious direction for improvement is that of including elements of tokamak geometry and physics in the Hamiltonian model. Obvious examples in this context would be toroidal geometry , magnetic, temperature and density equilibrium gradients (although the latter have already been partially treated in Ref. [12]). Perhaps even more fundamental is the problem of considering Hamiltonian closures at higher-order moments. Advanced fluid and gyrofluid models such as those of Refs. [1][2][3][4][5][6][7][8] reduced the gap between kinetic and fluid descriptions by evolving moments up to the heat fluxes, in both the parallel and perpendicular directions. A Hamiltonian derivation of three or four-moment models from kinetic theories is a natural step for the future, although the difficulty in identifying Hamiltonian closures might rapidly and nonlinearly increase when moving to higher order moments.

Acknowledgments

The author acknowledges useful discussions with P.J. Morrison and with the Nonlinear Dynamics team of the Centre de Physique Théorique. Financial support was received from the Agence Nationale de la Recherche (ANR GYPSI n. 2010 BLAN 941 03) and from the CNRS through the PEPS project GEOPLASMA2.


P Masci

O Bernard

F Grognard

E Latrille

J.-B B Sorba

J P Steyer

UR050 Lbe

Driving competition in a complex ecosystem: application to anaerobic digestion

Keywords: Anaerobic digestion, Bioreactors, Biological systems, Competition, Selection, Directed Evolution

Anaerobic digestion is a wastewater treatment process where bacteria degrade an organic substrate and produce methane, which can be used as a biofuel. The first task when starting up an anaerobic digester is the increase of its microbial population. It is a delicate phase, which is still not well understood, and its influence on the digester's future performance is not well known. During this phase, we show that a competition between the various species occurs and finally some species become dominant. In this paper, extending the competitive exclusion principle, we propose to drive the competition during this start-up phase, by regulating the volatile fatty acids concentration, with the aim of selecting species with good performance in the standard operating mode of the process. This new "selective" start-up strategy should lead to more efficient ecosystems.

I. INTRODUCTION AND MOTIVATION

Anaerobic digestion is a more and more widely used bioprocess for wastewater treatment. This complex ecosystem involves several hundreds of bacterial and achaeal species [START_REF] Delbès | Bacterial and archaeal 16s rdna and 16s rrna dynamics during an acetate crisis in an anaerobic digestor ecosystem[END_REF] that progressively degrade organic matter into methane and CO 2 . It has many advantages compared to the more widespread activated sludge process: it can handle concentrated substrates, produces few sludge, and methane can be recovered and used as a biofuel. However this process is difficult to manage since the steady state associated to the operating mode is not globally stable [START_REF] Hess | Design and study of a risk management criterion for an unstable anaerobic wastewater treatment process[END_REF]. As a corollary, the startup of the digester is a long and risky phase [START_REF] Wang | Accelerated start-up and enhanced granulation in upflow anaerobic sludge blanket reactors[END_REF] during which the digester loading is progressively increased in order to let the bacterial population adapt, grow and settle the reactor. Reactor startup is then a long procedure (from one month to almost one year) that is essential to achieve a high treatment capacity at steady-state [START_REF] Michaud | Methane yield as a monitoring parameter for the strat-up of anaerobic fixed film reactors[END_REF], and therefore this phase should be better understood and controlled.

Despite its key role, this phase did not receive so far much attention, perhaps because it is quite difficult to properly characterise since biomass measurements are not available. In [START_REF] Ramirez | Modeling and monitoring of microbial diversity in ecosystems -application to biological wastewater treatment processes[END_REF], it was shown that disturbing the microbial community in a "bang bang" like approach leads to better diversity and thus improves the bioprocess when facing a toxicant.

In this paper, we propose a new model for the start-up phase, including biodiversity. We assume that, at the beginning, N species are competing for the substrate. We show that, depending on the way the digester is started up, some populations will be enhanced and will preferentially remain in the digester. The underlying idea relies on the competitive exclusion principle, which has been theoretically studied and experimentally demonstrated in various conditions [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF], [START_REF] Hansen | Single-nutrient microbial competition[END_REF]. We exploit the idea proposed in [START_REF] Masci | Continuous selection of the fastest growing species in the chemostat[END_REF] to drive the competition during the start-up phase, by regulating the volatile fatty acids (VFA) concentration around a fixed value, which is determined by the normal operating mode of the digester. By doing so, we act upon the ecosystem in order to enhance the process efficiency. Also, bacterial biomass extinction, which is one of the main risks of the start-up phase, is avoided, and the obtained ecosystem presents an interesting steadiness property. Simulations and a first experiment are presented for supporting this start-up strategy.

II. PROCESS MODELLING

A. Presentation of the model

We propose here a simple macroscopic model (derived from [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]) in order to be able to handle the involved mathematics. We assume that, in a first acidogenesis step, the dissolved organic substrate, of concentration s 1 , is degraded by a population made of N a species of acidogenic bacteria (y j , with j ∈ [1; N a ]) into volatile fatty acids (VFA, denoted s 2 ). The growth rate of these bacteria is γ j (s 1 ):

k 1 s 1 γj (s1)yj -----→ y j + k 2 s 2 + k 5 CO 2
where k 1 , k 2 and k 5 are pseudo-stoichiometric coefficients which represent the transfer from substrate to acidogenic biomass, VFA and CO 2 . They are assumed to be constant between the various acidogenic species.

In a second step (methanogenesis), the VFA are degraded into CH 4 and CO 2 by N m methanogenic archaeabacteria (x i , with i ∈ [1; N m ]) with growth rate µ i (s 2 ):

k 3 s 2 µi(s2)xi -----→ x i + k 4 CH 4 + k 6 CO 2
where k 3 , k 4 and k 6 represent the transfer from substrate to methanogenic biomass, methane and CO 2 . They are also assumed to be constant between all the methanogenic species.

We assume that the concentrations are perfectly distributed throughout the reactor, and the dilution rate for the dissolved components is D. Since a part of the bacteria are attached on a fixed biofilm within the digester, the dilution rate for the y j and x i biomasses are β j D and α i D (i.e. the biomasses have a retention time 1/(β j D) or 1/(α i D)). The dynamical mass-balance model in a continuous stirred tank reactor is then straightforwardly derived [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF], [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]:

                 ṡ1 = D(s in 1 -s 1 ) -k 1 Na j=1 γ j (s 1 )y j ẏj = γ j (s 1 )y j -β j Dy j ṡ2 = D(s in 2 -s 2 ) + k 2 Na j=1 γ j (s 1 )y j -k 3 Nm i=1 µ i (s 2 )x i ẋi = µ i (s 2 )x i -α i Dx i q m = k 4 Nm i=1 µ i (s 2 )x i
(1) where s in 1 and s in 2 are respectively the concentration of the influent organic substrate and influent VFA, and q m is the methane flow rate.

We consider a Monod kinetics (2) for the substrate-saturated growth rate of acidogenic bacteria γ j (s 1 ) = γj

s 1 s 1 + h 1 j (2)
and an Haldane function [START_REF] Wang | Accelerated start-up and enhanced granulation in upflow anaerobic sludge blanket reactors[END_REF] for the methanogenesis to represent the possible inhibition by an accumulation of VFA [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]:

µ i (s 2 ) = μi s 2 s 2 + h 2 i + (s 2 /h 2 * i ) 2 (3)
where γ is the maximal growth rate of the acidogenic bacteria and μ is the potential maximum growth rates of the methanogenic bacteria, h 1 j and h 2 i the half-saturation constants associated to substrates s 1 and s 2 , and h 2 * i the inhibition constants associated to substrate s 2 .

III. THE "SELECTIVE" START-UP STRATEGY

A. Standard start-up strategy

The main danger during the start-up phase is imposing too high a dilution rate D, such that the VFA concentration becomes high and inhibitory, and all the methanogenic biomasses are washed out of the digester (see [START_REF] Hess | Design and study of a risk management criterion for an unstable anaerobic wastewater treatment process[END_REF] for a study of the two stable equilibria attraction bassins, in the mono-specific case). In order to avoid such a dramatic scenario, the plant operator often starts the digester with a very low dilution rate. Under such a control, the outcome of competition between bacterial species can be predicted in the case where ∀j ∈ {1,

• • • , N a }, β j = β and ∀i ∈ {1, • • • , N m }, α i = α for some α, β > 0.
Any acidogenic species whose maximum growth rate γj is lower than the dilution rate βD will be excluded from the digester, because for such a species ẏj < -K j y j , with K j = -γ j + βD > 0, and then y j (t) < y j (0)e -Kj t tends to zero. The same exclusion applies for all methanogenic species such that max s2 µ i (s 2 ) < αD.

For the other acidogenic species, it is possible to define equilibrium substrate concentration s * 1j such that γ j (s * 1j ) = βD where the s * 1j concentration depends on D:

∂ ∂D s * 1j ≥ 0
For the other methanogenic species, it is possible to define s * 2i and s

† 2i ≥ s * 2i such that µ i (s * 2i ) = µ i (s † 2i ) = αD where the s * 2i and s † 2i concentrations depend on D: ∂ ∂D s * 2i ≥ 0 and ∂ ∂D s † 2i ≤ 0
In order to prove the next property, we can then define a dilution rate D * low enough such that:

min i (s † 2i ) = max i (s * 2i )
and, for any S > 0, there exists

D S < D * such that ∀D < D S , min i (s † 2i ) > S > max i (s * 2i ) (4)
The following property determines the outcome of the competition when the dilution is small enough and when there is enough substrate in the input to avoid the straightforward wash-out of all species: Property: If ∀j, β j = β, ∀i, α i = α, the dilution satisfies D < D * , and s in 1 > min j (s * 1j ) then • the species with lowest s * 1j is selected among all acidogenic species (with lim t→∞ y j = y * j ), if, moreover, S = s in 2 + k2 k1 (s in 1 -s ⋆ 1j ) satisfies (4) • the species with lowest s * 2i is selected among all the methanogenic species. Proof : Under such a constant low dilution rate, the theory of competitive exclusion predicts (with Monod-like growth functions) that the most efficient species for the growth βD, i.e. the acidogenic species with lowest s * 1j , is selected [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF]. This is also true with Haldane-like growth functions when (4) is satisfied [START_REF] Wolkowicz | j.-h. accelerated start-up and enhanced granulation in upflow anaerobic sludge blanket reactors. water research 2004, 38, 2293[END_REF]: the methanogenic species with lowest s * 2i wins the competition and excludes all others.

Then if a constant dilution rate D < D * is imposed, two selection processes occur, which do not depend on the s in 1 and s in 2 input concentrations. Note that these two processes are independent, and that the trait which is optimized is the lowest substrate requirement. It is not possible to choose the s * 1j and s * 2i concentrations: selection occurs but it is not controlled. That is a major problem of this start-up strategy: species which would grow fast in the standard operating mode can be washed out; uncontrolled competition can lead to the selection of species with low growth rate in the standard operating mode, thus slowing the pollution removal process.

We propose an alternative strategy which permits to control the competition and select the most efficient species in the standard operating mode.

B. Objective of the proposed strategy

A reasonable objective for an anaerobic digester is to consume the pollutant substrate and obtain a constant VFA output concentration s2 (VFA is necessary for the process, but it is also a pollutant), with the highest possible dilution rate D, which corresponds to the wastewater treatment rate.

The objective of the strategy that we want to explore is to preserve and select the species whose equilibrium µ i (s 2 )/α i (denoted "relative growth rate") is maximal, thus optimizing the wastewater treatment rate (at equilibrium D = µ i (s 2 )/α i ). More than that, we propose to drive this selection since the beginning of the start-up phase, which may have a positive impact on the biofilm composition, thus resulting in better performances. For achieving this goal we regulate s 2 and show that it leads to the desired competition outcome.

C. Control design for the regulation of s 2

Here we assume that a control strategy has been set-up which ensures the convergence of s 2 towards the setpoint s2 < s in 2 . Because of space limitation, we won't detail the considered control law and neither prove its convergence. However we refer the reader to other works aiming at the regulation of substrates in an anaerobic digestion process [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF], [START_REF] Bernard | Modelling, monitoring and control of anaerobic digestion processes[END_REF], [START_REF] Mailleret | Robust nonlinear adaptive control for bioreactors with unknown kinetics[END_REF]. In the example section, the used controller is presented and its convergence is shown on a real plant.

D. Selection of archaeal species with maximal

µ i (s 2 )/α i
The principle of driving the competition for selecting species or individuals which maximize a chosen criterion is presented in [START_REF] Masci | Continuous selection of the fastest growing species in the chemostat[END_REF], where a general law is given for this selection, and specific selection criteria are given for Monod and Droop models. This theory can also be adapted and applied to system (1) if the VFA concentration is regulated to a constant value.

Main Theorem : In system (1), if a controller achieves the regulation of s 2 to s2 < s in 2 , then the methanogenic species x k with highest relative growth rate µ k (s 2 )/α k is selected and all the others are excluded.

Proof In a first step we will prove that the regulation of s 2 to s2 < s in 2 leads to the boundedness of the total methanogenic biomass Nm i=1 x i (t) (denoted X T (t)) in [X 0 ; X m ] with X 0 > 0. Therefore the total methanogenic biomass X T cannot diverge, and its complete washout cannot occur (at least one methanogenic species will remain in the chemostat). In a second step, we will show that the difference in relative velocity between species lead to the survival of only one species, with maximal relative growth rate.

First, s 2 regulation causes the total methanogenic biomass boundedness. Indeed, defining the variable z = s 1 + k1 k2 s 2 + k3k1 k2 X T and z in = s in 1 + k1 k2 s in 2 , we have

ż = Dz in -Ds 1 - k 1 k 2 Ds 2 -D k 3 k 1 k 2 i α i x i so that D(z in -z) ≤ ż ≤ D(z in -αz) (5)
where α = min i α i . We directly conclude that an upperbound on X T is deduced from an upperbound on z, so that

X m = k 2 k 3 k 1 max z in α , z (0)
We can deduce from this upperbound that ∞ 0 D(τ )dτ is unbounded. Indeed, if it was not the case and knowing that s 2 converges to s2 , any ẋi equation aymptotically yields ẋi = (µ(s 2 ) -D(τ ))x i whose solution is unbounded when the aforementionned integral is bounded, which is in contradiction with the existence of X m .

Also, the unboundedness of ∞ 0 D(τ )dτ and (5) yield lim inf t→∞ z ≥ z in , so that, having lim t→∞ s 2 = s2 and using the fact that s 1 (t) ≤ s in 1 for all times

lim inf t→∞ X T ≥ 1 k 3 (s in 2 -s2 )
so that there is a lower-bound X 0 on X T . Secondly, X T being lower and upper bounded, let us show that the difference in relative growth rates between the species lead to competitive exclusion. We denote x k the species with maximal µ k (s 2 )/α k , such that µ k (s 2 )/α k > µ i (s 2 )/α i , ∀i = k and t s a time such that µ k (s 2 (t))/α k > µ i (s 2 (t))/α i , ∀i = k, ∀t ≥ t s . Let us define

r i = ln x 1/α k k x 1/αi i which leads to ṙi (t) = ẋk α k x k - ẋi α i x i = µ k (s 2 (t)) α k - µ i (s 2 (t)) α i > 0
for all t ≥ t s . Then r i (t) will increase, so that lim t→+∞ r i (t) = +∞. This means that the ratio

x 1/α k k x 1/α i i
will tend to infinity.

As x k is upper bounded by X m , it cannot diverge, which implies that

lim t→+∞ x i (t) = 0, ∀i = k
As x k is also lower bounded by X 0 , then lim inf t→+∞ x k (t) > 0 and we obtain that only methanogenic species x k is not excluded from the chemostat.

Thus, if one can regulate the VFA concentration s 2 so that it remains constant at a given s2 value, then the selection of the species with maximal relative growth rate µ k (s 2 )/α k will occur, thus shaping the ecosystem and optimizing the wastewater treatment rate at equilibrium D = µ k (s 2 )/α k .

IV. SIMULATIONS

Even though, as previously said, several hundreds of microbial species are present in anaerobic digesters, we have simulated system (1) with only 3 acidogenic species and 4 methanogenic species to keep the system understandable. Species-specific parameters are taken from [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF] with ±50% values. The substrate and VFA inputs were s in 1 = 5g/L and s in 2 = 5mmol/L. The relative growth functions µ i (s 2 )/α i of the various methanogenic species is presented in Fig. 1. Species 1 grows better with low VFA concentration, species 3 grows faster under high VFA concentration, and species 2 is better for medium ones. Species 4 is less efficient than the others for any VFA concentration, so that it should be excluded whatever the start-up strategy. Relative growth functions µ i (s 2 )/α i of the 4 methanogenic species used in simulation. Species 1 (Solid), species 2 (Dashed), species 3 (Thick), species 4 (Dashdot)

A. Interest of the "selective" start-up strategy

In this section we show the interest of the selective start-up strategy, by comparing it with two other strategies.

The first strategy is the standard start-up strategy where D is chosen low (D = 0.05) to avoid bacterial wash out. The result is presented in Fig. 2, where we can see that species 1, which grows slowly for high VFA concentrations, is selected. The digester was started safely but the obtained ecosystem is such that the digester's performance will not be optimal. Fig. 2. Standard start-up strategy with low D = 0.05, to avoid archaeal washout. This strategy leads to the selection of species 1, which is not optimal for VFA concentrations higher than 2mmol/L (see Fig. 1)

The second strategy is the selective start-up strategy where the species which grows faster at VFA concentration s2 is favored. In Fig. 3 we use a controller to regulate s 2 to s2 = 3mmol/L, and species 2 is selected. We can see on Fig. 1 that it is the most efficient species at concentration 3mmol/L. For this strategy we used the following pseudolinearizing controller

D = λ 1 (s 2 -s 2 ) + k3 k4 q m s in 2 -s 2 + λ 2 t 0 (s 2 -s 2 (s))ds (6)
(λ 1 and λ 2 are two gains) where the non-measured acidogenic activity is replaced by a PI correction.

The last one is a naive start-up strategy where the value of D is constant and equal to µ 2 (s 2 ). The objective of this strategy is to obtain the same result as in the selective strategy, but the simulation of Fig. 4 shows that it leads to process failure: the VFA were accumulated such that they became inhibitory (≈ 100mmol/L), and the methanogenic biomasses were washed out. These simulations show that the selective startup strategy avoids bacterial extinction, while selecting efficient species for the digester in the standard operating mode.

B. Selecting a population optimal for a given VFA concentration

Let us suppose now that we want an anaerobic digester whose VFA concentration in the standard operating mode is s2 = 8mmol/L. Using linearizing controller (6) to regulate s 2 to s2 , we obtain the simulation of Fig. 5 where species 3 (the most appropriate at concentration 8mmol/L) is selected. This simulation emphasizes that depending on the digester VFA concentration in the standard operating mode, different species should be selected during the start-up phase, so that the obtained ecosystem will be more efficient.

C. The selective process consolidates the ecosystem

On Fig. 6 we can see that the dilution rate, during the s2 selective start-up, rises until selection occurs and then becomes almost constant and equal to µ 2 (s 2 )/α i . Then, after time t = 100days, the controller is turned off and the dilution rate is kept constant with no consequence for the digester, as the VFA concentration stays almost constant: the ecosystem is settled and composed of the most efficient population in the standard operating mode, such that it became naturally stable around the desired equilibrium (s 2 = s2 ). This result is of great importance and must be verified experimentally : the selective start-up strategy gives birth to a stable ecosystem involving bacterial species such that the natural steady state is exactly the requested one, without needing a control feedback anymore.

V. APPLICATION TO A REAL ANAEROBIC DIGESTER START-UP PHASE

In order to verify the proposed selection principle and to assess the effect of a controller that regulates the VFA during the start-up phase, experiments have been carried out at the LBE-INRA Laboratory in Narbonne (France).

A. Experiment design

The process is an up-flow anaerobic fixed bed reactor with a useful volume of 0.548 m3. The reactor is highly instrumented and many variables were measured during the experiments [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]. More details about the process and evaluation of its online instrumentation are available in [START_REF] Steyer | Evaluation of a four year experience with a fully instrumented anaerobic digestion process[END_REF].

B. Experimental results and discussion

Control (6) was applied, with command s2 = 4 g/L. The experimental results are presented in Fig. [START_REF] Hansen | Single-nutrient microbial competition[END_REF]. We can see on this figure that our control objective was attained : the VFA concentration was well driven towards s2 despite many disturbances related to failures and disturbances in the process. It is worth noting that the more the ecosystem is selected, the more efficient the regulation. At the end, when the interesting species have probably been selected in the digester, the regulation is more and more efficient.

Finally, the proposed strategy enabled us to reach a high dilution rate D = 0.78 while avoiding start up failure.

VI. CONCLUSIONS

We proposed a start-up strategy for an anaerobic digester, whose aim is to influence the structure of the ecosystem, so that the obtained ecosystem in the end of start-up will be optimal for a given VFA concentration. For implementing such a selective process, we drive the competition by regulating the VFA concentration at a fixed value. We showed how the obtained ecosystem should also be consolidated because of the selection that occurs during the start-up phase.

According to model (1), a control law was proposed and a mathematical demonstration was given for predicting the competition outcome. Simulations were used to verify and explain the utility of the method, and a real experiment was presented where the start-up strategy was tested.

We emphasize that the aim of this strategy is not to minimize the start-up duration, but rather to build an efficient ecosystem. This kind of startup strategy, where the objective is to mould the structure of the bacterial community for optimizing a particular process, was never (to our knowledge) proposed before for anaerobic digestion, and further investigations in this direction are presently performed.

Fig. 1 .

 1 Fig. 1.Relative growth functions µ i (s 2 )/α i of the 4 methanogenic species used in simulation. Species 1 (Solid), species 2 (Dashed), species 3 (Thick), species 4 (Dashdot)

Fig. 3 .

 3 Fig. 3. Selective start-up strategy : species 2, which is the most efficient under VFA concentration s2 = 3mmol/L, is selected.

Fig. 4 .

 4 Fig. 4. Naive start-up strategy : by choosing D high (D = µ 2 (s 2 )) since the beginning of start-up, the methanogenic biomasses became inhibited and were washed out.

Fig. 5 .

 5 Fig. 5. Start-up selection, under the closed loop control, of species which maximize µ 2 (s 2 )/α i , with s2 = 8mmol/L. Depending on the VFA concentration in the standard operating mode, different populations should be selected to optimize the process.

Fig. 6 .

 6 Fig. 6. The dilution rate (top) in the s2 experiment was kept constant after time t = 100days with negligible consequence on the VFA concentration (bottom), showing that selection leads naturally to a consolidated ecosystem where the desired equilibrium (s 2 = s2 ) is stable.

Fig. 7 .

 7 Fig. 7. Start-up experiment where s 2 is regulated at s2 = 4g/L. The influent flow rate (Q in = D/V ) and the VFA concentration are plotted.

VII. ACKNOWLEDGMENTS

This work has been carried out with the support provided by the CODA ARC project funded by the INRIA.