Hopf Bifurcation in a Model for Biological Control
Jorge Sotomayor, Luis Fernando Mello, Danilo Braun Santos, Denis de Carvalho Braga
aa r X i v : . [ m a t h . D S ] A ug Hopf Bifurcation in a Model forBiological Control
Jorge Sotomayor
Instituto de Matem´atica e Estat´ıstica, Universidade de S˜ao PauloRua do Mat˜ao 1010, Cidade Universit´aria,CEP 05.508-090, S˜ao Paulo, SP, Brazil e–mail: [email protected]
Luis Fernando Mello
Instituto de Ciˆencias Exatas, Universidade Federal de Itajub´aAvenida BPS 1303, Pinheirinho, CEP 37.500-903, Itajub´a, MG, Brazil e–mail: [email protected]
Danilo Braun Santos
Centro de Ciˆencias Sociais e Aplicadas, Universidade MackenzieRua Itamb´e, 45, Consola¸c˜ao, CEP 01302-907, S˜ao Paulo, SP, Brazil e–mail: [email protected]
Denis de Carvalho Braga
Instituto de Sistemas El´etricos e Energia, Universidade Federal de Itajub´aAvenida BPS 1303, Pinheirinho, CEP 37.500-903, Itajub´a, MG, Brazil e–mail: [email protected]
Abstract
In this paper we study the Lyapunov stability and Hopf bifurcation in abiological system which models the biological control of parasites of orangeplantations.
Key-words : Hopf bifurcation, stability, periodic orbit, biological control.
MSC : 70K50, 70K20. Introduction of the Mathematical Model
In this work we study a system of four coupled differential equations (1)which models the interaction between two biological species, each presentingtwo stages in their metamorphosis, living in a common habitat with limitedresources.The differential equations analyzed here are P ′ = dPdt = φ (cid:18) − Mc (cid:19) M − ( α + β ) P − k P GM ′ = dMdt = α P − µ M (1) L ′ = dLdt = φ (cid:18) − Gc (cid:19) G − ( α + β ) L + k P GG ′ = dGdt = α L − µ G. This model —an elaboration of Lotka-Volterra equations, taking into ac-count the stages or compartments in the biological populations— was pro-posed by Yang and Ternes [1, 2] and Ternes [3] for a study of the biologicalcontrol of orange plantations leaf parasites P , which is a pre-adult stage for M , by their natural enemies L , which is an early stage for G .Other differential equations have been proposed as models for interactingpopulations partitioned in compartments, representing several situations ofbiological interest. See, among many others, Hethcote et al. [4], Jacquez andSimon [5] and Godfray and Waage [6].In [1, 2] and [3] P and M are the densities of pupae and female adultsof Phyllocnistis citrella (which in its larva stage is the citrus leafminer ), L and G are the densities of larvae and female adults of its native parasitoid http://en.wikipedia.org/wiki/Biological control http://en.wikipedia.org/wiki/Pupa http://en.wikipedia.org/wiki/Larva aleopsomyia fausta (whose larvae feed on the pupae of M ).The meaning of the parameters in (1), where the notation of [1, 3] hasbeen preserved, is as follows: α is the rate of pupae that give rise to adults M , β is the mortality rate of pupae, µ is the mortality rate of adults M , φ is the rate of eggs that give rise to pupae, c is the carrying capacity ofthe population M , α is the rate of larvae that, evolving through pupae, giverise to adults G , β is the mortality rate of larvae and pupae, µ is the rateof mortality of adults G , φ is the oviposition rate of the parasite and c isthe carrying capacity of the population G . Here we assume that the pupa(respectively larva) population decreases (respectively increases) at a rateproportional to P − G encounters that is k P G (respectively k P G ).This model represents the evolution of female populations. If necessarythe male populations can be estimated using the sexual ratio of each species.
Remark 1.1
All the parameters α , β , µ , φ , c , k , α , β , µ , φ , c , k arepositive. As the damage to the P population must be larger than the benefitto the L population it is natural to assume that k ≥ k . Here will be established the location and the stability character of theequilibria of (1), four in number. Also is determined the bifurcation varietyin the space of parameters, representing the transition from asymptoticallystable to saddle type at the equilibrium point with positive coordinates,representing the coexistence of the two species. See Theorem 2.4 and itsCorollary 2.5.Fixing the all the parameters in (1) to biologically feasible values, takenfrom [3] and [7], but letting the interaction coefficients k and k vary in apositive quadrant, the nature of the bifurcation phenomenon in this planeby crossing the bifurcation curve is established. See Theorem 3.1 and Figure1. This is done by means of a computer assisted calculation of the first en.wikipedia.org/wiki/Oviposition c is studied in Theorem 3.3 and illustrated in Figure 3.In Section 4 the implications of the results in this paper are discussed andinterpreted from a wider perspective. Assume the following notation: R = α φ µ ( α + β ) , R = α φ µ ( α + β ) . (2)The differential equations (1) have four equilibria A = ( P , M , L , G ) = (0 , , , , (3) A = ( P , M , L , G ) = (cid:18) c µ α (cid:18) − R (cid:19) , c (cid:18) − R (cid:19) , , (cid:19) , (4) A = ( P , M , L , G ) = (cid:18) , , c µ α (cid:18) − R (cid:19) , c (cid:18) − R (cid:19)(cid:19) , (5)and A = ( P , M , L , G ) , (6)where P = c µ φ α φ φ + µ c c k k (cid:18) α φ (cid:18) − R (cid:19) − µ c k (cid:18) − R (cid:19)(cid:19) ,M = c α φ α φ φ + µ c c k k (cid:18) α φ (cid:18) − R (cid:19) − µ c k (cid:18) − R (cid:19)(cid:19) ,L = c µ α φ α ( α φ φ + µ c c k k ) (cid:18) c µ k (cid:18) − R (cid:19) + α φ (cid:18) − R (cid:19)(cid:19) ,G = c α φ α φ φ + µ c c k k (cid:18) c µ k (cid:18) − R (cid:19) + α φ (cid:18) − R (cid:19)(cid:19) . emark 2.1 If R > and R > , then the equilibria A , A and A , haveonly non-negative coordinates. If k < k max , where k max = α φ (cid:16) − R (cid:17) c µ (cid:16) − R (cid:17) , (7) then the coordinates of the equilibrium A are also non-negative. The Jacobian matrix of (1) at x = ( P, M, L, G ) ∈ R has the form J ( x ) = − α − β − k G φ − φ Mc − k Pα − µ k G − α − β φ − φ Gc + k P α − µ , (8)while its characteristic polynomial is given by p ( λ ) = det ( J ( x ) − λI ) = Θ Θ + α ( µ + λ ) k k P G, (9)where Θ = ( µ + λ )( α + β + k G + λ ) − α (cid:18) φ − φ Mc (cid:19) and Θ = ( µ + λ )( α + β + λ ) − α (cid:18) φ − φ Gc + k P (cid:19) . Recall that an equilibrium point x is said to be a saddle of type n − p ifthe Jacobian matrix J ( x ) has n eigenvalues with negative real parts and p eigenvalues with positive real parts. Theorem 2.2 If R > , R > and k < k max then:1. The equilibrium A is a saddle of type 2-2;2. The equilibrium A is a saddle of type 3-1;3. The equilibrium A is a saddle of type 3-1. roof. From (9) the eigenvalues of J ( A ) are given by λ = − ( α + β + µ ) + q ( α + β + µ ) + 4 α φ [1 − R ] ,λ = − ( α + β + µ ) − q ( α + β + µ ) + 4 α φ [1 − R ] ,λ = − ( α + β + µ ) + q ( α + β + µ ) + 4 α φ [1 − R ] ,λ = − ( α + β + µ ) − q ( α + β + µ ) + 4 α φ [1 − R ] , and satisfy: λ > λ < λ > λ <
0. This proves the firstassertion.From (9) the eigenvalues of J ( A ) are given by λ = −
12 ( α + β + µ ) + 12 r ( α + β + µ ) − α φ [1 − R ] ,λ = −
12 ( α + β + µ ) − r ( α + β + µ ) − α φ [1 − R ] ,λ = − ( α + β + µ )+ q ( α + β + µ ) + 4 α φ [1 − R ] + 4 c α µ α [1 − R ] k ,λ = − ( α + β + µ ) − q ( α + β + µ ) + 4 α φ [1 − R ] + 4 c α µ α [1 − R ] k . Is immediate to see that λ > λ <
0. If φ > α [( α + β + µ ) + 4 µ ( α + β )]then λ and λ are complex with negative real parts and if φ ≤ α [( α + β + µ ) + 4 µ ( α + β )]then λ < λ <
0. This proves the second assertion.6rom (9) the eigenvalues of J ( A ) are given by λ = − ( α + β + µ + c k [1 − R ])+ q ( α + β − µ + c k [1 − R ]) + 4 α φ ,λ = − ( α + β + µ + c k [1 − R ]) − q ( α + β − µ + c k [1 − R ]) + 4 α φ ,λ = −
12 ( α + β + µ ) + 12 r ( α + β + µ ) − α φ [1 − R ] ,λ = −
12 ( α + β + µ ) − r ( α + β + µ ) − α φ [1 − R ] . It follows that λ > λ <
0. If φ > α [( α + β + µ ) + 4 µ ( α + β )]then λ and λ are complex with negative real parts and if φ ≤ α [( α + β + µ ) + 4 µ ( α + β )]then λ < λ <
0. This proves the last assertion. (cid:4)
For the sake of completeness we state the following lemma which is aparticular case of the Theorem of Routh–Hurwitz. See [8], p. 62.
Lemma 2.3
The polynomial L ( λ ) = a λ + a λ + a λ + a λ + a , a > ,with real coefficients has all roots with negative real parts if and only if thenumbers a , a , a , a are positive and the inequality ∆ = a a a − a a − a a > is satisfied. heorem 2.4 If R > , R > and k < k max then all the coefficients ofthe characteristic polynomial of J ( A ) are positive. Therefore, if ∆ = a a a − a − a a > , (10) where a = α + β + µ + α + β + µ + k G ,a = α φ c M + α φ c G + ( α + β + µ + k G )( α + β + µ ) ,a = ( α + β + µ + k G ) α φ c G + ( α + β + µ ) α φ c M + α k k P G ,a = α α φ (cid:18) k (cid:18) − R (cid:19) P + φ c (cid:18) − R (cid:19) M (cid:19) , then the differential equations (1) have an asymptotically stable equilibriumpoint at A . If ∆ < then A is unstable. Proof.
From (9) the characteristic polynomial of J ( A ) is given by[ λ + ( α + β + µ + k G ) λ + α φ c M ][ λ + ( α + β + µ ) λ + α φ c G ]+ α µ k k P G + α k k P G λ, which can be written as λ + λ [ α + β + µ + α + β + µ + k G ]+ λ [ α φ c M + α φ c G + ( α + β + µ + k G )( α + β + µ )]+ λ [( α + β + µ + k G ) α φ c G + ( α + β + µ ) α φ c M + α k k P G ]+ α φ c M α φ c G + α µ k k P G . Now it is simple to see that the coefficients of the characteristic polynomialare given by a , a , a , a above. From the hypotheses these coefficients arepositive. The stability at A follows from Lemma 2.3.8 The following corollary is immediate from the fact that a i > Corollary 2.5
The Jacobian matrix J ( A ) has a pair of complex eigenvalueswith zero real part if and only if a − a a a + a a = 0 , (11) where a i are defined in Theorem 2.4. In next section we study the stability of A under the condition (11),complementary to the range of validity of Theorem 2.4. The study outlined below is based on the approach found in the book ofKuznetsov [9], pp 175-178.Consider the differential equations x ′ = f ( x , µ ) , (12)where x ∈ R and µ ∈ R m is a vector of control parameters. Suppose (12)has an equilibrium point x = x at µ = µ and represent F ( x ) = f ( x , µ ) (13)as F ( x ) = A x + 12 B ( x , x ) + 16 C ( x , x , x ) + O ( || x || ) , (14)where A = f x (0 , µ ) and B i ( x , y ) = X j,k =1 ∂ F i ( ξ ) ∂ξ j ∂ξ k (cid:12)(cid:12)(cid:12) ξ =0 x j y k , (15)9 i ( x , y , z ) = X j,k,l =1 ∂ F i ( ξ ) ∂ξ j ∂ξ k ∂ξ l (cid:12)(cid:12)(cid:12) ξ =0 x j y k z l , (16)for i = 1 , , ,
4. Here the variable x − x is also denoted by x .Suppose ( x , µ ) is an equilibrium point of (12) where the Jacobian matrix A has a pair of purely imaginary eigenvalues λ , = ± iω , ω >
0, and noother critical (i.e., on the imaginary axis) eigenvalues.Let p, q ∈ C be vectors such that Aq = iω q, A ⊤ p = − iω p, h p, q i = X i =1 ¯ p i q i = 1 . (17)The two dimensional center manifold can be parameterized by w ∈ R = C , by means of x = H ( w, ¯ w ), which is written as H ( w, ¯ w ) = wq + ¯ w ¯ q + X ≤ j + k ≤ j ! k ! h jk w j ¯ w k + O ( | w | ) , with h jk ∈ C , h jk = ¯ h kj .Substituting these expressions into (12) and (14) one has H w ( w, ¯ w ) w ′ + H ¯ w ( w, ¯ w ) ¯ w ′ = F ( H ( w, ¯ w )) . (18)The complex vectors h ij are to be determined so that equation (18) writesas follows w ′ = iω w + 12 G w | w | + O ( | w | ) , with G ∈ C .Solving the linear system obtained by expanding (18), the coefficients ofthe quadratic terms of (13) lead to h = − A − B ( q, ¯ q ) , (19) h = (2 iω I − A ) − B ( q, q ) , (20)where I is the unit 4 × w ¯ w , whose coefficient satisfies a singular system for h ( iω I − A ) h = C ( q, q, ¯ q ) + B (¯ q, h ) + 2 B ( q, h ) − G q, (21)which has a solution if and only if h p, C ( q, q, ¯ q ) + B (¯ q, h ) + 2 B ( q, h ) − G q i = 0 . Therefore G = h p, C ( q, q, ¯ q ) + B (¯ q, (2 iω I − A ) − B ( q, q )) − B ( q, A − B ( q, ¯ q )) i , (22)and the first Lyapunov coefficient l – which decides by the analysis of thirdorder terms at the equilibrium its stability, if negative, or instability, if pos-itive – is defined by l = 12 ω Re G . (23)A Hopf point ( x , µ ) is an equilibrium point of (12) where the Jacobianmatrix A has a pair of purely imaginary eigenvalues λ , = ± iω , ω > transversal if the curves of complex eigenvaluescross the imaginary axis with non-zero derivative.In a neighborhood of a transversal Hopf point with l = 0 the dynamicbehavior of the system (12), reduced to the family of parameter-dependentcontinuations of the center manifold, is orbitally topologically equivalent tothe complex normal form w ′ = ( γ + iω ) w + l w | w | , (24) w ∈ C , γ , ω and l are smooth continuations of 0, ω and the first Lyapunovcoefficient at the Hopf point [9], respectively. When l < l >
0) a familyof stable (unstable) periodic orbits can be found on this family of centermanifolds, shrinking to the equilibrium point at the Hopf point.11 .2 Hopf Bifurcation in the Biological Model
In this subsection we analyze the stability at A given by (6) under thecondition (11). From (12) write the Taylor expansion (14) of f ( x ). Thus A = − ( α + β ) − k G φ (1 − M c ) 0 − k P α − µ k G − ( α + β ) φ (1 − G c ) + k P α − µ (25)and, with the notation in (14) to (16), one has F ( x ) − A x = (cid:18) − φ M c − k P G, , − φ G c + k P G, (cid:19) . (26)From (14), (15), (16) and (26) one has B ( x , y ) = ( B ( x , y ) , , B ( x , y ) , , (27)where B ( x , y ) = − φ c x y − k ( x y + x y ) ,B ( x , y ) = − φ c x y + k ( x y + x y ) , and C ( x , y , z ) ≡ . (28)To pursue the analysis consider the following table of specific parameters α = 0 . β = 0 . µ = 0 . φ = 2 . c = 400000 α = 0 . β = 0 . µ = 0 . φ = 4 c = 100 (29)taken from [7] and [3], where their biological feasibility in Brazilian fields isdiscussed.With the above parameter values the differential equations (1) are in facta two parameter system of differential equations where the parameters are k and k and can be written equivalently as x ′ = f ( x , k , k ) , (30)12ith f ( x , k , k ) defined by the right-hand sides of (1).With the parameter values of table (29), the equilibrium point A (6) hasthe following coordinates P = 800000 (1 . − . k )4 .
508 + 1 . · k k , M = 933333 .
333 (1 . − . k )4 .
508 + 1 . · k k ,L = 444 .
444 (1 .
216 + 85550 . k )4 .
508 + 1 . · k k , G = 333 .
333 (1 .
216 + 85550 . k )4 .
508 + 1 . · k k , while R , R and k max , given by (2) and (7), have the form R = 3 . , R = 9 . , k max = 0 . . (31)From the above equation and the Remark 1.1, the set of admissible parame-ters is given by (see Fig 1) S = { ( k , k ) | < k < k max = 0 . < k ≤ k } . (32)In this set S the curve Σ = ∆ − (0) is well-defined (see (11)), where ∆ isgiven by1699 . − . k − . · k − . · k − . · k − . · k k − . · k − . · k k − . · k k − . · k k − . · k k − . · k k − . · k k − . · k k − . · k k + 2 . · k k − . · k k +1 . · k k + 6 . · k k − . · k k + 1 . · k k , representing the parameters where J ( A ) has a pair of purely imaginaryeigenvalues λ , = ± iω with ω = 1 . h . . − . · − + k )(1 . · − + k ) k (3 . · − + k k ) +8 . · − k + (cid:16) . · − + k k ) k k (cid:16) (1 . · − + (33)5 . · − k ) k − . · − (1 . · − k )(6 . · − + k ) +0 . k (9 . · − + k )(8 . · − + k ) (cid:17)(cid:17) / i / . Set of admissible parameters S and Hopf curve Σ .Thus one has (see Fig. 1) S = S + ∪ Σ ∪ S − . For parameter values in the region S + the equilibrium A is unstable sincethe Jacobian matrix J ( A ) has two complex eigenvalues with positive realparts and two other real negative eigenvalues. For parameter values in theregion S − the equilibrium A is asymptotically stable since J ( A ) has twocomplex eigenvalues with negative real parts and two other real negativeeigenvalues. The curve Σ is the curve of admissible parameters where theequilibrium A is a Hopf point. Theorem 3.1
Consider the differential equations (1) with the parametersgiven in the table (29). If ( k , k ) ∈ Σ then the two parameter family ofdifferential equations (1) has a transversal Hopf point at A . This Hopfpoint at A is unstable and for each ( k , k ) ∈ S − , but close to Σ , there existsan unstable periodic orbit near the asymptotically stable equilibrium point A .See Fig 1. omputer Assisted Proof. The proof follows the steps outlined in Sub-section 3.1. However, all the expressions in the proof are too long to be putin print. For this reason, in the site [10] have been posted the main stepsof the long calculations involved in the proof. This has been done in theform of a notebook for MATHEMATICA 5 [11]. A sufficient condition forbeing a Hopf point is that the first Lyapunov coefficient l = 0. In fact, itcan be shown numerically that l ( k , k ) > k , k ) ∈ Σ. Aparticular case and other related calculations are considered below for thesake of illustration.Take the particular point Q = ( k = 0 . , k = 0 . ∈ Σ with fivedecimal round-off coordinates [7]. For these values of the parameters A = (18543 . , . , . , . . The Jacobian matrix J ( A ) has eigenvalues λ = − . , λ = − . , λ , = ± . i, and thus ω = 2 . . (34)From (17) the eigenvectors q and p have the form q = . . i . − . i . . i . − . i ,p = . . i − . . i . . i . − . i . One has B ( q, q ) = − . . i . . i B ( q, ¯ q ) = . − . . · − i . From (19) and (20) the complex vectors h and h have the form − h = − . . · − i − . . · − i − . − . · − i − . − . · − i ,h = . . i . − . i . − . i − . − . i . From (22) the complex number G is given by G = 0 . − . i, (35)and from (23), (34) and (35) one has the first Lyapunov coefficient at Ql ( Q ) = 0 . > . (36)The above calculations have also been checked with 10 decimals round-offprecision performed using the software MATHEMATICA 5 [11]. See [10].Some other values of pairs ( k , k ) ∈ Σ, the values of the complex eigen-values of J ( A ) as well as the corresponding values of l ( k , k ) are listed thetable below. The calculations leading to these values can be found in [10].16 k complex eigenvalues of J ( A ) l ( k , k )0.0004813 0.0004812 ± . i . · − ± . i . · − ± . i . · − ± . i . · − ± . i . · − ± . i . · − ± . i . · − ± . i . · − ± . i . · − ± . i . · − ± . i . · − (cid:4) Remark 3.2
The value of the first Lyapunov coefficient l does depend onthe normalization of the eigenvectors q and p , while its sign is invariantunder scaling of q and p obeying the relative normalization. See [9], p. 99.The values l ( Q ) in (36) and l ( k , k ) in the above Table are obtained withtwo different choices of the eigenvectors q and p , see [10]. This explains thedifference in the order of magnitude of the numbers involved. As a consequence of Theorem 3.1 there are no Hopf points of codimension2 on Σ since the sign of the first Lyapunov coefficient does not change. InFig. 2 is illustrated the bifurcation diagram for a typical point on the curveΣ. Assuming the same values in the table (29) in next theorem we studythe behavior of the Hopf curve Σ in the set of admissible parameters S (seeequation (32)) as the parameter c increases. In fact, the carrying capacity,representing several other factors, has a determinant role on the populationsunder study. Theorem 3.3
The one parameter family of curves Σ c = ∆ − c (0) has onlyone point of tangency T with the line k = k for c = 650 . . For values Bifurcation diagram for a typical point on the curve Σ . c > . the curve Σ c does not intersect the set S . Therefore forvalues c > . the set S + is empty, S = S − and the equilibrium A is asymptotically stable for all values ( k , k ) ∈ S . See Fig. 3. Proof.
The surface of Hopf points, or equivalently the one parameter familyof Hopf curves, where J ( A ) has a pair of purely imaginary eigenvalues isdefined by Σ c = { ∆( k , k , c ) = 0 } where ∆( k , k , c ) (see (11)) is given by1699 . − . c k − . · k − . c k − . · k − . · c k k − . c k − . · c k k − . · c k k − . · c k k − . · c k k − . · c k k − . · c k k − . · c k k − . · c k k +2 . · c k k − . · c k k + 1 . · c k k +6 . · c k k − . · c k k + 1 . · c k k . Curve Σ intersects S at one point T .The intersection of the surface Σ c with the plane k = k determines thecurve C , given implicitly by N ( k , c ) = 1699 . − (2337 . c + 6 . · ) k − (1114 . c +2 . · + 4 . · c ) k − (214 . c + 4 . · c +1 . · c ) k − (7 . · c + 4 . · c + 1 . · c ) k − (4 . · c + 1 . · c + 6 . · c ) k + (2 . · c − . · c + 1 . · c ) k + (6 . · c − . · c ) k +1 . · c k = 0 . Differentiating implicitly the above expression with respect to k one has dc dk = − ∂N∂k ∂N∂c = 0 , d c dk < , at k = 0 . c = 650 . C is a graph nearthe point ( k = 0 . , c = 650 . k = 0 . dc /dk has no other zero. It is easy to verify through a calculation that the point19 = ( k , k ) = (0 . , . c for c = 650 . c at T for c = 650 . − . · , . · ) , which is parallel to the vector ( − , k = k . (cid:4) Remark 3.4
Since k max does depend on the parameter c , according to Eq.(7), so does the admissible region S = S c . For c = 650 . , a calculationgives k max = 0 . . This is compatible with position of T at k = k =0 . , as illustrated in Figure 3. In this paper we studied the system (1) of interest as a mathematical modelfor biological control, proposed by Yang and Ternes [1, 3, 2] and studied alsoby Santos [7]. Valuable field data are provided in [3], valid for the citrusleafminer and its native and imported enemies in the region of Limeira (S˜aoPaulo, Brazil). An extensive, enlightening discussion of the economic andagricultural interest of the problem, other pertinent differential equationsmodels as well as extensive bibliography, are also presented there.Under conditions made explicit in Remark 2.1 we determine the uniqueequilibrium point ( A ) with positive coordinates and establish necessary andsufficient conditions for its (Lyapunov) stability (Theorem 2.4). It can beseen however that this condition ∆ >
0, when expressed in terms of theparameters is a rational function whose denominator does not vanish and itsnumerator is a polynomial of too many terms to be put in print, but stillamenable to numerical calculations. For this reason the treatment of thestability of ( A ) in Subsection 3.2 is computer assisted. The conclusion ofthis study, made precise in Theorem 3.1, is the existence of periodic orbitsobtained by Hopf bifurcation, on the side (of ∆ = 0) where A is an attractor.20he study of the general analytic and geometric properties of the bound-ary of the stability region, given by the Hopf variety ∆ = 0, so as to includeparameter values of biological interest as proposed here as well as othersappearing in the work of Ternes [3], remain at the present moment as amathematical challenge. Theorem 3.3 gives only a thin slice of the geometry.The reports in [12] and [13], among many others, show that the interestfor the combat of the citrus leafminer extends to most regions where citrustrees grow.The Mathematica notebooks [10], with the table (29), used in the com-puter assisted arguments for the proofs of Theorems 3.1 and 3.3, can beadapted to tables with data pertinent to other geographic and climatic re-gions and involving different host–parasitoid interactions.In [2] Ternes and Yang discuss, with pertinent documentation, the intro-duction of a foreign parasitoid, Ageniaspis citricola to add the native
Galeop-somya Fausta in the combat with the leafminer, Phyllocnistis citrella . Theypropose a model with eight differential equations for the three species andtheir immature stages. In [2] and [3] are given starting steps for an analysisof the stability of the equilibria in this extended eight – dimensional system.Based in a numerical study of a complex equilibrium point they recommendthat the biological control of the leafminer be implemented with both the na-tive and foreign parasitoids. Meanwhile, the Hopf bifurcation analytic andcomputer algebra study of the complex equilibria of the eight equations, withthe methods used in the present paper, seems unsurmountable at the presentmoment, due to the large number of parameters involved.
Acknowledgement : The first and second authors developed this work un-der the projects CNPq Grants 473824/04-3 and 473747/2006-5. The firstauthor is fellow of CNPq. The fourth author is supported by CAPES. Thiswork was finished while the second author visited Universitat Aut`onoma deBarcelona, supported by CNPq grant 210056/2006-1.21 eferences [1] H. M. Yang and S. Ternes, Estudo dos efeitos de dinˆamica vital nummodelo de controle biol´ogico de pragas (Study of the effect of the vitaldynamics in a model for the biological control of plagues).
Revista deBiomatem´atica , 58-72 (1999) (in Portuguese). [2] H. M. Yang and S. Ternes, Um modelo determin´ıstico para avalia¸c˜aodo controle biol´ogico de praga de citros (A deterministic model for theevaluation of the biological control in plagues of citrus). Boletim dePesquisa e Desenvolvimento, Embrapa , 1-25 (2002) (in Portuguese). [3] S. Ternes, Modelagem e simula¸c˜ao da dinˆamica populacional da larva-minadora-da-folha-de-citros em intera¸c˜ao com seus inimigos naturais(Modelling and simulation of population dynamics of the leafminer ininteraction with its natural enemies) , Tese de Doutorado, Faculdadede Engenharia El´etrica e de Computa¸c˜ao, Unicamp, Campinas, Brazil(2001) (in Portuguese). http://libdigi.unicamp.br/document/?code=vtls000228128 [4] H. W. Hethcote, Y. Li and Z. Jing, Hopf bifurcation in models for per-tussis epidemiology,
Math. Comput. Modelling , 29-45 (1999).[5] J. A. Jacquez and C. P. Simon, Qualitative theory of compartmentalsystems, SIAM Review , 43-79 (1993).[6] H. C. J. Godfray and J. K. Waage, Predictive modelling in biologicalcontrol: the mango mealy bug ( Rastrococcus invadens ) and its para-sitoids,
J. Appl. Ecol. , 434-453 (1991).[7] D. B. Santos, Bifurca¸c˜ao de Hopf num modelo de controle biol´ogico (HopfBifurcation in a biological control model) , Disserta¸c˜ao de Mestrado, In-22tituto de Matem´atica e Estat´ıstica, Universidade de S˜ao Paulo, S˜aoPaulo, Brazil (2004) (in Portuguese). [8] L. S. Pontryagin,
Ordinary Differential Equations , Addison-Wesley Pub-lishing Company Inc., Reading (1962).[9] Y. A. Kuznetsov,
Elements of Applied Bifurcation Theory , second edi-tion, Springer-Verlag, New York (1998).[10] Site with the files used in computer assited arguments in this work: [11] S. Wolfram,
The Mathematica Book , fifth edition, Wolfram Media Inc.,Champaign (2003).[12] Workshop on the Citrus Leafminer and its Control in the Near East,F.A.O., Syria, Oct. (1996).