Tahani Mtar

Radhouane Fekih-Salem
email: radhouene.fekihsalem@isima.rnu.tn

Tewfik Sari

Bruce Watson

Eric Badouel

Eds. Oumar Niang

Effect of the mortality on a density-dependent model with a predator-prey relationship

Keywords: Coexistence, Limit cycle, Predator-prey relationship, Supercritical Hopf bifurcation MOTS-CLÉS : Coexistence, Cycle limite, Relation proie-prédateur, Bifurcation de Hopf super-critique

In this work, we consider an inter-specific density-dependent model of two species competing on a single nutrient in a chemostat, taking into account the prey-predator relationship. With different dilution rates, we give a complete analysis of the existence and local stability of all the steady states. Indeed, under the joint effect of competition with the prey-predator relationship and mortality, we show that the positive steady state can be destabilized by a supercritical Hopf bifurcation with the appearance of a stable limit cycle. However, without mortality, there is a stable persistence of the coexistence steady state where there cannot be periodic orbits.

RÉSUMÉ.

Dans ce travail, nous considérons un modèle densité-dépendant inter-spécifique de deux espèces en compétition sur un seul nutriment dans un chémostat, en tenant compte de la relation proie-prédateur. Avec des taux de dilution différents, nous donnons une analyse complète de l'existence et de la stabilité locale de tous les points d'équilibre. En effet, sous l'effet joint de la compétition avec la relation proie-prédateur et de la mortalité, nous montrons que l'équilibre positif peut se déstabiliser par une bifurcation de Hopf supercritiques avec l'apparition d'un cycle limite stable. Cependant, sans mortalité, il y a persistance stable de l'équilibre de coexistence où il ne peut pas y avoir des orbites périodiques.

Introduction

In the chemostat, the coexistence of species could be explained through several mechanisms, such as the intra and interspecific competition between the species [START_REF] Borsali | Contribution to the study of the effect of the interspecificity on a two nutrients competition model[END_REF][START_REF] Fekih-Salem | Sur un modèle de compétition et de coexistence dans le chémostat[END_REF], the flocculation of the species [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF][START_REF] Fekih-Salem | Emergence of coexistence and limit cycles in the chemostat model with flocculation for a general class of functional responses[END_REF], and the density-dependence of the growth functions [START_REF] Fekih-Salem | A density-dependent model of competition for one resource in the chemostat[END_REF][START_REF] Lobry | A new hypothesis to explain the coexistence of n species in the presence of a single resource[END_REF]. Also, predator-prey interaction can prove that species can coexist in the chemostat, see [START_REF] Boer | Food chain dynamics in the chemostat[END_REF][START_REF] Li | Simple food chain in a chemostat with distinct removal rates[END_REF]. Here, we will focus on the predator-prey relationship which is characterized by the fact that the one species (the prey) promotes the growth of the other species (the predator) which in turn inhibits the growth of the prey species. In this context, in [START_REF] Hajji | How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?[END_REF][START_REF] Mtar | Interspecific density-dependent model of predatorprey relationship in the chemostat[END_REF], it was considered an interspecific density-dependent model in the chemostat with the same dilution rate, taking into account the predator-prey relationship, when it was proved that the system may exhibit the coexistence and the bistability with a multiplicity of positive steady states. Moreover, in [START_REF] Mtar | Interspecific density-dependent model of predatorprey relationship in the chemostat[END_REF], the operating diagram according to the control parameters S in and D shows all steady states that can appear or disappear only through saddle-node or transcritical bifurcations. In the present work, we proposed to study an extension of the model studied in [START_REF] Hajji | How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?[END_REF][START_REF] Mtar | Interspecific density-dependent model of predatorprey relationship in the chemostat[END_REF], considering the predator-prey relationship and distinct dilution rates. The model takes the form    Ṡ = D(S in -S) -f 1 (S, x 2 )x 1 -f 2 (S, x 1 )x 2 , ẋ1 = (f 1 (S, x 2 ) -D 1 )x 1 , ẋ2 = (f 2 (S, x 1 ) -D 2 )x 2 , [START_REF] Boer | Food chain dynamics in the chemostat[END_REF] where S(t) is the concentration of the substrate at time t; x 1 (t) and x 2 (t) are, respectively, the concentrations of prey and predator at time t; D and S in are, respectively, the dilution rate and the concentration of the substrate in the feed device; f 1 represents the growth rate of the prey which depends only on S and x 2 and f 2 represents the growth rate of the predator which depends only on S and x 1 ; D i represents the removal rate of the species x i and can be modeled by

D i = D + a i , i = 1, 2
where a i represents the nonnegative mortality rate parameter of species x i . In the existing literature, it is well known that the addition of mortality terms of the species in a predatorprey model in the chemostat can cause destabilization of the system where stable limit cycles [START_REF] Boer | Food chain dynamics in the chemostat[END_REF][START_REF] Kooi | Chaotic behaviour of a predator-prey system in the chemostat[END_REF] and multiple chaotic attractors [START_REF] Kooi | Chaotic behaviour of a predator-prey system in the chemostat[END_REF] are found.

In our knowledge, model [START_REF] Boer | Food chain dynamics in the chemostat[END_REF] has not yet been studied in the literature. Thus, our study, in this paper, is the first attempt to analyze the effect of the mortality on the behavior of a predator-prey model in the chemostat with interspecific density-dependent growth rates. Using the nullcline method, we present a geometric characterization that describes all steady states of the model (1) and shows the multiplicity of positive steady states. Considering specific growth rates, we succeeded in finding a set of parameters such that the coexistence can hold around a stable limit cycle, in contrast to the case without mortality where the coexistence may occur only around a positive steady state (see [START_REF] Hajji | How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?[END_REF][START_REF] Mtar | Interspecific density-dependent model of predatorprey relationship in the chemostat[END_REF]).

This paper is organized as follows. The next section presents general assumptions for the growth functions of the model [START_REF] Boer | Food chain dynamics in the chemostat[END_REF] and the analysis of the existence of steady states. In Section 3, the asymptotic behavior analysis of model (1) was done. Considering specific growth rates, numerical simulations are presented in Section 4 to show the emergence of a limit cycle through a supercritical Hopf bifurcation. Finally, conclusion is drawn in Section 5. Most of the proofs are reported in Appendix A and the maximal number of solutions of an equation to determine the nullcline in this particular case with specific growth rates is given in Appendix B. Finally, all the parameter values used in simulations are provided in Appendix C.

Assumptions on the model and steady states

In what follows, we study model (1) using the following general assumptions on the growth rates f i , for i = 1, 2, j = 1, 2, i = j : (H0) f i : R 2 + -→ R + is continuously differentiable. (H1) f i (0, x j ) = 0, for all x j ≥ 0.

(H2)

∂f i ∂S (S, x j ) > 0, for all S ≥ 0, x 1 > 0 and x 2 ≥ 0.

(H3) ∂f 1 ∂x 2 (S, x 2 ) < 0 and ∂f 2 ∂x 1 (S, x 1 ) > 0, for all S > 0, x 1 ≥ 0 and x 2 ≥ 0.

(H4) f 2 (S, 0) = 0, for all S > 0. Assumption (H1) means that the substrate is necessary for the growth of the two species. Assumption (H2) means that the growth rate of each species increases with the concentration of the substrate. Assumption (H3) means that the growth of the first species x 1 is inhibited by the second species x 2 , while the growth of second species x 2 increases with the presence of first species x 1 . Assumption (H4) means that the prey species x 1 is necessary for the growth of the predator species x 2 . One has the following preliminary result on positivity and boundedness of solutions.

Proposition 2.1 For any non-negative initial condition, the solutions of system (1) remain non-negative and positively bounded. In addition, the set

