MODELLING COFFEE LEAF RUST DYNAMICS TO CONTROL ITS SPREAD

Coffee leaf rust (CLR) is one of the main diseases that affect coffee plantations worldwide. It is caused by the fungus Hemileia vastatrix. Damages induce severe yield losses (up to 70%). Its control mainly relies on cultural practices and fungicides, the latter having harmful ecological impact and important cost. Our goal is to understand the propagation of this fungus in order to propose a biocontrol solution, based on a mycoparasite that inhibits H. vastatrix reproduction. We develop and explore a spatio-temporal model that describes CLR propagation in a coffee plantation during the rainy and dry seasons. We show the existence of a solution and prove that there exists two threshold parameters, the dry and rainy basic reproduction numbers, that determine the stability of the equilibria for the dry and rainy season subsystems. To illustrate these theoretical results, numerical simulations are performed, using a non-standard finite method to integrate the pest model. We also numerically investigate the biocontrol impact. We determine its efficiency threshold in order to ensure CLR eradication. Mathematics Subject Classification. 35K57, 93D05, 65M06, 92D30. Received March 2, 2020. Accepted March 11, 2021.

The life cycle of H. vastatrix in favorable humidity and temperature conditions lasts 5 weeks on average [41]. Urediniospores are dispersed by rain and wind. They germinate and penetrate through stomates on the underside of the leaf. This infection process requires 24 to 48 hours [25]. The first symptom is a pale yellow lesion that appears 1 to 3 weeks after infection. Sporulation, i.e. the production of urediniospores and teliospores, occurs 2 weeks to several months after infection [41]. A single lesion produces 4 to 6 crops of spores, releasing 300,000 to 400,000 urediniospores over a period of 3 to 5 months [25].
As most rust fungi, H. vastatrix produces urediniospores for asexual reproduction and teliospores for sexual reproduction. Teliospores do not infect coffee leaves and have no known host [30], so it is not established that there is an effective sexual reproduction cycle.
During the dry season, H. vastatrix survives primarily as mycelium in the living tissues of the coffee leaves. As infected leaves drop prematurely during the dry months, a large amount of potential inoculum for the next rainy season disappears. However, a few green leaves always persist. Moreover, urediniospores can survive about 6 weeks, so there is always some viable inoculum to infect the newly formed leaves at the start of the next rainy season [2].

Model formulation
Based on this biological knowledge, we propose a spatio-temporal coffee-CLR interaction model within a coffee plantation. We represent space and urediniospore diffusion inside a bounded domain. At each point in space, we describe the evolution of the number of coffee branches, according to their epidemiological state, based on the CLR life cycle: healthy, latent with spores germinating or non sporulating lesions, infectious with sporulating lesions, and leafless at the end of the infection process. The infection of a branch corresponds to the infection of its leaves. We consider two periods: the rainy season during which berries are produced, and the dry season. This implies that most parameters are time-dependent, as they take two values according to the season.
We obtain the following model, that is illustrated in Figure 1: ∂ t S(t, x) = Λ(t) − ω(t)νU (t,x) N (t,x) S(t, x) − µ(t)S(t, x), ∂ t L(t, x) = ω(t)νU (t,x) N (t,x) S(t, x) − (θ(t) + µ(t))L(t, x), ∂ t I(t, x) = θ(t)L(t, x) − (α(t) + µ(t) + d)I(t, x), ∂ t J(t, x) = α(t)I(t, x) − µ(t)J(t, x), ∂ t U (t, x) = ε∆U (t, x) + γ(t)I(t, x) − (ν + µ U )U (t, x), where S(t, x), L(t, x), I(t, x), J(t, x), U (t, x) and B(t, x) represent the densities of healthy branches, latent branches, infectious branches, leafless branches, urediniospores and berries respectively, at time t and location x defined in Ω, a simply connected bounded domain (of R or R 2 ) with smooth boundary ∂Ω. The densities are /m or /m 2 depending on the choice of Ω. The total density of branches is N (t, x) = S(t, x) + L(t, x) + I(t, x) + J(t, x). The recruitment of healthy branches occurs at rate Λ(t). Urediniospores are deposited on leaves of all branches at rate ν and a fraction S/N lands on healthy branches. These spore covered healthy branches become latent branches with rate ω(t), which represents the germination efficacy, that is the number of healthy branches which become latent branches, per deposited spore. The latent branches become infectious at rate θ(t), where 1/θ(t) corresponds to the latency period, and in turn become leafless branches at rate α(t), where 1/α(t) corresponds to the sporulation period. All branches undergo natural mortality with baseline rate µ(t) and the infectious branches have an additional constant mortality rate d due to the disease. Urediniospores are produced by infectious branches at rate γ(t) and lose their ability to infect coffee branches at constant rate µ U . Parameter ε is the urediniospore diffusion coefficient. Berries are produced by all types of branches at different rates: They have a constant mortality rate µ B . All parameters are assumed to be non negative.
Having two seasons, let us consider T the yearly period and τ the duration of the dry season. We set the initial time t = 0 at the beginning of a dry season. Then, for n years, we assume that time-dependent parameters m(t), where m(t) ∈ {Λ(t), ω(t), θ(t), α(t), γ(t), µ(t)}, and δ i (t), with i ∈ {S, L, I, J}, are constant during each season and can be written as follows: where D represents the dry season and R the rainy season. Table 1 summarizes the biological meaning of parameter values for system (2.1).
with B(0, x) = 0 ∀x ∈ Ω. We assume that the harvest occurs instantaneously at the end of each rainy season. Thus B(t, x) is brought back to 0 at t = nT , with other variables unchanged. We suppose that the domain Ω is isolated, which implies that spores do not get in or out and that there is no interaction with the outside environment. To ensure this property, we assume that there is a large enough buffer zone without coffee plants around the coffee plantation (hence around the domain Ω), so that no urediniospores are introduced from outside the plantation. Moreover, we neglect the spores that could potentially get out of the plantation. This absence of transfer is classically represented through the homogeneous Neumann boundary: where η denotes the unit outward normal on ∂Ω.

Mathematical analysis
Herein, we present the basic properties of the subsystems of system (2.1), defined over each season. For any parameter m, we denote by m p the seasonal constant value of the parameter m where p = R for the rainy season and p = D for the dry season. Replacing the value of every parameter m(t) of system (2.1), with the corresponding constant value m p with respect to the season, and removing the ∂ t B equation since B is not present in the other equations of system (2.1), we obtain: (3.1)

Basic properties
In this section, we will first establish positivity of the solutions if they exist in Lemma 3.2, their boundedness in Lemma 3.3, and finally their existence and uniqueness in Lemma 3.6. Through these, we obtain positivity, boundedness and global existence of solutions of subsystem (3.1) in Theorem 3.7.
In order to show positivity, we will be using the following lemma of [26].
Since the initial values are non-negative and the growth functions on the right-hand side of system (3.1) are assumed to be sufficiently smooth in R 5 + , we have the following result over the positivity of a solution of subsystems (3.1), if it exists: Using the definition of supremum, there does not exist t ∈ [0, t 1 ) and x ∈ Ω such that any variable is equal to zero.
Among S, L, I, J, U , let us first consider S to be the variable such that there exists x with S(t 1 , x) = 0, and let us define: One has: where Λ p is obtained from the expression of ∂ t S(t, x) in subsystems (3.1). By integration of this expression over [0, t 1 ], we obtain: Since λ(t 1 , x) > 0 and φ S (x)λ(0, x) > 0, we have S(t 1 , x) > 0, which contradicts the hypothesis that S(t 1 , x) = 0.
Also, the ∂ t L, ∂ t I and ∂ t J equations are linear scalar equations with forcing terms ωpνU N S, θ p L and α p I that are positive over the interval [0, t 1 ); hence L(t 1 , x), I(t 1 , x) and J(t 1 , x) are positive.
Hence, if t 1 < T max as we supposed, we must have, for some x, U (t 1 , x) = 0. As I(t, x) > 0 over [0, t 1 ], the last equation of subsystems (3.1) then gives rise to the following inequality system: × Ω, it follows from Lemma 3.1 that U (t, x) > 0 in [0, t 1 ] × Ω. This contradicts U (t 1 , x) = 0, and hence the assumption t 1 < T max . Therefore all variables are positive for all times smaller than T max .
This achieves the proof. with bounded initial conditions. Then this solution is bounded, T max = +∞, 0 < N (t, x) = S + L + I + J ≤ N m and 0 < U (t, x) ≤ U m , where: The N -dynamics then satisfy: upon which we build the upper-bounding system The comparison principle then yields: x) stays bounded and is hence defined for all times: This upper-bound also holds for S, L, I and J which are integral parts of N .
From the upper-bound on I and the last equation of subsystems (3.1), we can deduce that: and the upper-bounding system: It then follows that for all (t, x) ∈ [0, T max ) × Ω: Hence, as was done for N , we have: This completes the proof Now we prove existence and uniqueness of a local solution of subsystems (5). The idea is to write the system in the abstract form on a suitable Hilbert space and then to use the use the Hille-Yosida Theorem to conclude We therefore introduce the linear space H = L 2 (Ω) 5 endowed with the usual inner product: H is a Hilbert space. Now define on H the linear operator A such that: We will need the following result.