Ω = (S, x 1 , x 2 ) ∈ R 3 + : S + x 1 + x 2 ≤ S in
is positively invariant and is a global attractor for the dynamics [START_REF] Boer | Food chain dynamics in the chemostat[END_REF].

The first step is to determine the steady states of (1). A steady state of (1) must be a solution of the system

     0 = D(S in -S) -f 1 (S, x 2 )x 1 -f 2 (S, x 1 )x 2 , 0 = (f 1 (S, x 2 ) -D 1 )x 1 , 0 = (f 2 (S, x 1 ) -D 2 )x 2 .
(

) 2
By (H4) and from the third equation of (2), when x 1 = 0, it follows that x 2 = 0, that is, system 1 cannot have a steady state where only the predator exists. Therefore, system (1) has the following types of steady states:

• the washout E 0 = (S in , 0, 0), that always exists, where both populations are extinct.

• the extinction of the predator E 1 = ( S, x1 , 0), where second population x 2 is extinct.

• the coexistence steady state E * = (S * , x * 1 , x * 2 ), where both populations survive. The components S = S and x 1 = x1 of a boundary steady state E 1 are the solutions of (2) with x 1 > 0 and x 2 = 0. Therefore, S and x1 are the solutions of equations

D(S in -S) = D 1 x1 , (3) f 1 ( S, 0) = D 1 . (4)
From ( 3), we obtain

S = S in - D 1 D x1 . (5)
Replacing S by this expression in (4), we see that x 1 = x1 must be a solution of

f 1 S in - D 1 D x 1 , 0 = D 1 . (6)
S is positive if, and only if, x1 < DS in /D 1 , that is to say, ( 6) has a solution in the interval ]0,

DS in /D 1 [. The function x 1 → f 1 (S in -D 1 x 1 /D, 0) is decreasing from f 1 (S in , 0) for x 1 = 0 to 0 for x 1 = DS in /D 1 .
Thus, there exists a solution 6) if, and only if,