Lemma 3.5 ([28]
). Let ε be a positive real number. Let f ∈ L 2 (Ω) and g the trace on ∂Ω of an element of H 1 (Ω). Assume that Ω is a bounded subset of R N , N = 1, 2 with smooth boundary and consider the following stationary boundary value problem: Then, the problem (3.9) has a unique solution u ∈ H 1 (Ω). Moreover u belongs to H 2 (Ω) and there exists a universal constant M > 0 that only depends on ε, N and Ω such that We have the following result about the existence of solution of subsystems (3.1).
Lemma 3.6. Assume that the initial conditions (2.3) holds, then there exists a unique local solution of problem (3.1) defined on [0, T max ) × Ω. More precisely: 3) can be written abstractly in the Hilbert space H in the following form: Ay is defined as in equation (3.7) and: Note that, since Ω is bounded, for y = (S, L, I, J, U ) ∈ H, one has F (y) ∈ H. Now let us show that the linear operator A is a maximal monotone operator on H.
With this is mind, one has: Therefore, for i = 1, . . . , 4, u i = 1 1+µ P v i ∈ L 2 (Ω) and one obtains u 5 by solving the following elliptic problem: According to Lemma 3.5, there exists a unique u 5 ∈ H 2 η (Ω) solution of (3.11) and A is maximal. By the Hille-Yosida theorem [28], we conclude that −A is the infinitesimal generator of a C 0 semigroup of contractions on H.
-We now show that F is Lipschitz continuous in both variables. Let y = (S, L, I, J, U ),ỹ = (S,L,Ĩ,J,Ũ ) ∈ H. Therefore, one has: , we obtain: where K is a constant that depends on the constant ω P , ν, θ P , α P , d, γ P and µ U . Thus F is uniformly Lipschitz continuous on H. Using the fact that the solution are positive, we can now apply Theorem 1.6 page 189 of [28] and conclude that system (3.10) has a unique local strong solution on [0, T max ) in the sense that (S, L, ; H), we can also use the method present in [43] to conclude prove. This completes the proof.

Equilibria and their stability
Here, we compute the equilibria of subsystems (3.1) with p = D or R and study their stability. Let The expression R (p) 0 is derived in (A.6) in the proof of the local stability of the disease free equilibrium. It corresponds to the basic reproduction number in the dry (p = D) or rainy (p = R) season, since (i) γp (ν+µ U ) represents the mean number of new urediniospores generated by a single infectious branch and (ii) measures the average number of new infectious branches generated by a single urediniospore introduced in a completely susceptible field. We prove the following results for the stability of equilibria of subsystems (3.1).
0 > 1 and close to one, the DFE is unstable and there exists a unique positive endemic equilibrium which is locally asymptotically stable.
Proof. The proof is given in Appendix.
Using Lyapunov theory and LaSalle's principle [21], we prove the global stability of the DFE, which implies that the CLR will dwindle until extinction, whatever the initial number of urediniospores and infectious branches.
Proof. Consider the function where a 1 , b 1 , c 1 and d 1 are positive constants to be chosen later. Consider then the following Lyapunov function candidate: where Ω is the spatial domain defined above. Let us consider the set V = C (0, ∞); (3.16) Now, using the fact that S N 1, equation (3.16) becomes: (3.17) Then, the positive constants a 1 , b 1 , c 1 and d 1 are chosen such that: Note that R Then equation (3.15) becomes: From the Neumann bounded condition ∂U ∂η = 0, one has that: Thus, R  [35,42]. Replacing of the value U = 0 into the first equation of subsystems (3.1) yields: The above equation (3.22) has a unique equilibrium S 0 p = Λp µp , which is GAS. Hence, since the solutions of subsystems (3.1) are bounded, S(x, t) → S 0 p . Thus, one can conclude that the DFE Q 0 p = (S 0 p , 0, 0, 0, 0) is globally asymptotically stable for subsystems (3.1). This concludes the proof.

Comparison of the dynamics in the dry and rainy seasons
The objective of this section is to compare the dynamics of subsystems (3.1) during the rainy and dry seasons. To this end, we first compare the basic reproduction numbers R Simple biological hypotheses give the relation between the parameter values during the dry and rainy seasons. They mainly rely on the fact that some mechanims require humid conditions: coffee tree growth is faster during the rainy season (Λ D Λ R ), germination occurs mostly during the rainy season (ω D ω R ), the lesion progresses faster during the rainy season (θ D θ R ) and the mortality rate of the branches is higher during the dry season (as a consequence of harvest, which always damages some branches) i.e (µ R µ D ). Hence, we have: (3.23)

Comparison of the basic reproduction numbers
Herein we compare the basic reproduction numbers of the subsystems: Using (3.23), one has: This means that the GAS of the DFE Q 0 R during the rainy season implies the GAS of the DFE Q 0 D during the dry season and the existence and local asymptotic stability of the endemic equilibrium Q * D during the dry season implies the existence and local asymptotic stability of the endemic equilibrium Q * R during the rainy season. Biologically speaking, if the disease is present during the dry season, it will be also present during the rainy season while, if the disease is absent during the rainy season, it will be also absent during the dry season.

Comparison of the infectious branches at the endemic level
To compare the disease at the endemic level, both basic reproduction numbers need to be greater than one, 1 < R (p) 0 for p = D and R. The expression of I * p depends on the parameters identified in (3.23). The partial derivatives of I * p with respect to parameters Λ p , ω p , θ p and µ p are: This shows that I * p is an increasing function if Λ p , ω p , θ p and a decreasing function of µ p . Now, using (3.23), we can deduce that: Therefore, if it reached equilibrium, the CLR would be more severe during the rainy season than during the dry season.

Numerical simulations
Herein, we present the results of numerical simulations of system (2.1) using a non standard finite difference method [23]. We take T = 365 days and τ = 120 days which correspond to the durations of the year and the dry season in Cameroon, respectively. We set Ω = [0, 100] meters. We consider that initially all branches are healthy, that urediniospores are concentrated in the middle of the plantation and that there are no berries, so   Table 1 summarizes the parameter values for system (2.1) used for numerical simulations. We suppose that γ p and α p have the same value during the dry and rainy seasons.