x 1 = x1 ∈ ]0, DS in /D 1 [ satisfying (
f 1 (S in , 0) > D 1 . (7)
If such an x1 exists then it is unique. Therefore, we can state the following result:

Proposition 2.2 Under assumptions (H0)-(H3), the boundary steady state E 1 of (1) with x1 > 0 and x 2 = 0 exists if, and only if, condition [START_REF] Fekih-Salem | Sur un modèle de compétition et de coexistence dans le chémostat[END_REF] holds where x1 is the solution of ( 6) and S is given by [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF]. If it exists, then it is unique.

The components S = S * , x 1 = x * 1 and x 2 = x * 2 of a positive steady state E * must be the solutions of (2) with x 1 > 0 and x 2 > 0, that is, S * , x * 1 and x * 2 are the solutions of the set of equations

D(S in -S) = D 1 x 1 + D 2 x 2 , (8) f 1 (S, x 2 ) = D 1 , (9)
f 2 (S, x 1 ) = D 2 . (10)
From ( 8), it follows that S * is given by

S * = S in - D 1 D x 1 - D 2 D x 2 . ( 11
)
Replacing S * by this expression in equations ( 9) and ( 10), we see that

(x 1 = x * 1 , x 2 = x * 2 ) must be a solution of f 1 (x 1 , x 2 ) = 0, f 2 (x 1 , x 2 ) = 0, (12)
where

f i (x 1 , x 2 ) := f i (S in - D 1 D x 1 - D 2 D x 2 , x j ) -D i , i = 1, 2, j = 1, 2, i = j. ( 13
)
One see that S * is positive if, and only if,

D 1 x * 1 /D + D 2 x * 2 /D < S in
, that is to say, (12) has a positive solution in the set

M := (x 1 , x 2 ) ∈ R 2 + : D 1 D x 1 + D 2 D x 2 < S in .
To solve (12) in M , we need the following results:

Lemma 2.1 Assume that ( 7) holds. The equation f 1 (x 1 , x 2 ) = 0 defines a smooth decreasing function

F 1 : [0, x1 ] -→ [0, x2 ] x 1 -→ F 1 (x 1 ) = x 2 , such that F 1 (x 1 ) = 0, F 1 (0) = x2 and -D 1 /D 2 < F 1 (x 1 ) < 0 for all x 1 ∈ [0, x1 ],
where x1 is the solution of ( 6) defined by f 1 (x 1 , 0) = 0 and x2 is the solution of the equation f 1 (0, x 2 ) = 0. Furthermore, the graph To determine the function defined by the equation f 2 (x 1 , x 2 ) = 0, we will see the number of solutions of the equation f 2 (x 1 , 0) = 0 in the following result.

γ 1 of F 1 lies in M (see Fig. 1(b)). More precisely, (x 1 , F 1 (x 1 )) ∈ M for all x 1 ∈ [0, x1 ].
Lemma 2.2 The equation f 2 (x 1 , 0) = 0 has a solution in [0, DS in /D 1 ] if, and only if, max x1∈[0,DSin/D1] f 2 S in - D 1 D x 1 , x 1 ≥ D 2 . ( 14
)
Generically, this equation admits an even number of solutions in [0,

DS in /D 1 ].
In what follows, without loss of generality, we can assume that the equation f 2 (x 1 , 0) = 0 admits two solutions x 1 1 and x 2 1 in [0, DS in /D 1 ] (see Fig. 1(a,c)). Appendix B shows that this property is fulfilled when the function x 1 → f 2 (x 1 , 0) is deduced from the specific growth rate (21) satisfying our hypotheses. For more general multimodal functions, the study can be treated similarly. Then, we can state the following result:

Lemma 2.3 Under condition (14), the equation f 2 (x 1 , x 2 ) = 0 defines a function F 2 : [x 1 1 , x 2 1 ] -→ [0, DS in /D 2 [ x 1 -→ F 2 (x 1 ) = x 2 , such that F 2 (x 1 1 ) = F 2 (x 2 1 ) = 0 and -D 1 /D 2 < F 2 (x 1 ) for all x 1 ∈ [x 1 1 , x 2 1 ]
, where x 1 1 and x 2 1 are the solutions of the equation f 2 (x 1 , 0) = 0. Moreover, the graph γ 2 of F 2 lies in M and (x 1 , F 2 (x 1 )) ∈ M for all x 1 ∈ [x 1 1 , x 2 1 ] (see Fig. 1(c)). The following result is a consequence of the previous lemmas.

Proposition 2.3 A positive steady state E * = (S * , x * 1 , x * 2 )
of (1) exists if, and only if, the curves γ 1 and γ 2 have a positive intersection, where S * is given by ( 11) and (x * 1 , x * 2 ) is a positive solution of equations

x 2 = F 1 (x 1 ) and x 2 = F 2 (x 1 ). ( 15
)
Note that x1 , x 1 1 and x 2 1 represent the coordinates of the intersections of the curves γ 1 and γ 2 with the axis x 2 = 0. To determine the number of positive steady states, three cases must be distinguished according to their relative positions:

Case 1 : x1 < x 1 1 < x 2 1 , Case 2 : x 1 1 < x1 < x 2 1 , Case 3 : x 1 1 < x 2 1 < x1 . ( 16
)
Note that Case 1 or Case 3 is equivalent to the condition f 2 (x 1 , 0) < 0, while Case 2 is equivalent to f 2 (x 1 , 0) > 0 (see Fig. 1(a)). The number of positive steady states of ( 1) is given by the following proposition.

Proposition 2.4 Assume that ( 7) and ( 14) hold.

• If Case 1 holds, then the system (1) can have no positive steady state.

• If Case 2 holds, then there exists at least one positive steady state. Generically, there is an odd number of positive steady states.

• If Case 3 holds, then the system can have no positive steady state or generically an even number of positive steady states.

Local stability of steady states

In this section, we focus on the study of local asymptotic stability of each steady state of system [START_REF] Boer | Food chain dynamics in the chemostat[END_REF]. For convenience, we shall use the abbreviation LES for Locally Exponentially Stable steady state. It is useful to use the change of variables z = S + x 1 + x 2 . Using the variables (z, x 1 , x 2 ), system (1) can be written

   ż = D(S in -z) -a 1 x 1 -a 2 x 2 , ẋ1 = (f 1 (z -x 1 -x 2 , x 2 ) -D 1 )x 1 , ẋ2 = (f 2 (z -x 1 -x 2 , x 1 ) -D 2 )x 2 .
(17) Let J denote the Jacobian matrix of (17) evaluated at (z, x 1 , x 2 ):

J =   -D -a 1 -a 2 Ex 1 f 1 -D 1 -Ex 1 -(E + G)x 1 F x 2 (H -F )x 2 f 2 -D 2 -F x 2   (18)
where

E = ∂f 1 ∂S , F = ∂f 2 ∂S , G = - ∂f 1 ∂x 2 , H = ∂f 2 ∂x 1 . (19)
which are nonnegative. The stability of the boundary steady states E 0 and E 1 is given as follows.

Proposition 3.1

• E 0 is LES if,
J * =   -D -a 1 -a 2 Ex * 1 -Ex * 1 -(E + G)x * 1 F x * 2 (H -F )x * 2 -F x * 2  
where E, F , G and H are given by ( 19) and evaluated at E * . The characteristic polynomial is given by:

P (λ) = -λ 3 + c 1 λ 2 + c 2 λ + c 3 ,
where

c 1 = -(D + Ex * 1 + F x * 2 ), c 2 = -D 1 Ex * 1 -D 2 F x * 2 + (F G -GH -EH)x * 1 x * 2 , c 3 = (D 1 GF -DGH -D 2 EH)x * 1 x * 2 .
Since c 1 < 0, according to the Routh-Hurwitz criterion, E * is LES if, and only if,

c 3 < 0 and c 4 = c 1 c 2 + c 3 > 0. (20)
In the next proposition, we will show that the sign of c 3 is given by the position of curves of functions F 1 (•) and F 2 (•).

Proposition 3.2 One has c 3 = D 2 F D2 D E + G (F 1 (x * 1 ) -F 2 (x * 1 ))x * 1 x * 2 .
Since the condition c 4 > 0 of the Routh-Hurwitz criterion (20) could be unfulfilled where c 4 can change its sign by varying the control parameter S in , the next section is devoted to show numerically how the positive steady state E * could be destabilized with the emergence of stable limit cycles.

Numerical simulations

In order to show that the condition c 4 > 0 evaluated at E * = (S * , x * 1 , x * 2 ) could be unfulfilled, all biological parameters were fixed. Then, the control parameter S in was varied. To see the change of sign of the function S in → c 4 (S in ) evaluated at E * , it is useful to illustrate the curve of this function. To this end, we assumed that the growth functions satisfying the assumptions (H0)-(H4) are given by :

f 1 (S, x 2 ) = m 1 S K 1 + S 1 1 + x2 Ki , f 2 (S, x 1 ) = m 2 S K 2 + S x 1 L 2 + x 1 , (21)
where m i , K i , i = 1, 2, and L 2 are the maximum growth rates and the Michaelis-Menten constants, respectively, while K i is the inhibition factor due to x 2 for the growth of the species x 1 . Indeed, we succeeded in finding a set of parameters (see Table 1) such that c 4 can change its sign as S in increases (see Fig. 2(a)). In fact, the solution S cr in 6.265 of equation c 4 (S in ) = 0 represents the critical value of S in for which the positive steady state loses its stability.

To understand and analyze this change of local behavior of E * in S cr in without any bifurcation with other steady states, we determine numerically the eigenvalues of the Jacobian matrix J * at the positive steady state E * . Indeed, this Jacobian matrix J * has one negative eigenvalue and one pair of complex-conjugate eigenvalues :

λj = µ(S in ) ± iν(S in ), j = 1, 2.
Increasing the control parameter S in , Fig. 2(b) shows that this pair crosses the imaginary axis at the critical value S in = S cr in from negative half plane to positive half plane, that is, it becomes purely imaginary for S cr in such that µ(S cr in ) = 0, with ν(S cr in ) = 0. Thus, E * becomes unstable with occurrence of stable limit cycle via a supercritical Hopf bifurcation. For S in 6.27 > S cr in , the numerical simulations (see Fig. 2(c)) show the occurrence of limit cycle when the oscillations are sustained.

Conclusion

In this work, we have analyzed the effect of the mortality of species on the behavior of an interspecific density-dependent model with a predator-prey relationship. We give a complete analysis of the existence and local stability of all steady states where we have shown the multiplicity of positive steady states by using the nullcline method. Considering specific growth rates, we succeeded in finding a set of parameters such that a condition of the Routh Hurwitz criterion (20) can change its sign as S in increases. More precisely, the positive steady state E * can destabilize through a supercritical Hopf bifurcation with the occurrence of a stable limit cycle. This feature cannot occur without mortality, where generically the coexistence can be only around a positive steady state. is a solution [START_REF] Boer | Food chain dynamics in the chemostat[END_REF]. By uniqueness of solutions, the set Γ 0 cannot be reached in finite time by trajectories for which x 1 (0) > 0 or x 2 (0) > 0. Furthermore, for i = 1, 2, if x i = 0 then ẋi = 0 and the set

Γ i = (S, x 1 , x 2 ) ∈ R 3 + : S ≥ 0, x i = 0, x j ≥ 0 , j = 1, 2, i = j
is invariant under the system (1) because the function

t -→ S in + (S(0) -S in )e -Dt , 0, x 2 (0)e -D2t (resp. t -→ (S(t), x 1 (t), 0))
is a solution of [START_REF] Boer | Food chain dynamics in the chemostat[END_REF]. By uniqueness of solutions, the set Γ i cannot be reached in finite time by trajectories for x i (0) > 0. Therefore, the solutions remain non-negative. Let z = S + x 1 + x 2 . From system (1), we have

ż = D(S in -z) -a 1 x 1 -a 2 x 2 .
Consequently, ż ≤ D(S in -z).

Using Gronwall's lemma, we obtain z(t) ≤ S in + (z(0) -S in )e -Dt for all t ≥ 0.

We deduce that z(t) ≤ max(z(0), S in ) for all t ≥ 0.

Therefore, the solutions of (1) are positively bounded and are defined for all t ≥ 0. From (22), we deduce that the set Ω is positively invariant and is a global attractor for (1).

Proof of Lemma 2.1. Let l 1 be a fixed line defined by x 1 = b 1 that intersects the line δ defined by

D 1 D x 1 + D 2 D x 2 = S in , at point x 2 = c 2 = D D 2 S in - D 1 D 2 b 1 ,
(see Fig. 1(b)). From hypotheses (H1)-(H3), it follows that the function

x 2 → f 1 (b 1 , x 2 ) is decreasing from f 1 (S in -D 1 b 1 /D, 0)-D 1 for x 2 = 0 to -D 1 for x 2 = c 2 . Therefore, there exists a unique b 2 ∈ [0, c 2 [ such that f 1 (b 1 , b 2 ) = 0 if, and only if, f 1 S in - D 1 D b 1 , 0 ≥ D 1 , which is equivalent to b 1 ≤ x1 since the function x 1 → f 1 (S in -D 1 x 1 /D, 0) is de- creasing and
- D 1 D 2 < F 1 (x 1 ) = - ∂ f1 ∂x1 ∂ f1 ∂x2 = D1 D ∂f1 ∂S -D2 D ∂f1 ∂S + ∂f1 ∂x2 < 0. ( 23
)
If x 2 = 0, then the equation F 1 (x 1 ) = 0 is equivalent to f 1 (x 1 , 0) = 0 which has a unique solution x 1 = x1 using Proposition 2.2 when condition (7) holds. Similarly, if x 1 = 0, then the equation F 1 (0) = x 2 is equivalent to f 1 (0, x 2 ) = 0 which has a unique solution x 2 = x2 ∈]0, DS in /D 2 [ when condition [START_REF] Fekih-Salem | Sur un modèle de compétition et de coexistence dans le chémostat[END_REF] holds, since the function

x 2 → f 1 (0, x 2 ) is decreasing from f 1 (S in , 0) -D 1 for x 2 = 0 to -D 1 for x 2 = DS in /D 2 .
Proof of Lemma 2.2. Under hypothesis (H0), the function

x 1 → f 2 (x 1 , 0) is contin- uous on [0, DS in /D 1 ] with f 2 (0, 0) = f 2 (DS in /D 1 , 0) = -D 2 .
Then, the equation f 2 (x 1 , 0) = 0 has a solution in [0, DS in /D 1 ] if, and only if,

max x1∈[0,DSin/D1] f 2 (x 1 , 0) ≥ 0,
that is, condition (14) holds.

Proof of Lemma 2.3. Let l 1 be a fixed line defined by x 1 = b 1 that intersects the line δ defined by

D 1 D x 1 + D 2 D x 2 = S in , at point x 2 = c 2 = D D 2 S in - D 1 D 2 b 1 , ( see
f 2 S in - D 1 D b 1 , b 1 ≥ D 2 , which is equivalent to b 1 ∈ [x 1 1 , x 2
1 ] in the particular case where the equation f 2 (x 1 , 0) = 0 admits only two solutions x 1 1 and x 2 1 (see Fig. 1(a,c)). Thus, we have shown that for all b

1 ∈ [x 1 1 , x 2 1 ] there exists a unique b 2 ∈ [0, c 2 [ such that f 2 (b 1 , b 2 ) = 0. Hence, we define the function F 2 by b 2 = F 2 (b 1 )
where the graph γ 2 of this function lies in M . Using the implicit function theorem and assumptions (H2)-(H3), we obtain

- D 1 D 2 < F 2 (x 1 ) = - ∂ f2 ∂x1 ∂ f2 ∂x2 = - D 1 D 2 + D D 2 ∂f2 ∂x1 ∂f2 ∂S . (24)
When x 2 = 0, we have F 2 (x 1 ) = 0, that is, f 2 (x 1 , 0) = 0. From Lemma 2.2, this equation has two solutions x 1 1 and x 2 1 in [0, DS in /D 1 ] since condition (14) holds. Proof of Proposition 3.1. At E 0 , the Jacobian matrix J defined by (18) is written as follows:

J 0 =   -D -a 1 -a 2 0 f 1 (S in , 0) -D 1 0 0 0 -D 2   .
The eigenvalues are negative if, and only if, f 1 (S in , 0) < D 1 . The Jacobian matrix at E 1 is given by

J 1 =   -D -a 1 -a 2 E x1 -E x1 -(E + G)x 1 0 0 f 2 (S in -D1 D x1 , x1 ) -D 2  
where E and G are evaluated at E 1 . The eigenvalues are

f 2 S in - D 1 D x1 , x1 -D 2
which is negative if, and only if, x1 < x 1 1 or x1 > x 2 1 , together with the eigenvalues of the upper left square matrix

A = -D -a 1 E x1 -E x1 . Since tr A = -(D + E x1 ) < 0 and det A = DE x1 + a 1 E x1 > 0,
the real part of the eigenvalues of A are negative. Therefore, E 1 is LES if, and only if, Case 1 or Case 3 holds.

Proof of Lemma 3.2. Using definition (19) and the expressions of F 1 (x 1 ) and F 2 (x 1 ) in ( 23) and ( 24), it follows that

F 1 (x 1 ) = D1 D E -D2 D E -G and F 2 (x 1 ) = - D 1 D 2 + D D 2 H F .
Then, we obtain

F 1 (x 1 ) -F 2 (x 1 ) = D 1 F G -D 2 EH -DGH D 2 F D2 D E + G .
Consequently, at E * , we have

c 3 = D 2 F D 2 D E + G (F 1 (x * 1 ) -F 2 (x * 1 ))x * 1 x * 2 .

B. The particular case for growth functions (21)

In this particular case with the specific growth rates (21), we determine the maximal number of solutions of the equation f 2 (x 1 , 0) = 0 in the following lemma :

Straightforward calculation shows that

φ (x 1 ) = m 2 (-DD 1 K 2 + D 2 1 L 2 )x 2 1 -2DD 1 L 2 m 2 (S in + K 2 )x 1 (K 2 D + S in D -D 1 x 1 ) 2 (x 1 + L 2 ) 2 + m 2 D 2 L 2 S in (S in + K 2 ) (K 2 D + S in D -D 1 x 1 ) 2 (x 1 + L 2 ) 2 .
Hence, the solutions of φ (x 1 ) = 0 are given by those of the following equation

m 2 (-DD 1 K 2 + D 2 1 L 2 )x 2 1 -2DD 1 L 2 m 2 (S in + K 2 )x 1 + m 2 D 2

C. Parameters used in numerical simulations

D. References

Figure 1 .

 1 Figure 1. With the specific growth functions (21): (a) Number of solutions of equation f2(x1, 0) = 0, (b) definition of F1, (c) definition of F2.

Figure 2 .

 2 Figure 2. (a) Change of sign of c4(Sin). (b) Variation of a pair of complex-conjugate eigenvalues as Sin increases. (c) Stable limit cycle when the oscillations are sustained.

  x1 is a solution of[START_REF] Fekih-Salem | Emergence of coexistence and limit cycles in the chemostat model with flocculation for a general class of functional responses[END_REF]. Thus, we have shown that for all b1 ∈ [0, x1 ] there exists a unique b 2 ∈ [0, c 2 [ such that f 1 (b 1 , b 2 ) = 0, that is, each line l 1 meets the set f 1 (x 1 , x 2 ) = 0 exactly once if 0 b 1 x1. Thus, we define the function F 1 by b 2 = F 1 (b 1 ). The graph γ 1 of this function lies in M (see Fig.1(b)). Using the implicit function theorem and from (H2)-(H3), the function F 1 is smooth and decreasing since

Fig. 1 (

 1 c)). From hypotheses (H1)-(H2), the function x 2 → f 2 (b 1 , x 2 ) is decreasing from f 2 (S in -D 1 b 1 /D, b 1 ) -D 2 for x 2 = 0 to -D 2 for x 2 = c 2 . Hence, there exists a unique b 2 ∈ [0, c 2 [ such that f 2 (b 1 , b 2 ) = 0 if, and only if,

Lemma B. 1

 1 Consider the specific growth functions (21) and assume that (14) holds. Then, the equation f 2 (x 1 , 0) = 0 has exactly two solutions x 1 1 and x 2 1 on [0, DS in /D 1 ]. Proof of Lemma B.1. Assume that the growth functions are given by (21). Let φ be the function defined by φ(x 1 ) := f 2 S in -D 1 D x 1 , x 1 .

L 2 S

 2 in (S in + K 2 ) = 0 which has at most two solutions on ]0, DS in /D 1 [ since it is an algebraic equation of degree two in x 1 . Moreover, under (H0)-(H4), φ is positive, continuous and φ(0) = φ(DS in /D 1 ) = 0. Consequently, φ must has a unique positive extremum on ]0, DS in /D 1 [. Therefore, when (14) holds, the equation f 2 (x 1 , 0) = 0 has exactly two solutions x 1 1 and x 2 1 on [0, DS in /D 1 ] (see Fig. 1(a)).

Table 1 .

 1 Parameter values used for (1) when the growth rates f1 and f2 are given by (21).

	Parameter	m1 K1 Ki m2 K2 L2	a1	a2	D	Sin
	Figure 1	4	2	3	8	1.8	1	0.3 0.2	1	3
	Figure 2(a,b)	4	2	3	8	0.1 0.2 0.3 0.2	1	variable
	Figure 2(c)	4	2	3	8	0.1 0.2 0.3 0.2	1	6.27

Proceedings of CARI 2020

Effect of the mortality on a density-dependent model with a predator-prey relationship

Acknowledgments. The first author thanks the Tunisian Ministry of Higher Education, Scientific Research and Technology for financial support and INRAE for hosting her in Montpellier during the preparation of this work. This work was supported by the Euro-Mediterranean research network TREASURE (http://www.inra.fr/treasure).

A. Proofs

Proof of Proposition 2.1. If S = 0 then Ṡ = DS in > 0. Thus no trajectory can leave the positive octant R 3 + by crossing the boundary face S = 0. In addition, if x 1 = 0 and x 2 = 0 then ẋ1 = 0 and ẋ2 = 0 and the set

is invariant under the system (1) because the function t -→ (S(t), x 1 (t), x 2 (t)) = S in + (S(0) -S in )e -Dt , 0, 0


Franklin Platini Agouanet
email: agouanetf@yahoo.com

Jean Jules Tewa B -Magloire

Remy Etoua

Bruce Watson

Eric Badouel

Eds. Oumar Niang

Mathematical modelling and optimal control of banana black sigatoka disease

Keywords: Mathematical modelling, black leaf streak disease, optimal control, sexual reproduction, asexual reproduction, delay MOTS-CLÉS : Modélisation mathématiques, Maladie des raies noires, contrôle optimal, reproduction sexuée, reproduction asexuée, retard

In this work, we propose a mathematical pathogen-host model which focused on the epidemiological cycle of banana black leaf streak disease and which takes into account control strategies to reduce the spores production ,this mainly thanks to an differential equations system with time delay in wich, the host dynamics is coupled to that of the fungus responsible for the disease. By using optimal control theory, we show that an optimal control exists for this problem and Pontryagin's maximum principle is used to characterize an optimal control. Numerical simulations are provided to illustrate our results.

RÉSUMÉ.

Dans ce travail, nous proposons un modèle mathématiques hôte-pathogène axé sur le cycle épidémiologique de la maladie des raies noires du bananier et prenant en compte les stratégies de contrôle pour réduire la production des spores, ceci à partir d'un système d'équations différentielles avec retard dans lequel la dynamique de la plante hôte est couplé avec celle du champignon responsable de la maladie. En utilisant la théorie du contrôle optimal, nous montrons qu'un contrôle optimal existe pour ce problème et le principe du maximum de Pontryagin est utilisé pour caractériser ce contrôle optimal. Des simulations numériques sont faites pour illustrer nos résultats.

Introduction

Banana is grown in more than 120 countries around the world, primarily in the tropical, rain-fed regions of Africa, Asia and Latin America, where there is very important in the food security of more than 400 million people and constitute a employment and income source for local populations [START_REF] Ludivine Lassois | La banane : de son origine à sa commercialisation[END_REF]. However, both banana and banana-plantain cultures are hampered by plant parasitic nematodes, insect pests and foliar fungi(IRAD). Black sigatoka disease also called black leaf streak disease (BLSD) is caused by a pathogen fungi , Mycosphaeralla fijiensis, which is the most costly and damaging banana leaf disease worldwide([1]). BLSD causes significant decrease of photosynthetic surface by a general drying of the leaf system resulting the loss production of ranging from 20 to more than 50% [START_REF] Didy | Caractérisation des populations de Mycosphaerella fijiensis et épidémiologie de la cercosporiose noire du bananier dans la région de Kisangani, RDC[END_REF]. Mycosphaerella fijiensis Morelet which reproduces sexually and asexually. Ascospores are the product of sexual reproduction, while conidia are produced through asexual reproduction and these spores are responsible for the spread of the disease. This fungus performs its entire life cycle on the banana tree and their disease cycle consists of four distinct stages that include spore germination, penetration of the host, symptom development and spore production.([1]) The susceptibility of the crop necessitates the cultural practices and use of multiple fungicides at relatively high frequencies. Such applications are potentially detrimental, not only to the environment, but also to those who live and work in areas in which banana is treated to control the disease([1]). In this paper, we propose a mathematical pathogenhost model which focused on the epidemiological cycle of the disease and which takes into account control strategies to reduce the spores production.

The model formulation

We formulate an optimal control pathogen host model for black sigatoka disease in order to derive optimal control strategies with minimal implementation cost. The control function, u represent time dependent efforts made to reduce the ascospores(sexual spore) and conidia(asexual spore) reproduction practiced on a time interval [0, T ] where T is the cropping season duration. Known practices of control efforts include cultural practices and treating by fungicides. We divide the host population(N) suppose constant into two compartments using variables S and I to describe the total number of susceptible leaf and infected leaf, respectively. On the other hand, Y and Z describe respectively the sexual and asexual spore groups in the pathogen population. Note that the term 'leaf' is a shorthand for the "leaf part that a lesion occupies", so that multiple infections cannot occur in the model. In the Moungo region (Cameroon), the incubation times are the shorter (13-16 days) were obtained between May and November, a period more favorable for illness, while the longest was obtained in February (26 days)( [START_REF] Fourée | Contibution à l'étude épidémiologique de la cercosporiose noire dans la zone bananiaire du Mungo de 1987 à 1989[END_REF]), so the delay (τ ) represents the incubation period of the disease. We assumed that, all the leaves are destroyed and the fungus also disappears at the end of cropping season, and that the spore dispersal rate is equal to the spore deposite rate [START_REF] Frédéric | Mate Finding, Sexual Spore Production, and the Spread of Fungal Plant Parasites[END_REF].

let β and α the rate of deposit respectively sexual and asexual spore upon leaves, and 0 < p < 1 be their infectivity(we assume that conidia infectivity is equal to the ascospores infectivity [START_REF] Landry | Modélisation des dynamiques de maladies foliaires de cultures pérennes tropicales à différentes échelles spatiales: cas de la cercosporiose noire du bananier[END_REF]), i.e. the probability that a spore in contact with a susceptible leaf succeeds to infect it. Let σ be the number of asexual spores produced per infected leaf per unit time, and γ be the number of sexual spores (ascospores) produced per infected leaf per unit time.

-Infestation function is given by p(βX + αY ) S N . Hence, use the fact that the density of susceptible leaves can be deduced from the infected leaf density through the delayed equation s(t -τ ) = N -x(t -τ ), where 0 < τ < T stands for the time delay, we have the dynamics of infected leaves given by:

İ(t) = p(βX(t) + αY (t)) 1 - I(t -τ ) N -mI(t)
where m is the infectious leaves mortality rate.

-The asexual spore dynamics is given by Ẏ (t) = σI(t) -βY (t) where σI(t) is the asexual spore production function.

-Sexual spore production thus results from the interaction between two individuals with compatible mating types (+ and -). However, sexual spore production is conditioned to the local presence of a mating partner, whose probability is (I/2)/N , assuming a balanced mating-type ratio and the sexual sopre production function is given by γ 2N I 2

[5]; hence the sexual spores dynamics is given by

Ẋ(t) = γ 2N I(t) 2 -βX(t).
Hence, our model is given for t ∈ [0, T ] by:

               Ẋ = p(βY + αZ) 1 - X(t -τ ) N -mX Ẋ = (1 -u) γ 2N X 2 -βY Ż = (1 -u)σX -αZ
Let us use the following variable change:

x = I N , y = Y N and z = Z N , so for t ∈ [0, T ] our model is gigen by:      ẋ = p(βy + αz) [1 -x(t -τ )] -mx ẏ = (1 -u) γ 2 x 2 -βy ż = (1 -u)σx -αz (1) with initial conditions y(θ) = y 0 ≥ 0 , z(θ) = z 0 ≥ 0 et x(θ) = ϕ(θ) ≥ 0 for θ ∈ [-τ, 0] where ϕ ∈ C([-τ, 0], R), the Banach space of continuous functions mapping the interval [-τ, 0] into R + equipped with the norm ||ϕ|| = sup -τ ≤t≤0 |ϕ(t)|.
We given the parameters values in the following table.

Our problem consists in maximizing the yield at the end of the cropping season, while controled the champignon propagation and minimizing the control cost. We propose the following objective function:

J(u) = T 0 (u(t) 2 + Bx(t))dt.
where B is the weight constants of control of the infected group. The set of admissible controls is defined as follows

Γ = {u ∈ L 1 ([0, T ])/0 ≤ a ≤ u ≤ b ≤ 1; ∀t ∈ [0, T ]}

Preliminary result

In this section, we establish the existence of a solution to the delay differential system using results from R.D. Driver's text [START_REF] Driver | Ordinary ans delay Differential Equations[END_REF]. For notational purposes, we use :

f (t, X) =    p(βy + αz)(1 -x(t -τ )) -mx (1 -u) γ 2 x 2 -βy (1 -u)σx -αz   
where X = t (x, y, z). we can define a new function

U t : [-τ, 0] -→ R θ -→ U t (θ) = U (t + θ) where U ∈ {x, y, z}.
So the cauchy problem associed on systm (1) is given by

Ẋ(t) = F (X t ) X 0 = ϕ where F (x t (0), y t (0), z t (0), x t (-τ )) = f (x(t), y(t), z(t), x(t -τ )).
By [START_REF] Driver | Ordinary ans delay Differential Equations[END_REF], Using the fact that the right hand side of (1) has continuous partial derivatives, we obtain a local Lipschitz condition of F and we conclude for the existence of solution. Now, we show that all these variables are always positive when time evolves.

Theorem 3.1. For the initial conditions (x(θ), y(θ), z(θ)) ∈ R 3 + , the solution (x(t), y(t), z(t)) of system (1) are non negative and bounded for all time t > 0.

Proof. see Appendix 1.

For u = 0; let R 0 = pσ m and S 0 = p( γ 2 + σ) 2 2γm
. Let us make the following realist assumption:

γ 2 > σ Theorem 3.2.
1) if R 0 < 1 and S 0 < 1, the system (1) has only disease free equilibrium E 0 .

2) if R 0 < 1 and S 0 > 1, the system (1) has three equlibria point: E 0 and two endemic equilibria E + and E -.

3) if R 0 < 1 and S 0 = 1, the system (1) has two equlibria point: E 0 and one endemic equilibria E * .

4) if R 0 ≥ 1, the system (1) has two equilibria: DFE E 0 and unique endemic equilibrium E.

Proof. see Appendix 2.

Existence of an optimal control

We note that the existence of an optimal control can be proved by using results from Fleming and Rishel, [(theorem 4.1 pp 68-69 [START_REF] Fleming | Deterministic and Stochastic optimal control[END_REF]), for this, we will check the following properties..

1) The set of controls and corresponding state variables is non-empty.

2) The control set Γ is convex and closed.

3) The right hand side of the state system is bounded by a linear function in the state and control.

4) The integrand of the objective functional is convex on Γ and there exist nonnegatif constants c 1 , c 2 and β > 1 satisfying L(u(t), X(t)) ≥ c 1 |u| β/2 -c 2 .

Theorem 4.1. There exists an optimal control u * and a corresponding solution (x * , y * , z * ) of the initial value problem (1) that minimizes the cost function J in Γ such that J(u * ) = min a≤u≤b J(u).

Proof. An existence result in Lukes ([ 10], Theorem 9.2.1) for the state system (1) with bounded coefficients on the finite interval time is used to give condition 1). By definition, the control set Γ is convex and closed, so condition 2) is satisfied. The right hand side of system (1) is continuous and can be written as a linear function of u with coefficients depending on time and state. Furthermore, by Γ and theorem 3.1, all variables x, y, z and u are bounded on [0, T ], so condition 3 hold. For 4., the integrand function L of the objective functional is clearly convex in the controls since it is quadratic. Then, for c 1 = 1, and β = 4, there exist c 2 such that

L(t, X, u) = u 2 + Bx ≥ |u| β/2 ≥ c 1 |u| β/2 -c 2 .

Characterization of the optimal control

With the existence of the optimal control established, we now present the optimality system using the Pontryagin's minimum principle for delayed control problem from Göllman and al [START_REF] Göllmann | Optimal ccontrol problems with delays in state and control variables subject to mixed control state constraints[END_REF]. We discuss the theorem that relates to the characterization of the optimal control. The optimality system can be used to compute candidates for optimal control. To do this, we begin by defining a augmented Hamiltonian.

For the adjoint variables Λ = (λ 1 , λ 2 , λ 3 ) and penalty multipliers, w i (t), i=1,2 ; that are attached to the control constraints; We define the augmented Hamiltonian as follows:

H(u(t), X(t), x(t -τ ), Λ(t)) = u(t) 2 + Bx(t) + λ 1 (t) [p(βy(t) + αz(t))(1 -x(t -τ )) -mx(t)] +λ 2 (t) (1 -u(t)) γ 2 x(t) 2 -βy(t) + λ 3 (t) [(1 -u(t))σx(t) -αz(t)] +w 1 (a -u(t)) + w 2 (u(t) -b),
where w 1 et w 2 are non-negatif and satifiyng w 1 (a -u * (t)) = w 2 (u * (t) -b) = 0 for optimal control u * .

Theorem 5.1. Given optimal control u * and the corresponding solutions x * , y * , z * of system (1), there exist adjoint variables λ 1 , λ 2 and λ 3 verifiyng the system:

         λ1 (t) = -B + λ 1 m -λ 2 (t)(1 -u * ) γ 2 x * (t) -λ 3 (t)(1 -u * )σ +X [t0,T -τ ] (t)λ 1 (t + τ )p(βy * (t + τ ) + αz * (t + τ )) λ2 (t) = -λ 1 (t)βp(1 -x * (t -τ )) + βλ 2 (t) λ3 (t) = -λ 1 (t)αp(1 -x * (t -τ )) + αλ 3 (t) (2)
with transversality conditions λ i (T ) = 0; i = 1, 2, 3.

Furthermore, u * are represented by u * = max{a, min{b, 1 4

x * (γx * λ 2 + 2σλ 3 )}} Proof. Using the Pontryagin's minimum principle for delayed control problem from Göllman and al [START_REF] Göllmann | Optimal ccontrol problems with delays in state and control variables subject to mixed control state constraints[END_REF], adjoint differential equation is giving by

                   λ1 (t) = - ∂H ∂x(t) (t) -X [t0,T -τ ] (t) ∂H ∂x(t -τ ) (u(t + τ ), X(t + τ ), x(t), Λ(t + τ )) λ2 (t) = - ∂H ∂y(t) (u(t), X(t), x(t -τ ), Λ(t)) λ3 (t) = - ∂H ∂z(t) (u(t), X(t), x(t -τ ), Λ(t))
with traversality conditions λ i (T ) = 0; i = 1, 2, 3. where X [t0,T -τ ] is the "indicatrice" function on [t 0 , T -τ ], hence we obtain system(2).

Also, the optimal control u * can be solved from the optimality condition

∂H ∂u = 0, that give u * = 1 4 x * (γx * λ 2 + 2σλ 3 ) + w 1 -w 2 .
To determine an explicit expression for the optimal control without w 1 and w 2 a standard optimality technique is utilized. We consider the following three cases. i) On the set {t : a < u * (t) < b}, we have w 1 = w 2 = 0 and u * (t) = 1 4

x * (t)(γx * (t)λ 2 (t) + 2σλ 3 (t)).

ii) On the set {t : u * (t) = a}, we have w

2 (t) = 0. Hence a = u * = 1 4 x * (γx * λ 2 + 2σλ 3 ) + w 1 . Therefore 1 4 x * (γx * λ 2 + 2σλ 3 ) ≤ a since w 1 (t) ≥ 0. iii) On the set {t : u * (t) = b}, we have w 1 (t) = 0. Hence b = u * = 1 4 x * (γx * λ 2 + 2σλ 3 ) -w 2 . Therefore 1 4 x * (γx * λ 2 + 2σλ 3 ) ≥ b because w 2 ≥ 0.
Thus i), ii) and iii) can resume at