CLR dynamics
We first choose γ p = 1.5 for p = D and R, so that R   Table 1. Subplots represent healthy branches S, latent branches L, infectious branches I, leafless branches J, urediniospores U and berries B.
Then we choose γ p = 8 for p = D and R, so that R (R) 0 = 5.09 and R (D) 0 = 2.77. In this case, we have shown in Lemma 3.8 that the endemic equilibria exist and is locally asymptotically stable during both seasons. Numerical results are shown in Figure 3. One can observe that the peaks of the infectious branches and urediniospores flatten very quickly compared to the other variables. The convergence towards the endemic stationary solution is clear at the beginning of the third year (after 1000 days approximately). Oscillations are observed, as the endemic equilibrium is higher during the rainy season (subsystem (3.1) with p = R) than the dry season (subsystem (3.1) with p = D). Note that urediniospores hardly reach the edges of the plantation during the first year, so that healthy branches are produced in large quantity, contrary to the middle of the plantation where urediniospores directly infect branches, which limits growth. Hence, when the urediniospores reach the edges during the second year, there is a large infection peak at the edges, which contrasts with the middle of the plantation. Figure 4 is a temporal representation of Figure 3 for three different values of x taken in the centre of the domain (x = 50 m, yellow curve), on the border (x = 1 m, blue curve) and in-between (x = 25 m, red curve). Convergence towards the stationary endemic solution can also be observed here, as all curves overlap for the three x values at the end of the third year and on. At the centre of the domain, the system very rapidly converges towards this stationary solution, but it takes a year for x = 25 m. On the border of the domain, branches remain susceptible during almost two years. Then the number of latent L and infectious I branches peak, as urediniospores U , initially in the centre of the domain, reach its borders. These peaks occur during the  Fig. 3). All other parameter values are given in Table 1. x is chosen in the centre of the domain (yellow curve), on the border (blue curve) and in-between (red curve). Subplots represent healthy branches S, latent branches L, infectious branches I, leafless branches J, urediniospores U and berries B.
second rainy season. However, no such peak is observed for leafless branches J. Indeed, the number of leafless branches starts to increase at the end of the second rainy season (around time t = 600 days), but the arrival of the dry season reduces the number of branches before it can peak very high. Solutions then oscillate between the endemic equilibria of the dry and rainy seasons (purple and green dashed lines, respectively). Infection has a negative impact on the production, which is reduced to less than a third of its value without disease.
These simulations allow to conclude that the dynamics of system (2.1), can be inferred from the dynamics of subsystems (3.1) during the dry and rainy seasons.
Biological control is a potentially powerful tool for managing coffee leaf rust that allows organic certification. Threfeore, we introduced in this study a mycoparasite (a parasitic fungus whose host is another fungus), Lecanicillium lecanii, known to hamper the reproduction of H. vastatrix. In our model, this biocontrol agent reduces the production of urediniospores by infectious branches with an efficiency q. The production rate of urediniospores γ p hence becomes (1 − q)γ p and system (3.1) becomes: The basic reproduction numbers of system (4.2), for p = D and R, are then: .
The higher the biocontrol efficiency q, the lower the reproduction numbers. So high enough values of q should allow to control the disease.
We simulated system (4.2) for different values of biocontrol efficiency q over 8 years. Figure 5 presents the temporal evolution of the state variables (S, L, I, J, U, B) integrated over domain Ω. Table 2 provides the reproduction numbers defined in equation (4.3), as well as the berry production and production loss during the 8th year, the latter being computed by reference to the disease-free case q = 1.
-q = 0 corresponds to no control (red curves) and R > 1, the number of spores U drops as expected, as the mycoparasite hampers the production of spores. So the number of healthy branches S increases. However, the number of latent L, infectious I and leafless J branches, barely change. Indeed, the production of latent branches depends on the product SU . The production of berries B then increases, as it is mostly driven by the healthy branches, but the yield loss is still high at 68.2%. -For a 70% efficiency (magenta curves), corresponding to R 0 , these observations are amplified: spores seriously drop; healthy branches and berries largely increase; latent, infectious and leafless branches decrease. The yield loss is 33.5%. -For a 75% efficiency (black curves), also corresponding to R 0 , infection remains very low, leading to an acceptable 6.0% yield loss.
-An 76% efficiency (blue curves), still corresponding to R 0 , almost yields the same results as a perfect efficiency (q = 1). Infection is almost negligible, so the yield loss is as low as 1.8%. -q = 1 (green curves) corresponds to a perfect (and unrealistic) efficiency, with R Hence a 75% biocontrol efficiency is enough to sustain the berry production in the plantation, with a negligible yield loss. Higher efficiencies, moreover, achieve disease eradication.  Table 1. Subplots represent healthy branches S, latent branches L, infectious branches I, leafless branches J, urediniospores U and berries B, these variables being integrated over domain Ω. Table 2. Impact of biocontrol efficiency q on reproduction numbers, defined in (4.3), and on berry production during the 8th and last year.

Conclusion
In this paper, we have proposed and analysed a PDE model that describes the dispersal of CLR in a coffee plantation during the rainy and dry seasons and its behaviour over time. Furthermore, we computed the diseasefree and endemic equilibria of the two subsystems defined during the rainy and dry seasons. We showed that the basic reproduction numbers during the two seasons can determine the dynamics of global model: when the basic reproduction number is less than one during the rainy season, then CLR globally decreases till extinction; when it is greater than one for the dry season, then CLR persists.
We implemented a biocontrol in our model, corresponding to a mycoparasite such as Lecanicillium lecanii, which hampers CLR reproduction at all times. This solution was tested in Mexico [11] but is still under development. A rather high biocontrol efficiency (75% at least) is necessary in our model to control the disease, but lower efficiencies still improve coffee production notably. Moreover, the mycoparasite is applied all year round, so it is not easily implemented in practice and it involves important costs. It would be interesting to study when to deploy the mycoparasite in a cost-efficient way. An ideal mycoparasite should sustain the dry season and efficiently control CLR, so that the coffee plantation during the rainy season would suffer reasonable yield losses. In further work, we will also include cultural management and other biocontrol agents, in particular natural endophytes which are affordable for growers starting a new plantation.
Several extensions to this work are considered: (i) adding a stage structure on the coffee branches; (ii) simplifying the model, using an impulsive formalism for the dry season as in [22,45], in order to obtain analytical results on the global model behaviour; (iii) solving an optimal control problem, consisting in maximising coffee production while minimising the control costs.

Appendix A.
Proof of Lemma 3.8.
-Local stability of the disease-free equilibrium Let (S, L, I, J, U ) be a solution of subsystems (3.1). Then, according to Kiehöfer [20], this solution can be written in the following form: With this in mind, subsystems (3.1) can be written in the following compact form: where D = diag(0, 0, 0, 0, 0, ε). The linearisation of subsystems (A.2) in the neighbourhood of Q 0 p is: where L(Q 0 p ) is the Jacobian matrix at the DFE Q 0 p of subsystems (3.1) in the absence of diffusion, that is: Let g j , j ∈ N, be the j th eigenfunction of operator −∆ with Neumann boundary conditions, so that: where λ j are the associated eigenvalues verifying 0 = λ 0 < λ 1 < λ 2 < . . .
According to [20], the expanded expression of W in equation (A.3) can be written as: where each Y j (t) ∈ R 5 . Substituting this expression into equation (A.3) yields: Thus, the DFE Q 0 p is stable if and only if each Y j (t) → 0 when t → ∞, that is, if and only if all the eigenvalues of matrix H j = L(Q 0 p ) − λ j D have negative real parts. Matrix H j can be written in the following form: The characteristic polynomial of matrix H j is given by Firstly one can observe that a 2 > 0, a 1 > 0 and Hence a 0 > 0 if γpθpωpν (θp+µp)(αp+µp+d)(ν+µ U +ελj ) < 1, which is satisfied for all (non negative) λ j if it is satisfied for λ 0 = 0. The condition a 0 > 0 then leads to: Secondly, the expression Using the fact that a 2 > 0, a 1 > 0, a 0 > 0 and a 2 a 1 > a 0 , the Routh-Hurwitz stability criterion indicates that all the eigenvalues of matrix H j have negative real parts. Hence the DFE Q 0 p of subsystems (3.1) is LAS if and only R (p) 0 < 1.
-Existence and local stability of the endemic equilibrium Q * p Suppose that R (p) 0 > 1. The expression (3.12) of the endemic equilibrium can easily be established as the solution of a set of linear equations, derived from equating the right-hand side of equation (3.1) to zero. We first prove that I * p and S * p are positive. Recall that: .
Note that: Now using the fact that θpd (θp+µp)(αp+µp+d) < 1, one has: This concludes the existence of the endemic equilibrium. Second, we investigate the local stability of the endemic equilibrium, using the following theorem.
Theorem A.1 (Castillo-Chavez and Song [12]). Consider the following ordinary differential equations, with a parameter ψ: Without loss of generality, it is assumed that 0 is an equilibrium for system (A.7) for all values of the parameter ψ, that is f (0, ψ) ≡ 0 for all ψ. Assume A 1 : A = D x f (0, 0) = ( ∂fi ∂xj (0, 0)) is the linearization matrix of system (A.7) around the equilibrium 0 with ψ evaluated at 0. Zero is a simple eigenvalue of A and all other eigenvalues of A have negative real parts; A 2 : Matrix A has a nonnegative right eigenvector u and a left eigenvector v corresponding to the zero eigenvalue. Let f k be the k th component of f and The local dynamics of (A.7) around 0 are totally determined by a and b.
Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.
In order to apply this theorem, subsystems (3.1) can be written as follows, with z = (S, L, I, J, U ): Solving R (p) 0 = 1, we obtain the following bifurcation value for parameter ω p : ω * p = (µ p + θ)(α p + µ p + d)(ν + µ U ) νγ p θ p . We linearise this system at the DFE Q 0 p , as previously in equation (A.5), setting ω p to ω * p . We need to determine the eigenvalues of matrix L ω * p − λ j D, where λ j is an eigenvalue of the Laplacian operator −∆ (simplified notation L ω * p is used instead of L ω * p (Q 0 p )). With λ j = λ 0 = 0, the matrix admits β 0 = 0 as eigenvalue, the other eigenvalues still having a negative real part. For the other λ j , all eigenvalues have negative real parts. So assumption A 1 of theorem A.1 is verified. Let us now verify assumption A 2 . We need to compute the left and right eigenvectors of matrix L ω * p − λ j D associated with the eigenvalue β. The left eigenvector, denoted by v = (v 1 , v 2 , v 3 , v 4 , v 5 ), satisfies the following equation: where I and 0 are the identity matrix and null vector of dimension 5, respectively. For β = β 0 = 0 one has: = 0 0 0 0 0 which gives: Similarly the right eigenvector of matrix L ω * p − λ j D, denoted by u = (u 1 , u 2 , u 3 , u 4 , u 5 ) T , satisfies the following equation: For β = β 0 = 0, one has: which gives: Let us now compute a, defined in assumption A 2 of theorem A.1, for system (A.8): The only terms that are non null correspond to: ω p ν * S 0 for i, j = 2, 3, 4.