u * (t) = max{a, min{b, 1 4 x * (t)(γx * (t)λ 2 (t) + 2σλ 3 (t))}}

Numerical simulations

In this section, we present the numerical solution of our control problem and compare it with the solution in the absence of control. We adopted the forward-backward sweep method to solve numerically our optimal model with time delay [START_REF] Hattaf | Optimal Control of Treatment in a Basic Virus Infection model[END_REF]. After making an initial guess for the control functions, we first solved the initial valued state system forward in time. Then, using the same initial guess for the control functions, we solve the adjoint system with terminal conditions backward in time. The controls are updated in each iteration using the optimality conditions. The parameters used in solving the optimality system are estimates of data from BLSD as summarized in Table 1 with m = 0.08, τ = 15, a = 0 and b = 1. With this values, the threshold give R 0 = 2.5. The cropping season of banana variyng to 10 at 12 months [START_REF] Banana | Banana Planters[END_REF], so we take T = 365 days.

The simulations plots are given in figure1 for B = 2 (A, B, C, D) and figure(2) for B = 5 (a, b, c, d) where We compare the cases with control (red curves) and without control (blue curves). We observe that the control u(t) ( see (D) and d) ) reduces the disease but its effectiveness depends of parameter B which can represent the fungicides concentration on the leaf in the case of chemical control.

Figure 1: Simulation for B = 2. We remark that the control is less effective because the weight of control B on the infected leaves is law, nevertheless we can see that the control u(t) (see D)) grows weakly up to 160 days where it reaches its maximum and remains there up to 225 days (which corresponds to the flowering time [START_REF] Banana | Banana Planters[END_REF]) and begins to decrease.

Figure 2: Simulation of system (1) with B = 5. We can remark that the control is more effective because the weight of infected leaves is high. We can see that the control u(t) (see d)) grows rapidly and reaches its maximum from 50 days and remains so until 230 days( which corresponds to the flowering time [START_REF] Banana | Banana Planters[END_REF]) before starting to decrease slowly.

Conclusion

In this work, we have proposed a mathematical pathogen-host model with time delay of banana black leaf streak disease which takes into account control strategies to reduce the spores production. We have designed an optimal control problem that consists to maximizing the yield at the end of the cropping season, while controled the champignon propagation and minimizing the control cost we have defined. We have characterized the optimal control using the Pontryagin's maximum principle and we have solved numerically our system. For this numerical simulations, we saw that the control reduces the disease and its effectiveness depends of parameter B of our cost function which can represent the constant weight of control on the infected leaf.

[1] A.C.L. CHURCHILL, " Mycosphaerella fijiensis, the black leaf streak pathogen of banana:

progress towards understanding pathogen biology and detection, disease development, and the By equation(2) of system(1), we have ẏ(t) ≤ γx * 2 -βy(t) which implies that y(t) ≤ x * γ 2β

and the same by eqyation(3) of system(1), we deduce that z(t) ≤ σx * α , therefore x, y and z are bounded.

Appendix 2.Proof of theorem(3.2)

For u = 0 and p = q, the equilibria of system(1) are the solutions of system:

     p(βy + αz)(1 -x) -mx = 0
(1) γ 2

x 2 -βy = 0 (2)

σx -αz = 0 (3)
. From ( 2) and ( 3

), whe have y = γ 2 x 2 = 0 or pγ 2 x 2 -p( γ 2 -x 2 ±= γ 2 x * 2

 222222 and z = σ α x. Introduce in (1), we obtain x σ)x + m -pσ = 0. (E) For x = 0, we have y = z = 0 and whe obtain disease free equilibrium E 0 = (0, 0, 0). If m -pσ ≤ 0 or R 0 = pσ m ≥ 1, then using the Descarte sign ruler's, equation (E) has a unique positive solution x; then y = γ 2β x 2 , and z = σ α x, hence endemic equilibriumgiven for R 0 ≥ 1. If R 0 < 1, then the discriminant of equation (E) is given by ∆ = p 2 ( γ 2 +σ) 2 -2pγm.Hence, if S 0 > 1, i.e. ∆ > 0 then (E) has two positive solutions x + and x -and we obtain equilibria E + and E -with y ± = γ 2 and z ± = σ α x ± . Finally if S 0 = 1 that implies ∆ = 0, (E) has a unique positive solution x * and we have equilibrium E * with y * and z * = σ α x * .

Table 1 :

 1 Dimensional variables and their estimated values for M. fijiensisThe problem now is to find the control u * satisfying J(u

	Parameter	Meaning	value	Reference
	β	Sexual spores deposition rate	48 per day F Hamelin and al, 2017 [5]
	α	Asexual spores deposition rate	2	assumed
	p	Sexual spores infectivity	0.01	Landry , 2015 [9]
	γ	Sexual spore production rate	160 per day Virginie and al, 2017 [12]
	σ	Asexual spore production rate	20 per day	Virginie and al, 2017 [12]
	m	infectious mortality rate	0.04 to 0.1	assumed

* ) = min a≤u≤b J(u)

Proceedings of CARI 2020

Mathematical modelling and optimal control of banana blak sigatoka disease

Appendix 1.Proof of theorem(3.1)

The model system can be put into the matrix form X = g(X) where X = (x, y, z) T and

We have g 1 (0, z, y) ≥ 0 , g 2 (x, 0, z) ≥ 0 and g 3 (x, y, 0) ≥ 0. Hence, every solution which starts in R 3 + will remain in R 3 + . We prove now that the solutions of model system(1) are bounded. using the fact that x is continous over the finite interval [0, T ], we conclude that x has a maximum x * in [0, T ].


Nedjma Djezzar
email: n.djezzar@univ-batna2.dz

Noureddine Djedi

Yves Duthen

Open Archive Toulouse Archive Ouverte

All living organisms depend on highly complex networks of intercellular signals to coordinate cell-cell communication in various physiological functions. A large number of mechanisms that rely only on cell signaling are capable of development, collective behavior, self-repair and immunity. This abstract proposes the use of an ANN-based GRN model to orchestrate cell-cell communication in two contexts: unicellular organism via Quorum Sensing and multicellular organisms via Notch signaling pathway. A simple grammar inspired from membrane computing is proposed to describe the dynamic cell-like structure.

Introduction

Communication is vital to life. It represents the coordination of cells and tissues to accomplish various physiological functions as well as the adaptation to the environment. All forms of life on earth (eukaryotes, prokaryotes and metazoan organisms) need to communicate to survive.

Over hundreds of millions of years, evolution of metazoan organisms has produced a complexity both in phenotypes as well as the developmental mechanisms by which this complexity is generated. It is now known, in developmental biology, that this process is based on good coordination of intra and intercellular communication signals. Indeed, a small toolkit of cell signaling co-opted recurrently during development is sufficient for tissue morphogenesis, angiogenesis and metastasis.

Cell-cell communication occurs when a molecular signal is emitted by a cell called emitting cell and recognized by another cell called receiving cell. Basically, signals are composed of growth factors, hormones, cytokines, and neurotransmitters that bind to specific cell receptors. This ligand-receptor interaction then generates several types of intracellular signals that may be relayed by downstream signaling and by protein phosphorylation cascades. Finally, these intracellular signals lead to the activation of transcription factors that regulate target genes responsible for the expression of diverse biological outcomes: activation of cell cycle, proliferation, differentiation.

It does not come as a surprise that "Unicellular organisms including both eukaryotes (e.g. yeast) and prokaryotes (e.g. bacteria) utilize similar communication systems as multicellular organisms. For example: receptors for insulin are present in Drosophila which are similar to their vertebrate counterparts in binding properties, general structure, and insulin-stimulated enzyme activity." [START_REF] Leroith | Receptors for intercellular messenger molecules in microbes: similarities to vertebrate receptors and possible implications for diseases in man[END_REF] In the artificial life community, cellular interactions have been investigated using cellular automata [START_REF] Gerlee | Evolving homeostatic tissue using genetic algorithms[END_REF][START_REF] Nichele | Genotype Regulation by Self-modifying Instruction-Based Development on Cellular Automata[END_REF], agent-based models [START_REF] Ferreira | Dynamic Structure Discovery and Repair for 3D Cell Assemblages[END_REF], and GRN-based dynamics [START_REF] Kim | The relationship between microscopic and collective properties of gene regulatory network-based morphogenetic systems[END_REF]. In most of these models, cell-cell interactions are abstract and they do not rely on an advanced exploration of the underlying biological mechanisms. In fact, communications between cells are based on: a diffusion system, information pathway passed through neighboring cells, or the concentration of a molecule that has a direct effect in regulating gene expression and cell behavior. Moreover, communication between unicellular organisms like bacteria Quorum Sensing has not been well investigated by the Artificial Life community. Here, we assume that more realistic unicellular and multicellular approaches must involve more investigation of, the architecture of the chemical communication network, how chemical information is transduced to control gene expression; how intra-and interspecies cell-cell communication is accomplished; and the fascinating option of unicellular-multicellular crosstalk.

ANN Based GRN Model

Various modeling techniques have been developed to simulate the functioning of the GRN. They can be divided into two classes. The first class: logical models, including Boolean networks and their variants. Although effective in dealing with noise, they fail to consider temporal dynamic aspects that are important features of regulatory network [START_REF] Vijesh | Modeling of gene regulatory networks: A review[END_REF]. The second class: Continuous models based essentially on differential equations allow more detailed descriptions of GRN dynamics.

To describe the dynamics of a generic GRN supervising cellcell communication, we use an Artificial Neural Network (ANN) modeling technique introduced by Vohradsky [START_REF] Vohradsky | Neural network model of gene expression[END_REF]. Formally, an ANN is defined by: We extend the formula to take into account cellular processes (e.g., signal translation) that influence the gene sensitivity.

𝑑𝑋 𝑖 𝑑𝑡 ⁄ = 𝑣 𝑖 𝑓 (∑ 𝑤 𝑖𝑗 𝑋 𝑗 -𝜗 𝑖 𝑁 𝑗=1 ) -𝑘 𝑖 𝑋 𝑖 𝜗 𝑖 𝑣 𝑖 𝑘 𝑖 > 0  𝑋 𝑖 , with 𝑖 = 1, … ,

Dynamic cell-like structure

Based on the idea that a cell is a complex structure, we consider a cell as a multi-set of objects arranged in regions and delimited by another main object which is the membrane. Each object within the cell is composed of a multi-set of objects also. An object can be: a lipid, a protein, a receptor, a cytoplasm or a nucleus. So, the model of the cell-like structure is a hierarchical arrangement of multi-sets of objects.

To simulate this hierarchical arrangement of the cell structure, we introduce a simple cellular grammar inspired from formal grammars like in Membrane Computing [START_REF] Paun | Membrane computing: an introduction[END_REF]. However, rules of the grammar don't govern (in any way) development. Development occurs through cellular interactions that lead to gene expression and regulation. Thus, rules are just a formal description of the cell-like structure. This formal description can be used to observe the resulting changes in the cell structure and deduce the resulting grammar. Membrane computing is a computational model that solves NP complete problems. This can be used to investigate the computational power of the resulting digital organisms.

Application on unicellular organisms (Quorum Sensing)

Bacterial cell-cell communication, termed Quorum Sensing (QS), is a signal-response process that regulates gene expression according to the cell density in order to coordinate collective behaviors. In [START_REF] Ouannes | Emergent group behaviors from bacteria quorum sensing simulation[END_REF], we have introduced a molecular network that uses Ordinary Differential Equations (ODEs) to calculate the concentration of Gram negative bacteria QS components. The simulation shows that bacterial cells are able to form micro-colonies and behave like a multicellular organism (figure 1). The principal limitation of the model resides in the use of a proposed algorithm of QS to simulate bacteria cells communication instead of a regulation mechanism that describes the mutual enhancing-inhibiting influences in genes of a GRN. Moreover, bacterial response to QS is various and includes biofilm formation, bioluminescence, drug-resistance... Our work focuses on investigating the underlying mechanisms behind these phenomena; how chemical information is transduced to control gene expression, and observing in simulation the emerging phenotypic traits.

Figure 1: Micro-colonies of moving bacteria emerged from QS communication [START_REF] Ouannes | Emergent group behaviors from bacteria quorum sensing simulation[END_REF].

Application on multicellular Organisms (Notch signaling pathway)

Notch signaling plays a central role in cell differentiation decision via the lateral inhibition process, which is a common phenomenon during selection of neural and muscle precursor cells [START_REF] Bray | Notch signalling in context[END_REF]. Briefly, two or more equivalent cells express both the receptor Notch and the ligand. A negative feedback loop operates so that ligand expression became repressed in some cells which will be specialized on receiving cells, and this will lead to asymmetry. Despite the fact that the Notch pathway has a myriad of roles in different developmental and disease contexts, it is relatively simple in its functioning. This is why, we made the choice of using specifically this signaling pathway as instance of our model.

Conclusion

In the lineage of our previous work [START_REF] Djezzar | L-systems and artificial chemistry to develop digital organisms[END_REF], we use Generative Developmental Systems (GDS) that consist of a grammar inspired by membrane computing to simulate the dynamic structure of the cell. This is made to explore the computational power of the emerged organisms. Multicellular organisms are ecosystems that cohabit cells of the organism and other microorganisms that are bacteria. We assume that a link must be made for a future study of both the artificial ecosystem and the pathological context (QS virulence pathway) considering the two species.

  𝑁, is the i-th gene product concentration.  𝑣 𝑖 : is the maximal expression rate of gene 𝑖.  𝑤 𝑖𝑗 is the connection weight of control of gene j on gene i. Positive values of wij indicate enhancing influences while negative values define inhibiting influences.  𝜗 𝑖 : is the influence of external input on gene i, which modulates the gene's sensitivity of response to activating or repressing influences. 𝜗 𝑖 to keep the activation from growing without bounds.  𝑘 𝑖 Degradation of the i-th gene expression product.

 𝑓: Represents a non-linear sigmoid transfer function modifying the influence of gene expression products 𝑋 𝑖 and external input