- Email: [email protected]
Hopf bifurcation in a grazing system with two delays
Accepted Manuscript Hopf bifurcation in a grazing system with two delays A. Mendy, J.J. Tewa, M. Lam, P. Tchinda Mouofo
PII: DOI: Reference:
S0378-4754(19)30054-0 https://doi.org/10.1016/j.matcom.2019.02.006 MATCOM 4740
To appear in:
Mathematics and Computers in Simulation
Received date : 30 November 2017 Revised date : 26 December 2018 Accepted date : 11 February 2019 Please cite this article as: A. Mendy, J.J. Tewa, M. Lam et al., Hopf bifurcation in a grazing system with two delays, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.02.006 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript Click here to view linked References
Hopf bifurcation in a grazing system with two delays A. Mendya , J. J. Tewab , M. Lama , P. Tchinda Mouofob,1 a
University Cheikh Anta Diop, Department of Mathematics and Computer Science, Faculty of Science and Technic, P.O. Box 5005 Dakar b Department of Mathematics and Physics, National Advanced School of Engineering, University of Yaound´e I, P.O. Box 8390 Yaound´e, Cameroon
Abstract Life in the Sahelian zone is very difficult due to lack of resources on a permanent basis. This scarcity of resources leads the pastor, with their flock to organize themselves to better resist the difficult conditions. Pastoralism is of paramount importance for this region. There is an urgent need to understand this way of life in order to help pastors to better organize themselves. In this paper, we propose mathematical models that consist of three food trophic chains. We consider and formalize, in a given location, interactions between browsers, grazers and forage resources. We first analyze the resulting system without delays and in a second time take into account two positive and discrete time-delays. In our model, delays denote an average time required in order that food consumed by browsers and grazers are beneficial for them. Qualitative analysis of the delays-free model shows several equilibria as well as various situations of multistability. Considering delay as parameter, we investigate the effect of delay on the stability of equilibria. It is found that the delays can lead the system dynamic behavior to exhibit stability switches, and Hopf bifurcation occurs when the delay crosses some critical value. By applying the normal form theory and the center manifold theorem, the explicit formula which determine the stability and direction of the bifurcating periodic solutions are determined. Finally based on a nonstandard numerical scheme, we provide some numerical simulations in order to illustrate our qualitative results and support discussions. Keywords: Pastoralism, Browsers, Grazers, Hopf Bifurcation, Multistability. 1. Introduction Pastoralism describes the interdependent relationship between pastoralists, their ruminant herds and their habitat. The exploitation of livestock remains strongly dominated by extensive methods of herd management. The pastoral system is dominant along the Sahelian belt with relatively large herds of cattle and small ruminants (Report 2010 [26]). It is based on extensive exploitation of natural resources without recourse to zootechnical inputs, with the exception of years with critical forage deficit (see Report 2010 [26]). This type of livestock relies on the mobility or transhumance of herders and herds in search of water, pastures (during the long period of dry season, 8 to 10 months), and areas of salt cures (see also Report 2010 [26]). The dynamics of change in livestock activities are at the heart of Email addresses: [email protected] (A. Mendy), [email protected] (J. J. Tewa), [email protected] (M. Lam), [email protected] (P. Tchinda Mouofo) 1 Corresponding author: +(237) 679 44 12 29
Preprint submitted to Mathematics and Computers in Simulation
December 26, 2018
the challenges of sustainable development of many territories in the grassland and pastoral areas (Dedieu et al. 2009 [10]). The pastoral farming systems are those in which more than 90% of the dry matter consumed by livestock comes from grazing(Carriere 1996 [6]). Extensive pastoral production is practiced in 25% of the world’s land, from the arid zones of Africa (66% of the continent’s land) and the Arabian Peninsula to the tableland of Asia and Latin America (Nori et al. 2008 [24]). It provides 10% of the world’s meat production and supports some 200 million pastoralist households that raise almost 1 billion head of camels, cattle and small ruminants, about one third of them in sub-Saharan Africa, where pastoral farming represents about 20% of the Gross Domestic Product (GDP) of the different states (Nori et al. 2008 [24]). In the Sahel, the environmental consequences are less obvious (Smith 1996 [33]). The great majority of food consumed by ruminants in the sahelian region is still formed of natural pasture (Carriere 1996 [6], Research reports and studies 2009 [29]). Despite this important role, pastoral systems face many challenges related notably to the marginalization of pastoralists and to the increasing competition with other user groups especially farmers for accessing to natural resources (pasture, water points etc.) (see Report 2012 [27]). Recent studies indicate that the characteristics of the food and animal’s body size, rather than anatomical adaptations for grazing and browsing, have greater significance for the nutritional ecology of the herbivores (Shipley 1996 [30]). To better manage range livestock and habitats of free-ranging herbivores, ecologists and range managers often wish to estimate their forage intake rate. Morphological characteristics of grasses makes it easier to estimate the intake rate of grazers than that of browsers (Shipley 1996 [30]). The relative role of processes such as competition, facilitation and predation in structuring herbivore assemblages are inherently difficult to ascertain, and this has caused considerable controversy (McNaughton 1984 [23], Sinclair 1985 [31], Sinclair and Northon-Griffiths 1982 [32]). Inter-specific competition arises whenever one species reduces shared food resources below the level that can be exploited efficiently by another species (Bell 1971 [4]). However, by reducing plant biomass one species may also benefit another by facilitating access to forage of a suitable height or quality (Bell 1971 [4], McNaughton 1984 [23]). Long-term sustainable pastoral development requires a good knowledge of the dynamics of multiple factors underlining pastoralism, and here research has a crucial role to play. The links between research and policy also need to be strengthened, so that policy responds to the dynamics of pastoral livelihoods in an efficient manner (African Union 2010 [1]). The issue of livestock impacts on the dryland environment generally refers to a number of major themes, the main ones being desertification, drought, overgrazing, sedentarization of livestock breeders, rights of access to livestock resources (Carriere 1996 [6]). These highly mediatized subjects are the matter of intense controversy between specialists because, depending on whether an economic, sociological or ecological approach is adopted, the impacts are assessed in a different, even contradictory way (Carriere 1996 [6]). But it was also interesting to observe the evolution of pastoral areas according to the modalities of their exploitation by the animals. The demand for resources for herds depends on the number and type of animals and their nutritional requirements. The question of the resource is very important throughout the world and particularly in the Sahel, which has a wet season (june to october) and a dry season (november to may) (Tewa et al. 2014 [35]). An important issue in this part of Africa is the management of herds by pastoralists in order to maintain diversity between browsers and grazers in the sahel as well as manage the coexistence between resource and herbivores. 2
Mathematical models can help to better understand the action of pastoralists and to help them to better organize their work for an optimal yield. To authors knowledge, there is no mathematical model that addresses interactions between resource, browsers and grazers in Sahelian regions. In order to contribute to these problems related to Sahelian pastoralism, we therefore develop in this paper a model of interactions between the resource and a herd of herbivores in a given location. In our model, we make two important hypothesis. The first one is that as herbivory reduces the plant biomass of its maximum crop, the productivity of plants will increase. In fact, some empirical evidences show that when herbivory reduces plant biomass, for maintaining the ecosystem and avoid extinction, plants become resilient to herbivory and therefore increase their productivity (McNaughton 1984 [23], Lebon et al. 2014 [19]). The second hypothesis come from the fact that after consumption of resources, browsers and grazers need some time so that the consumed resource is beneficial to them. Taking the second hypothesis leads to delay differential systems and a rich dynamic of our model. This work is the first one which models interaction of herbivory and resource and taking time delays. The paper is organized as follows. We formulate a delayed mathematical model for herbivores and their resource interactions in section 2. In section 3, basics mathematical results are obtained and stability of equilibria are established when there is no delay in section 4. In section 5, the existence of Hopf bifurcation is established when delays are considered. Moreover, direction and stability of the Hopf bifurcation are also presented in section 5. Some numerical simulations that illustrate our analytical results are presented in section 6 and finally, a summary is presented in section 7. 2. Model construction 2.1. Dynamic of resource The resource includes grass, shrub and tree biomasses. The resource-herbivore interactions can be considered as a particular case of the well-mastered prey-predator system in which the logistic growth curve is most appropriate for writing primary productivity (NoyMeir 1975 [25], Caughley 1976 [7], Samuel et al. 1989 [38]). To understand the impact of herbivore pasture on the dynamics of the resource, several theories have been proposed including the Herbivores Optimization Hypothesis (HOH) (McNaughton 1984 [23], Hilbert et al. 1981 [17]). The HOH indicates that as the intensity of pasture increases, net primary productivity increases to optimum grazing intensity, then decreases to a lower level than the non-fattened resource (Samuel et al. 1989 [38]). This hypothesis of optimization of the grazing of herbivores was also used by Lebon et al. 2014 [19]. The resource growth function that models the HOH is defined by (see Lebon et al. 2014 [19]): ψ(B, H) = rB (H)B − δB 2 , with rB (H) =
H 1+ β
H 1− , γ
where B represents plant biomass and H herbivore (or pest) population. In this paper, we use this approach to model the growth of the resource in the presence of two populations of herbivores. We consider the following HOH-like function ψ(R, Hb , Hg ) = rφR (Hb , Hg )R − δR2 , 3
with φR (Hb , Hg ) =
αb Hb + αg Hg 1+ β
αb Hb + αg Hg 1− , γ
where αb , αg ∈ {0; 1}, γ > β and αb Hb + αg Hg < γ. The term δR2 is the density-dependent death rate of the resource as a consequence of local competition for water and nutrients. We note that in the herbivore-free case (i.e. ψ = 1), the resource follows a logistic growth (e.g. Tewa et al. 2014 [35], Zhilan et al. 2011 [11], Li et al. 2006 [20]). Primary production is heavily dependent on precipitation (Carriere 1996 [6]), which means that the intrinsic growth rate may be negative. The intrinsic growth rate has been defined as (see Tewa et al. 2014 [35]): r = r01 p + r02 − r2 , where p denotes the mean annual precipitation, r01 is the growth rate of the resource with precipitation, r02 is the growth rate of the resource without precipitation, r2 is the natural death rate. The particularity of this model, which is the difference with classic predator-prey models in the literature, is the fact that the intrinsic production rate can be negative (r < 0) (Tewa et al. 2014 [35]). This situation is frequently observed in Sub-Saharan Africa, particularly in Sahelian region when precipitations are scarce (Tewa et al. 2014 [35]). Intensification of this situation can conduct pastoralists to migrate with their herds: this is the transhumance phenomenon. 2.2. Dynamic of browsers Browsers include camels, beefs, giraffes, horses etc. (Johan et al. 1998 [8]). They feed on grass, shrubs and stress biomass. In the study of interactions between herbivores, and resource, it is important to determine which specific form of functional response, describing the amount of resource consumed by an animal in the herd per unit of time, is ecologically meaningful and provides a solid basis for theoretical development (Tewa et al. 2013 [37], Tewa et al. 2012 [36]). In the simplest case, such a function is a linear function of the forage resources density, which is integrated into the classic Lotka-Volterra model. The linear functional response is a limited case and can only be seen over a short period. However, the Holling type II or Holling type III functional response can be used as an improvement of the linear-like function. With the Holling type II function, at low resource biomass, resource experiences herbivory while with the Holling type III function herbivory is very low at low resource biomass. The Holling type III function is typical of generalist consumers moving from one resource to another and concentrating their activities in a region where the resource is abundant (Bates et al. 2016 [3]). This function is defined by: P (R) =
aR2 , b + R2
where a is the maximum consumption rate of browsers on the resource, b is the half-saturation parameter. As modelling hypothesis, we assume that there is no alternative food for herbivores. Consequently, in absence of resource biomass, both browsers and grazers go to extinction.
4
2.3. Dynamic of grazers Grazers include sheep, goats, etc. (Koppel et al. 1998 [8]). Grazers feed on grass and can not reach shrubs and trees. The interactions between the resource and grazers follow the functional response that takes into account the interspecific competition (competition between browsers and grazers). An interspecific competition is a competition between two different species for the same resource. The competition for the resource arises from the fact that the resource is limited and that the competitors interfere to obtain it. BeddingtonDeAngelis 1975 [9] proposed a function defined by: g(R, Hg , Hb ) =
c1 Hg R, 1 + c2 Hg + c3 Hb
where c1 , c2 and c3 are positive parameters. c1 is the catch rate, c2 is the manipulation time, and c3 is the interference amplitude between browsers and grazers. The term c1 Hg corresponds to the consumption of the resource by the small herbivores 1 + c2 Hg + c3 Hb taking into account its competition with the large herbivores (DeAngelis et al. 1975 [9]). The functional response of Beddington-DeAngelis suggests that the effects of interferences between herbivores become negligible when the density of the resource biomass is high. This function characterizes the fact that in addition to encountering the resource, herbivores likely collide with each other because the resource is limited. This function also assumes that herbivores spend their time in three activities: searching for the resource, consuming the resource and manipulating other herbivores. After consumption of the resource, browsers need a time so that the consumed resource is beneficial to them. This period is represented by τ1 . Similarly, for the consumption of the resource by grazers, this period is represented by τ2 . The following system of ordinary differential equations models interactions between forage resources, browsers and grazers : dR(t) dt
αb Hb (t) + αg Hg (t) αb Hb (t) + αg Hg (t) = r 1+ 1− R(t) − δR2 (t) β γ −
c1 Hg (t) aR2 (t) R(t) − Hb (t), 1 + c2 Hg (t) + c3 Hb (t) b + R2 (t)
dHb (t) aR2 (t − τ1 ) = e Hb (t) − (db + hb )Hb (t), b dt b + R2 (t − τ1 ) dHg (t) c1 Hg (t) = eg R(t − τ2 ) − (dg + hg )Hg (t), dt 1 + c2 Hg (t) + c3 Hb (t)
(1)
where, db and dg are the natural mortality rates of large herbivores and small herbivores respectively, hb and hg are the harvest rates of large herbivores and small herbivores respectively, eb and eg are the conversion rates of the resource consumed into large and small herbivorous biomass, respectively. System (1) is considered with the following initial values R(θ) = ϕ1 (θ) ≥ 0, Hb (θ) = ϕ2 (θ) ≥ 0, Hb (θ) = ϕ3 (θ) ≥ 0, θ ∈ [−τ, 0], where τ = max{τ1 , τ2 }, ϕ1 , ϕ2 , ϕ3 ∈ C [−τ, 0], R3+ , the Banach space of continuous functions mapping the interval [−τ, 0] into R3+ , with R3+ = {(R, Hb , Hg )/ R ≥ 0, Hb ≥ 0, Hg ≥ 0}. 5
3. Basic results 3.1. Boundedness of solutions We start by showing that solutions of System (1) that start in R3+ will remain there and are uniformly bounded. The following Lemma is valid. Lemma 1. 1) If r ≤ 0, then all the trajectories of the System (1) go to the zero solution. 2) If r > 0, (i) The nonnegative orthant R3+ is a positively invariant region for the System (1). (ii) Let ε > 0, the set Ωε defined by r γ 3 Ωε = (R, Hb , Hg ) ∈ R+ : 0 ≤ R ≤ 1+ + ε; 0 ≤ αb Hb + αg Hg ≤ γ δ β is a positively invariant and absorbing set for the System (1). Proof: See Appendix B.1. From Lemma 1, one deduces that System (1) is biologically well-posed, that is, solutions are non negative for non negative initial data and remain bounded. Moreover, we also deduce from Lemma 1 that System (1) is a dissipative system (Hale 1988 [15]). In the sequel we give the results concerning existence of equilibria for System (1). 3.2. Equilibria of the model We reach the step of characterizing equilibria of System (1). Let us set, ˜ ˜ rαb δR eb a β aR 4βγηa , ηs = , ηa = 1− , S01 = , , σR˜ = N0 = 2 2 ˜ βσR˜ γ r µb b+R 1 2 (γ − β) 1 − + 4βγ ηs δµg βγc2 + αg (γ − β) γ−β L0 = , L1 = and L2 = . 2βγδµg c2 βγc1 αg βγδµg c2 reg c1 + + reg c1 r c2 rαg eg c1 Proposition 1 (Existence of equilibria). The following results hold for System (1). r (1) Equilibrium E0 = (0, 0, 0) always exists while equilibrium E1 = , 0, 0 exists whenδ ever r > 0. ˜ H ˜ b , 0 of System (1), there exist two cases. (2) For the grazer-free equilibrium E˜ = R, ˜ H ˜ 0 , 0) with (a) When N0 = 1, we have a unique equilibrium E˜0 = (R, b s bµb ˜0 = 1 γ − β 1 − 1 . ˜= and H R b e b a − µb 2αb ηs E˜0 is ecologically meaningful whenever S01 > 1 and ηs > 1.
6
˜ H ˜ 1 , 0) and E˜2 = (R, ˜ H ˜ 2 , 0) with (b) When N0 < 1, we have two equilibria: E˜1 = (R, b b s 2 ˜ b1 = 1 (γ − β) 1 − 1 − (γ − β)2 1 − 1 H − 4βγ(ηa − 1) 2αb ηs ηs and
˜ b2 = 1 (γ − β) 1 − 1 H 2αb ηs
+
2 1 (γ − β)2 1 − − 4βγ(ηa − 1) . ηs
s
E˜1 is ecologically meaningful whenever N0 < 1, S01 > 1, ηa > 1 and ηs > 1 while E˜2 is ecologically meaningful whenever (S01 > 1 and ηa < 1) or (S01 > 1, ηa < 1 and ηs < 1). ¯ 0, H ¯ g ) of System (1) exists with (3) The browser-free equilibrium E¯ = (R, µ g ¯= ¯ g ) and H ¯ g is a positive solution of equation R (1 + c2 H eg c1 ¯ 3 + γ − β (1 − L2 ) H ¯ 2 + βγc2 + αg (γ − β) (1 − L1 ) H ¯ g + βγ (L0 − 1) = 0. H g g αg L 2 αg c2 L1 αg2 c2
(2)
(4) A coexistence equilibrium E ∗ = (R∗ , Hb∗ , Hg∗ ) exists if S01 > 1, N ∗ > 1, and Hg∗ being a positive solution of equation D2 Hg∗2 + D1 Hg∗ + D0 = 0,
(3)
r
eg c1 R∗ − µg (1 + c2 Hg∗ ) bµb , Hb∗ = and the expression of D2 , D1 e b a − µb c3 µ g and D0 are given in the proof. ∗
where, R =
Proof: See Appendix B.2. The complexity of results given in Proposition 1 calls for some remarks: first, concerning existence of browsers-free equilibrium, one may obtain up to two equilibria of that type thanks to the modelling option that consists on considering the HOH and BeddingtonDeAngelis’s function. A similar observation can be done for the grazer-free equilibrium as well as the coexistence equilibrium. Indeed, thank to the modelling choice of HOH and Holling type III functional response, on may obtain up to three grazer-free equilibria and up to two coexistence equilibria. This is not a common feature of predator-prey models. However in our case, it is desirable since it may lead to several cases of multistability i.e. more than one equilibrium is stable in the same time (see below). In the sequel, we deal with the problem of stability characterization of various equilibria of model (1). For that, we first consider the case where there is no delay and, secondly focus on the case where delays are considered. This process is motivated by the fact that if an equilibrium is unstable without delay, it will remain unstable when delays are considered (see Rebecca et al. 2000 [28] and Martin et al. 2001 [22]).
7
4. Stability analysis when τ1 = τ2 = 0 Let us set, ˜ eb a R˜ H˜b e g c1 R eb ar2 δµg 1 = = S0 = ; L = ; S ; S ; 0 0 0 µb (bδ 2 + r2 ) e g c1 r µb µg (1 + c3 H˜b ) ˜2H ˜b ¯2 aR eb aR ˜ H˜b ¯ R = ; S S0R = . ¯2) 1 ˜ b + δ(b + R ˜ 2 )2 µb (b + R abH Theorem 1 (Local stability results). The following results hold. (1) If r < 0, then E0 = 0, 0, 0 is locally asymptotically stable. r , 0, 0 is locally asymptotically stable. (2) If S0 < 1 and L0 > 1, then E1 = δ ˜ ˜ ˜ ˜ (3) Suppose that S0RHb < 1 and S1RHb < 1. We have two cases. ˜ ˜ ˜ (a) Grazer-free equilibrium E = R, Hb , 0 is locally asymptotically stable whenever it is unique. (b) If there exist two grazer-free equilibria, then one is locally asymptotically stable and the other is unstable. ¯ R ¯ (4) Suppose that S0 < 1, then we have two cases for the browser-free equilibria E = ¯ 0, H ¯g . R, (a) If we have one equilibrium, then it is locally asymptotically stable.
(b) If we have three equilibria, then two are locally asymptotically stable and the other is unstable. Proof: See Appendix B.3. Before going further, let us stress the fact that in Theorem 1 the case r = 0 is not considered. However, when r = 0, using a comparison argument (e.g. Yatat et al. 2014 [39]), solution of System (1) go to extinction. Furthermore, from Theorem 1, one deduces that multistability involving several equilibria of the same type or of different types may exist. This means that depending on initial conditions, the system may converge to different equilibria. From an ecological point of view, these multistabilities mean that transitions between herd, browsers and grazers amount may be more complex and confirm that these systems in Sahelian are very complex systems. In the sequel we provide some interpretations of Theorem 1 in terms of pastoral mechanisms. Pastoral interpretation
• The stability of equilibrium E0 = 0, 0, 0 means that there will be no more forage resources and no animals of the herd in the considered area. There are many explanations of this situation. Droughts can cause animals mortalities through starvation, emergency slaughtering, and sales, or definitive herds migrations (transhumance), which can create severe drops in the herd sizes. • When this removal concerns only herbivores, depending on some climatic changes, the forage resources growth towards the maximal quantity needed for livestock and the r equilibrium E1 = , 0, 0 is stable. δ 8
• If the disappearance concerns only grazers, we have a coexistence between the resource and the browsers. A situation that can be explained by the extinction of the grass resource and therefore the small herbivores may disappear or the pastor may decide to migrate with the grazers (transhumance). Consequently, the grazers-free equilibrium ˜ H ˜ b , 0 is stable. E˜ = R, • If the disappearance concerns only the browsers, we have a coexistence between the resource and the grazers. A situation that can be explained by the insufficiency of the available resource because browsers need a large availability of resources. Consequently, this situation can cause pastoralists to migrate with browsers and therefore, ¯ 0, H ¯ g is sable. the browsers-free equilibrium E¯ = R,
Theorem 2 (Hopf bifurcation). ˜ H ˜ b , 0) when δ The System (1) has a Hopf bifurcation around grazer-free equilibria E˜ = (R, passes through δh , where ˜ 2 − b) aH˜b (R . (4) δh = ˜ 2 )2 (b + R Proof: See Appendix B.4. The pastoral interpretation of the Hopf bifurcation is that herbivores will coexist with the resource, exhibiting oscillatory balance behavior. Theorem 3 (Local stability of coexistence equilibria). The coexistence equilibrium E ∗ = (R∗ , Hb∗ , Hg∗ ) is locally asymptotically stable if the conditions of the Routh-Hurwitz criterion are verified, that is: C1 < 0,
C2 < 0
and
C3 > 0,
(5)
where C1 = m∗11 + m∗33 , C2 = (m∗11 + m∗33 )(m∗11 m∗33 − m∗13 m∗31 − m∗12 m∗21 ) + (m∗12 m∗21 m∗33 − m∗13 m∗21 m∗32 ), C3 = m∗12 m∗21 m∗33 − m∗13 m∗21 m∗32 , with m∗11 , m∗33 , m∗12 , m∗21 , m∗13 and m∗31 as given in the proof. Proof: See Appendix B.5. 5. Stability analysis with positive delays In this section we study the stability behavior of equilibria in the presence of discrete delays τ1 > 0 and τ2 > 0. Let E ∗ = (R∗ , Hb∗ , Hg∗ ) be an equilibrium of the model, and let us define u(t) = R(t) − R∗ , v(t) = Hb (t) − Hb∗ and w(t) = Hg (t) − Hg∗ . Then the linearized system of (1) at E ∗ is given by
9
du(t) = a1 u(t) + a2 v(t) + a3 w(t); dt dv(t) (6) = m1 u(t − τ1 ) + b2 v(t); dt dw(t) = m u(t − τ ) + c v(t) + c w(t); 2 2 2 3 dt where αb Hb∗ + αg Hg∗ c1 Hg∗ αb Hb∗ + αg Hg∗ 2abR∗ 1− −2δR∗ − − a1 = r 1 + Hb∗ ; ∗ ∗ ∗2 2 β γ 1 + c2 Hg + c3 Hb (b + R ) ∗2 c1 c3 Hg∗ rαb aR a2 = γ − β − 2(αb Hb∗ + αg Hg∗ ) R∗ − R∗ ; + βγ b + R∗2 (1 + c2 Hg∗ + c3 Hb∗ )2 c1 R∗ (1 + c3 Hb∗ ) rαb 2eb abR∗ ∗ ∗ a3 = γ − β − 2(αb Hb∗ + αg Hg∗ ) R∗ − = ; m H ; 31 βγ (1 + c2 Hg∗ + c3 Hb∗ )2 (b + R∗2 )2 b eg c1 Hg∗ 2eb abR∗ ∗ eb aR2 m1 = H ; b2 = − µb ; m2 = ; (b + R∗2 )2 b b + R2 1 + c2 Hg∗ + c3 Hb∗ eg c1 c3 Hg∗ R∗ µg c2 Hg∗ c2 = − ; c3 = − . (1 + c2 Hg∗ + c3 Hb∗ )2 1 + c2 Hg∗ + c3 Hb∗ We then express system (6) in matrix form as follow u(t) u(t) u(t − τ1 ) u(t − τ2 ) d v(t) = A v(t) + A1 v(t − τ1 ) + A2 v(t − τ2 ) , dt w(t) w(t) w(t − τ1 ) w(t − τ2 ) where
a1 a2 a3 0 0 0 A = 0 b2 0 , A1 = m1 0 0 , 0 c2 c3 0 0 0
0 0 0 A2 = 0 0 0 . m2 0 0
The characteristic equation of system (6) at equilibrium E ∗ is given by ∆(ξ, τ1 , τ2 ) = |ξI − A − A1 e−ξτ1 − A2 e−ξτ2 | = 0, that is, P0 (ξ) + P1 (ξ)e−ξτ1 + P2 (ξ)e−ξτ2 = 0,
(7)
where, P0 (ξ) = ξ 3 − (a1 + c3 + b2 )ξ 2 + (a1 c3 + a1 b2 + c3 b2 )ξ − a1 c3 b2 , P1 (ξ) = −m1 a2 ξ + m1 a2 c3 − m1 a3 c2 and P2 (ξ) = −a3 m2 ξ + a3 m2 b2 . 5.1. Stability analysis when τ1 > 0 and τ2 = 0
0, 0, 0 ,
In the case where τ1 > 0 and τ2 = 0, the characteristic equations at E0 = r ¯ 0, H ¯ g are independent of the delay τ1 . Then equilibria E0 , E1 , 0, 0 and E¯ = R, E1 = δ and E¯ have the same stability as in the delays-free case.
10
˜ existence of a switch threshold 5.1.1. Stability analysis of E: ˜ When τ1 > 0 and τ2 = 0, the characteristic equation In the sequel, we only focus on E. ˜ at E is given by (8) ξ 3 + B1 ξ 2 + B2 ξ + (B3 ξ + B4 )e−ξτ1 = 0. where B1 = −(˜ a1 + c˜3 ), B2 = a ˜1 c˜3 , B3 = −m ˜ 1a ˜2 and B4 = a ˜2 c˜3 m ˜ 1. Substituting ξ = iω in the characteristic equation (8), the real and the imaginary part satisfy ω 3 − B2 ω = B3 ω cos ωτ1 − B4 sin ωτ1 ,
(9)
B1 ω 2 = B3 ω sin ωτ1 + B4 cos ωτ1 .
(10)
Eliminating τ1 by squaring and adding equations (9) and (10), we get the algebraic equation ω 6 + (B12 − 2B2 )ω 4 + (B22 − B32 )ω 2 − B42 = 0.
(11)
Substituting ω 2 = η in the above equation, we obtain the following cubic equation η 3 + θ20 η 2 + θ10 η + θ00 = 0,
(12)
where, θ20 = B12 − 2B2 , θ10 = B22 − B32 and θ00 = −B42 . Note that θ20 = B12 − 2B2 = (˜ a1 + c˜3 )2 − 2˜ a1 c˜3 = a ˜21 + c˜23 > 0 and θ00 = −B42 < 0. By Descartes’rule of sign, the equation (12) has at least one positive root. Remark 1. When τ1 = τ2 = 0, by Theorem 1 (3), E˜ is locally asymptotically stable if ˜ ˜ ˜ ˜ ˜1 < 0) and a ˜2 < 0. S0RHb < 1 (i.e. c˜3 < 0), S1RHb < 1 (i.e. a The following theorem gives a criterion for the switching in the stability behavior of E˜ in terms of the delay parameter τ1 . Theorem 4. Suppose that E˜ exists and is locally asymptotically stable for (1) with τ1 = τ2 = 0. Also let, ω02 = η0 be a positive root of (12). Then, there exists τ1 = τ10 such that E˜ is locally asymptotically stable for τ1 ∈ (0, τ10 ] and unstable for τ1 > τ10 . Furthermore, the system undergoes a Hopf bifurcation at E˜ when τ1 = τ10 . Proof: See Appendix B.6. ˜ direction and stability of the Hopf Bifurcation 5.1.2. Stability analysis of E: ˜ H ˜ b , 0) of system In the previous section, the stability of the positive equilibrium E˜ = (R, (1) has been analized as well as the existence of Hopf bifurcation at this equilibrium. In this section, we shall study the properties of the Hopf bifurcation obtained in Theorem 4 and the stability of bifurcated periodic solutions occurring through Hopf bifurcation by using the normal form theory and the center manifold reduction for retarded functional differential equations (RFDEs) due to Hassard et al. 1981 [16]. Throughout this section, we always assume that system (1) undergoes Hopf bifurcation at the positive equilibrium 0 ˜ H ˜ b , 0) for τ1 = τ1j E˜ = (R, , (j = 0, 1, 2....) and then ±iω0 is the corresponding purely imaginary roots of the characteristic equation (8). Theorem 5. There exist three real numbers µ2 , β2 and T2 which determine the qualities of bifurcating periodic solution in the center manifold at the critical value τ10 such that: 11
(a) When µ2 > 0, then the Hopf bifurcation is supercritical and the bifurcating periodic solution exists for τ1 > τ10 . When µ2 < 0, then the Hopf bifurcation is subcritical and the bifurcating periodic solution exists for τ1 < τ10 . (b) When β2 < 0, the bifurcating periodic solutions are stable and when β2 > 0, the bifurcating periodic solutions are unstable. (c) The period increases when T2 > 0, and decreases when T2 < 0. Proof: See Appendix B.7. 5.1.3. Stability analysis of E ∗ when τ1 > 0 and τ2 = 0 When τ1 > 0 and τ2 = 0, the characteristic equation at E ∗ is given by ∆(ξ, τ1 ) = M (ξ) + N (ξ)e−ξτ1 = 0,
(13)
where M (ξ) = ξ 3 + A2 ξ 2 + A1 ξ and N (ξ) = S1 ξ + S0 with A2 = −(a∗1 + c∗3 ), A1 = a∗1 c∗3 − a∗3 m∗2 , S1 = m∗1 a∗2 and S0 = m∗1 a∗2 c∗3 − m∗1 a∗3 c∗2 . Substituting ξ = iω in the characteristic equation (13), the real and imaginary parts satisfy ω 3 − A1 ω = S1 ω cos ωτ1 − S0 sin ωτ1 ,
(14)
A2 ω 2 = S1 ω sin ωτ1 + S0 cos ωτ1 .
(15)
Eliminating τ1 by squaring and adding equations (14) and (15), gives the algebraic equation ω 6 + (A22 − 2A1 )ω 4 + (A21 − S12 )ω 2 − S02 = 0.
(16)
Substituting ω 2 = η in the above equation, gives the following cubic equation η 3 + (A22 − 2A1 )η 2 + (A21 − S12 )η − S02 = 0.
(17)
By the intermediate value theorem, the equation (17) has at least one positive root. The following theorem gives a criterion for the switching in the stability behavior of E ∗ in terms of the delay parameter τ1 . Theorem 6. Suppose that E ∗ exists and is locally asymptotically stable for (1) with τ1 = τ2 = 0. Also suppose that ω12 = η1 is a positive root of (17). Then, there exists τ1 = τ100 such that E ∗ is locally asymptotically stable for τ1 ∈ (0, τ100 ] and unstable for τ1 > τ100 . Furthermore, the system undergoes a Hopf bifurcation at E ∗ when τ1 = τ100 . The proof of Theorem 6 is done in the same way as that of Theorem 4. 5.1.4. Stability analysis of E ∗ : direction and stability of the Hopf Bifurcation Similar to section 5.1.2, the direction of Hopf of E ∗ is given on the sequel Theorem 7. There exist three real numbers µ3 , β3 and T3 which determine the quality of the bifurcating periodic solution in the center manifold at the critical value τ100 such that: (a) When µ3 > 0, then the Hopf bifurcation is supercritical and the bifurcating periodic solution exists for τ1 > τ100 . When µ3 < 0, then the Hopf bifurcation is subcritical and the bifurcating periodic solution exists for τ1 < τ100 . (b) The bifurcating periodic solutions are stable when β3 < 0, and unstable if β3 > 0 (c) The period increases when T3 > 0 and decreases when T3 < 0. The proof of the Theorem 7 is done in the same way as that of Theorem 5. 12
5.2. Stability analysis when τ1 = 0 and τ2 > 0
0, 0, 0 ,
In the case where τ1 = 0 and τ2 > 0, the characteristic equations at E0 = r ˜ 0, H ˜ g are independent of the delay τ2 . Then, equilibria E0 , E1 E1 = , 0, 0 and E˜ = R, δ and E˜ have the same stability properties as in the delays-free cases. 5.2.1. Analysis of E¯ when τ1 = 0 and τ2 > 0 When τ1 = 0 and τ2 > 0, the characteristic equation at E¯ is given by ξ 3 + C1 ξ 2 + C2 ξ + C3 + (C4 ξ + C5 )e−ξτ2 = 0.
(18)
where C1 = −(¯ a1 + c¯3 + ¯b2 ); C2 = a ¯1 c¯3 + ¯b2 a ¯1 + ¯b2 c¯3 ; C3 = −¯ a1 c¯3¯b2 ; C4 = −¯ a3 m ¯ 2 and C5 = a ¯3 m ¯ 2¯b2 . Substituting ξ = iω in the characteristic equation (18), the real and imaginary parts satisfy ω 3 − C2 ω = C4 ω cos ωτ2 − C5 sin ωτ2 ,
(19)
C1 ω 2 − C3 = C4 ω sin ωτ2 + C5 cos ωτ2 .
(20)
Eliminating τ2 by squaring and adding equations (19) and (20), gives the algebraic equation ω 6 + (C12 − 2C2 )ω 4 + (C22 − C42 − 2C3 C1 )ω 2 + C32 − C52 = 0.
(21)
Substituting ω 2 = η in the above equation, gives the following cubic equation η 3 + θ2 η 2 + θ1 η + θ0 = 0,
(22)
where, θ2 = C12 − 2C2 ; θ1 = C22 − C42 − 2C3 C1 and θ0 = C32 − C52 . We have θ2 = C12 − 2C2 = (¯ a1 + c¯3 + ¯b2 )2 − 2(¯ a1 c¯3 + ¯b2 a ¯1 + ¯b2 c¯3 ) = a ¯21 + c¯23 + ¯b22 > 0 and θ0 = C32 − C52 = (C3 − C5 )(C3 + C5 ) = (−¯ a1 c¯3¯b2 − a ¯3 m ¯ 2¯b2 )(−¯ a1 c¯3¯b2 + a ¯3 m ¯ 2¯b2 ) = ¯b22 (¯ a1 c¯3 + a ¯3 m ¯ 2 )(¯ a1 c¯3 − a ¯3 m ¯ 2 ). We consider only the case when (22) possesses exactly one positive root. This occurs if ¯1 c¯3 + a ¯3 m ¯ 2 < 0. Under this condition, the equation (22) possesses a θ0 = C32 − C52 < 0 i.e. a 2 positive root η2 = ω2 and thus (18) has the pair of purely imaginary roots ±iω2 . Remark 2. When τ1 = τ2 = 0, by Theorem 1 (4) E¯ is locally asymptotically stable, if ¯ ¯3 < 1. S0R < 1 and a Using a similar approach as in Theorem 4, one deduces the following result. Theorem 8. Suppose that E¯ exists and is locally asymptotically stable for (1) when τ1 = τ2 = 0. Also suppose that, ω22 = η2 is a positive root of (22). Then, there exists a τ2 = τ20 such that E¯ is locally asymptotically stable for τ2 ∈ (0, τ20 ] and unstable for τ2 > τ20 . Furthermore, the system undergoes a Hopf bifurcation at E¯ when τ2 = τ20 . The proof of Theorem 8 is done in the same way as that of Theorem 4.
13
5.2.2. Nonexistence of switching stability around E¯ When τ1 = 0 and τ2 > 0, the characteristic equation at E¯ is given by ξ 3 + C1 ξ 2 + C2 ξ + C3 + (C4 ξ + C5 )e−ξτ2 = 0.
(23)
Equation (23) can be written as P (ξ) + Q(ξ)e−ξτ2 = 0,
(24)
where, P (ξ) = ξ 3 + C1 ξ 2 + C2 ξ + C3 and Q(ξ) = C4 ξ + C5 . Clearly P (ξ) and Q(ξ) are both analytic functions in Re (ξ) > 0. The following results hold: i) P (0) + Q(0) = C3 + C5 > 0, ¯ ii) P (−iy) = P¯ (iy) and Q(−iy) = Q(iy), Q(ξ) iii) lim sup = 0. |ξ|→+∞ P (ξ)
In addition
F (y) = |P (−iy)|2 − |Q(iy)|2 = y 6 + (C12 − 2C2 )y 4 + (C22 − C42 − 2C3 C1 )y 2 + C32 − C52 , is a cubic expression in y 2 . By using simple algebraic calculation we have, C12 − 2C2 = a ¯21 + c¯23 + ¯b22 > 0, C22 − C42 − 2C3 C1 = a ¯21¯b22 + ¯b22 c¯23 + (¯ a1 c¯3 + a ¯3 m ¯ 2 )(¯ a1 c¯3 − a ¯3 m ¯ 2 ), 2 2 2 C3 − C5 = ¯b2 (¯ a1 c¯3 + a ¯3 m ¯ 2 )(¯ a1 c¯3 − a ¯3 m ¯ 2 ), ¯1 c¯3 + a ¯3 m ¯ 2 > 0. C32 − C52 > 0 and C22 − C42 − 2C3 C1 > 0 if a Therefore F (y) = 0, has no positive real root if a ¯1 c¯3 +¯ a3 m ¯ 2 > 0. Hence by applying Theorem 4.1 in Kuang 1993 [18], the system (1) does not possess any stability switches around E¯ when a ¯1 c¯3 + a ¯3 m ¯ 2 > 0. Remark 3. Since the system (1) does not always undergo a Hopf bifurcation at the positive ¯ we are not interested by the direction of the Hopf bifurcation at equilibrium equilibrium E, E¯ 5.2.3. Stability analysis of E ∗ when τ1 = 0 and τ2 > 0 When τ1 = 0 and τ2 > 0, the characteristic equation at E ∗ is given by P0 (ξ) + Q0 (ξ)e−ξτ2 = 0,
(25)
where, P0 (ξ) = ξ 3 + R2 ξ 2 + R1 ξ + R0 and Q0 (ξ) = T ξ with R2 = −(a∗1 + b∗2 + c∗3 ), R1 = a∗1 b∗2 + c∗3 a∗1 + b∗2 c∗3 − a∗2 m∗1 , R0 = c∗3 a∗2 m∗1 − m∗1 a∗3 c∗2 and T = −a∗3 m∗2 . Clearly P0 (ξ) and Q0 (ξ) are both analytic functions in Re (ξ) > 0. The following results hold: i) P0 (0) + Q0 (0) = R0 6= 0, ¯ 0 (ix), ii) P0 (−ix) = P¯0 (ix) and Q0 (−ix) = Q 14
Q (ξ) 0 iii) lim sup = 0. |ξ|→+∞ P0 (ξ)
Then,
F (x) = |P0 (−ix)|2 − |Q0 (ix)|2 = x6 + (R22 − 2R1 )x4 + (R12 − T 2 − 2R0 R2 )x2 + R02 = 0, is a cubic expression in x2 . Therefore, by Descartes’rule of sign, the equation F (w) = 0 has no positive real root if R22 − 2R1 > 0 and R12 − T 2 − 2R0 R2 > 0.
(26)
By applying Theorem 4.1 in Kuang 1993 [18], when τ1 = 0, τ2 > 0 and conditions (26) hold, the system does not possess any stability switches around E ∗ . 5.3. Stability analysis when τ1 > 0 and τ2 > 0 Before going further, let us recall the following results. (a) When τ1 > 0 and τ2 > 0, the characteristic equation at E¯ can be written as P0 (ξ) + P1 (ξ) + P2 (ξ)e−ξτ2 = 0, where, P0 , P1 and P2 are some polynomials. Therefore, the delay τ1 has no influence ¯ on equilibrium E. (b) When τ1 > 0 and τ2 > 0, the characteristic equation at E˜ can be written as P0 (ξ) + P2 (ξ) + P1 (ξ)e−ξτ1 = 0, where, P0 , P1 and P2 are some polynomials. Therefore, the delay τ2 has no influence ˜ on equilibrium E. In what follows, we use the algebraic method developed in Lin et al. 2012 [21] to derive a closed form for stability switching curves of our delayed system when τ1 > 0 and τ2 > 0. We provide some properties of these curves and stability switching direction. To applied this method, let us write the following characteristic equation at equilibrium E ∗ , D(ξ, τ1 , τ2 ) = P0 (ξ) + P1 (ξ)e−ξτ1 + P2 (ξ)e−ξτ2 = 0,
(27)
where, P0 = ξ 3 + R1 ξ 2 + R2 ξ, P1 = S1 ξ + S2 and P2 = T1 ξ, with R1 , R2 , S1 , S2 and T1 given by equation (7). A- Stability crossing curves Lemma 2. As (τ1 , τ2 ) varies continuously in R2+ , the number of zeros (with multiplicity counted) of D(ξ, τ1 , τ2 ) on C+ can change only if a zero appears on or cross the imaginary axis. The proof of Lemma 2 can be found in any book on functional differential equations, for example in Kuang 1993 [18] and Smith 2010 [34]. Remark 4. From the above lemma, to study stability switch, we need to find the roots on the imaginary axis. 15
In the following, we investigate the existence of purely imaginary roots ξ = iω (w > 0) of equation (27). Substituting this into (27), gives −iωτ2 + P1 (iω)e−iωτ1 = 0. D(iω, τ1 , τ2 ) = P0 (iω) + P2 (iω)e Since |e−iωτ1 | = 1, we have
|P0 (iω) + P2 (iω)e−iωτ2 | = |P1 (iω)|, which is equivalent to
P 0 + P2 e
−iωτ2
(28)
iωτ2 ¯ ¯ = |P1 |2 P0 + P2 e
where P0 = −R1 ω 2 + i(−ω 3 + R2 ω), P1 = iS1 ω + S2 and P2 = iT1 ω. After simplification, (28) becomes |P0 |2 + |P2 |2 + 2Re (P0 P¯2 ) cos(ωτ2 ) − 2Im (P0 P¯2 ) sin(ωτ2 ) = |P1 |2 . Thus,
|P1 |2 − |P0 |2 − |P2 |2 = 2Re (P0 P¯2 ) cos(ωτ2 ) − 2Im (P0 P¯2 ) sin(ωτ2 ). (29) Since P0 6= 0 and P2 6= 0, we have [Re (P0 P¯2 )]2 + [Im (P0 P¯2 )]2 = |P0 P¯2 |2 > 0. Therefore, equation (29) can be written as |P1 |2 − |P0 |2 − |P2 |2 Re (P0 P¯2 ) Im (P0 P¯2 ) = cos(ωτ ) − sin(ωτ2 ). 2 2|P0 P¯2 | |P0 P¯2 | |P0 P¯2 |
(30)
Then, there exists a continuous function φ2 (ω) such that cos(φ2 (ω)) =
Re (P0 P¯2 ) Im (P0 P¯2 ) and sin(φ (ω)) = . 2 |P0 P¯2 | |P0 P¯2 |
Therefore, equation (30) becomes |P1 |2 − |P0 |2 − |P2 |2 = cos(φ2 (ω) + ωτ2 ). 2|P0 P¯2 |
(31)
Obviously, a sufficient and necessary condition for the existence of τ2 ∈ R+ satisfying equation (31) is 2 2 2 (32) |P1 | − |P0 | − |P2 | ≤ 2|P0 P¯2 |. Let us consider
n o Ω2 = ω ∈ R+ : |P1 |2 − |P0 |2 − |P2 |2 ≤ 2|P0 P¯2 | .
We have the following result.
Lemma 3. The set Ω2 is non empty.
16
Proof: See Appendix B.8. Let |P1 |2 − |P0 |2 − |P2 |2 ˜ , ψ2 ∈ [0, π]. 2|P0 P¯2 |
cos(ψ˜2 (ω)) = We obtain
ψ˜2 (ω) − φ2 (ω) + 2πn2 , n2 ∈ Z. (33) ω Substituting (33) into (27), gives the expression of τ1 (ω) unconditionally with each ω ∈ Ω2 i.e. P1 1 + 2πn1 , n1 ∈ Z. (34) τ1,n1 (ω) = arg − ω P0 + P2 e−iωτ2 Thus, the stability crossing curves are I = (τ1 (ω), τ2 (ω)) ∈ R2+ , ω ∈ Ω2 . (35) τ2,n2 (ω) =
Another way to find τ1 is to analyze τ1 similarly to the analysis of τ2 , which gives τ1,n1 (ω) =
ψ˜1 (ω) − φ1 (ω) + 2πn1 , n1 ∈ Z, ω
(36)
where, cos(ψ1 (ω)) =
|P2 |2 − |P0 |2 − |P1 |2 , ψ1 ∈ [0, π], 2|P0 P¯1 |
with the following condition on ω: 2 2 2 |P2 | − |P0 | − |P1 | ≤ 2|P0 P¯1 |. Let us consider
(37)
n o Ω1 = ω ∈ R+ : |P2 |2 − |P0 |2 − |P1 |2 ≤ 2|P0 P¯1 | .
By squaring both sides of the two conditions (32) and (37), one can show that (32) and (37) are equivalent. Thus, Ω1 = Ω2 ≡ Ω. We call Ω the crossing set.
Lemma 4. The set Ω consists of a finite number of intervals of finite length. Proof: See Appendix B.9. Remark 5. Since F (0) > 0, then F (ω) has roots 0 < a1 < b1 ≤ a2 < b2 < ... ≤ aN < bN < ∞, (N ≤ 6) and Ω=
N [
Ωk , where Ωk = [ak , bk ].
k=1
We allow bk−1 = ak with concern that ak may be a root with even multiplicity. 17
It is obvious that the stability crossing curves obtained in this way should be the same as I in (35). Since τ1 , τ2 obtained in both ways depend only on ω ∈ Ω, they should be the same respectively and can be defined as in (33) and (36) as well to get the crossing curves, i.e., n o I= τ1,n1 (ω), τ2,n2 (ω) ∈ R2+ , ω ∈ Ω, n1 , n2 = −1, 0, 1, 2... . Since Ω =
N [
Ωk , we have
k=1
I=
[
k=1,2,...,N
Ink1 ,n2 ,
(38)
n1 =−1,0,1,2... n2 =−1,0,1,2...
where Ink1 ,n2
n o 2 = τ1 (ω), τ2 (ω) ∈ R+ , ω ∈ Ωk , n1 , n2 = −1, 0, 1, 2... .
(39)
One can show that each Ink1 ,n2 is continuous in R2+ . In addition Ink1 ,n2 may connect with each other. Similar analysis can be found in Gu et al. 2005 [14]. B- Crossing direction Next, we will discuss the direction in which the solutions of the characteristic equation given by D(ξ, τ1 , τ2 ) = 0 cross the imaginary axis as (τ1 , τ2 ) deviate from a curve in Ink1 ,n2 . We will call the direction of the curve that corresponds to increasing ω, the positive direction. Observe that, the curve passes through the points corresponding to the end points of Ωk , the positive direction is reversed. For the purpose of discussing the direction of crossing, we need to consider τ1 and τ2 as functions of ξ = σ + iω, i.e., functions of two real variables σ and ω, and partial derivative notation needs to be adopted instead. Since the tangent of Ink1 ,n2 along the ∂τ ∂τ 2 1 , , the normal to Ink1 ,n2 pointing to the right region of the positive direction is ∂ω ∂ω ∂τ ∂τ1 2 positive direction is ,− . Also, as a pair of complex conjugate solutions of ∂ω ∂ω D(ξ, τ1 , τ2 ) = 0 cross the imaginary axis to the C+ , (τ1 , τ2 ) moves along the direction ∂τ ∂τ 1 2 , . We can therefore conclude that the inner product of these two vectors are ∂σ ∂σ positive, i.e., ∂τ ∂τ ∂τ ∂τ1 1 2 2 Γ(ω) = , . ,− = det ∆(ω) > 0, ∂σ ∂σ ∂ω ∂ω where ∂τ1 ∂τ1 ∂σ ∂ω ∆(ω) = ∂τ . ∂τ2 2 ∂σ ∂ω σ=0,ω∈Ω
The region at the right of Ink1 ,n2 has two more characteristic roots with positive real parts. On the other hand, if the previous inequality is reversed, then the region at the left of Ink1 ,n2 has two more characteristic roots with positive real parts. By the implicit function theory, we have ˆ1 R ˆ 2 −1 R0 −I0 R ∆(ω) = , I0 R0 Iˆ1 Iˆ2 18
where, R0 = I0 = ˆ1 = R ˆ2 = R Iˆ1 = Iˆ2 = Since det
∂Re D(ξ, τ1 , τ2 ) ∂σ ξ=iω ∂Im D(ξ, τ1 , τ2 ) ∂σ ξ=iω ∂Re D(ξ, τ1 , τ2 ) ∂τ1 ξ=iω ∂Re D(ξ, τ1 , τ2 ) ∂τ2 ξ=iω ∂Im D(ξ, τ1 , τ2 ) ∂τ1 ξ=iω ∂Im D(ξ, τ1 , τ2 ) ∂τ2 ξ=iω R0 −I0 I0 R0
) 2 X , = Re P00 (iω) + Pk0 (iω) − τk Pk (iω)e−iωτk k=1 ) ( 2 X , = Im P00 (iω) + Pk0 (iω) − τk Pk (iω)e−iωτk (
k=1
−iωτ1
= ω 2 S1 cos ωτ1 − ωS2 sin ωτ1 , = Re − iωP1 (iω)e −iωτ2 = ω 2 T1 cos ωτ2 , = Re − iωP2 (iω)e −iωτ1 = −ω 2 S1 sin ωτ1 − ωS2 cos ωτ1 , = Im − iωP1 (iω)e = Im − iωP2 (iω)e−iωτ2 = −ω 2 T1 sin ωτ2 .
ˆ 1 Iˆ2 − R ˆ 2 Iˆ1 ). = R02 + I02 > 0, we have sign(Γ(ω)) = sign(R
ˆ 1 Iˆ2 − R ˆ 2 Iˆ1 = −ω 2 T1 sin ωτ2 ω 2 S1 cos ωτ1 − ωS2 sin ωτ1 R 2 2 + ω T1 cos ωτ2 ω S1 sin ωτ1 + ωS2 cos ωτ1 , h i = ω 2 T1 ωS2 cos(ω(τ1 − τ2 )) + ω 2 S1 sin(ω(τ1 − τ2 )) , 1 −τ2 ) = ω 2 Im P¯1 P2 eiω(τ , o n = −ω 2 |P0 P¯1 |Im eiφ1 +iωτ1 , = −ω 2 |P0 P¯1 | sin φ1 + ωτ1 .
Hence, we have the following result.
Theorem 9. For any k = 1, 2, ..., N , (i) If sin φ1 + ωτ1 < 0, then the region at the right of Ink1 ,n2 has two more characteristic roots with positive real parts. (ii) If sin φ1 + ωτ1 > 0, then the region at the left of Ink1 ,n2 has two more characteristic roots with positive real parts. Remark 6. If we know the number of characteristic roots with positive real parts when τ1 = τ2 = 0, we can use the criterion in Theorem 9 to find the number of characteristic roots with positive real parts for any (τ1 , τ2 ) in R2+ . Hence, the stability of the characteristic equation is completely known. 6. Numerical Simulations Using a nonstandard numerical scheme as in Anguelov et al. 2012 [2], we provide the following numerical simulations where parameter values have been chosen in such a way that they obey the conditions for stability or bifurcation. 19
6.1. Numerical Simulations when τ1 = τ2 = 0 In this section, we present some numerical simulations to illustrate the analytical results observed in the previous sections. We first consider the following set of parameter values: r = 0.01, β = 0.001, γ = 2, αb = 0.1, αg = 0.15, δ = 0.73, a = 0.6, c1 = 0.32, c2 = 0.3, c3 = 0.15, eb = 0.45, db = 0.015, hb = 0.15, eg = 0.65, dg = 0.75, hg = 0.8. Figure 1 depicts monostability as well as two types of bistabilities situations when considering the browser’s half saturation parameter b as a bifurcation parameter. Namely, monostability of the equilibrium E1 = (K, 0, 0), bistability between limit cycle and equilibrium E1 for b ∈ [0.2, 0.9] and the bistability between equilibrium E1 = (K, 0, 0) and equilibrium E˜ = (R, Hb , 0) for b ∈ [0.9, 1.5]. Specifically, in Figure 2, we illustrate these two bistabilities situation for two particular values of the browser’s half saturation parameter b. L0 >1, S0<1
14
12
Browsers (Hb)
10
Stable equilibrium (R, 0, 0)
(1)
Stable equilibrium (R, Hb , 0)
(1)
Unstable equilibrium (R, Hb , 0)
8
Stable limit cycle in the 2D R−Hb plan
6
(2)
Unstable equilibrium (R, Hb , 0)
4
2
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
Bifurcation parameter: b
Figure 1: A bifurcation diagram of system (1) without delays and when considering the browser’s half saturation parameter b as a bifurcation parameter.
20
Bistability of E1=(R, 0, 0) and E2=(R, Hb, 0) with b=1.4
15
• Browsers
0.1
Browsers
10
0.08 0.06 0.04 0.02
Zoom
0 0
•0.02
0.04 0.06 0.08
0.1
Resource
5
0 0
5
10
15
Resource
Bistability of E =(R, 0, 0) and the limit cycle with b=0.4 1
15
Browser
0.8
Browsers
10
0.6 0.4 0.2 0
5
0 0
0
•
0.05
0.1
Zoom
1
0.15
0.2
0.25
0.3
0.35
Resource
2
3
4
5
6
Resource
Figure 2: Phase portrait of resource and browsers without delays. Illustration of bistability between equi˜ = (R, Hb , 0) (top panel). Phase portrait of resource and browsers librium E1 = (K, 0, 0) and equilibrium E without delays. Illustration of bistability between equilibrium E1 = (K, 0, 0) and the limit cycle (bottom panel).
Let us now assume that r = 1.4, β = 0.1, γ = 2, αb = 0.1, αg = 0.015, δ = 0.33, a = 0.1, b = 0.2, c1 = 0.9, c2 = 0.03, c3 = 0.15, eb = 0.15, db = 0.015, hb = 0.15, eg = 0.3, dg = 0.3, hg = 0.4. In Figure 3 we illustrate a (forward) bifurcation phenomenon involving equilibria E¯ and E1 when we consider the grazer’s catch rate c1 as a bifurcation parameter. Indeed, for L0 > 1 i.e. c1 < 0.55, equilibrium E1 is stable while for L0 < 1 i.e. c1 > 0.55, equilibrium E1 becomes unstable and equilibrium E¯ is stable (see Figure 3-(i)).
21
(i)
(ii)
Forward bifurcation in system (1)
Phase portrait in the Resource−Grazers 2D plan
0.8
4.5 4
0.7
Stable equilibrium (R, 0, Hg)
Stable equilibrium (R, 0, 0)
Unstable equilibrium (R, 0, 0)
3
0.5
Grazers
g
Grazers (H )
0.6
3.5
0.4 0.3
L0>1
2.5 2 1.5
L0<1
0.2
1
0.1 0 0
•
0.5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0 0
1
1
2
3
4
5
6
7
8
9
Resource
Bifurcation parameter: c1
(iii)
(iv)
Phase portrait in the Resource−Browsers 2D plan
Phase portrait in the Browsers−grazers 2D plan
2
6
1.8 5
1.6 1.4 Grazers
Browsers
4 1.2 1 0.8
2
0.6 0.4
1
•
0.2 0 0
3
1
•
2
3
4
5
6
7
8
0 0
9
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Browsers
Resource
Figure 3: (i) Illustration, without delays, of the stability of equilibrium E1 = (K, 0, 0) for the bifurcation ¯ = (R, 0, Hg ) for c1 > 0.55. (ii) Phase portrait of resource parameter c1 < 0.55 and stability of equilibrium E ¯ = (R, 0, Hg ). (iii) Phase portrait of and grazers without delays. Illustration of stability of equilibrium E ¯ = (R, 0, Hg ). (iv) Phase resource and browsers without delays. Illustration of stability of equilibrium E ¯ = (R, 0, Hg ). portrait of browsers and grazers without delays. Illustration of stability of equilibrium E
6.2. Numerical Simulations with positive delays In what follows, the numerical experiments are performed on system (1) to illustrate our theoretical findings when delays are considered. From the above discussion when τ1 > 0 and τ2 = 0, we may determine the direction of the Hopf bifurcation and the stability of the bifurcating periodic solution. We consider the following set of parameter values: r = 0.01, β = 0.001, γ = 2, αb = 0.1, αg = 0.015, δ = 0.728, a = 0.6, b = 1.2, c1 = 0.322, c2 = 0.3, c3 = 0.15, eb = 0.45, eg = 0.65, db = 0.015, dg = 0.75, hb = 0.15, hg = 0.8. For these parameter values, conditions textcolorredin Theorem 4, page 11 are satisfied and system (1) becomes conditionally stable around the boundary equilibrium point E˜ for τ1 ∈ [0, τ10 ) (see, Figure 4) and unstable for τ1 > τ10 (e.g. see Figure 5). Indeed, when τ1 passes through d the critical value τ1 = τ10 = 0.3587 and (Reξ(τ1 ))|τ1 =τ10 = 0.1570 > 0, the equilibrium dτ1 E˜ losses its stability and the system (1) experiences a Hopf bifurcation. From Sections 5 and 5.1.2, we can determine the nature of the stability and direction of the periodic solution bifurcating from the boundary equilibrium E˜ at the critical point τ10 . We can compute C2 (0) = 4.2289 + 1.5425, µ2 = −26.9288 < 0, β2 = 8.4578 > 0 and T2 = −31.5892 < 0. Hence, the bifurcating periodic solution exists and the corresponding periodic solution is subcritical as illustrate in Figure 6. The period of the periodic solution is 31.5892 and its period decreases with τ1 since T2 < 0. 22
Figure 4: Dynamics of system (1) when τ1 = 0.1793 < τ10 .
Figure 5: Dynamics of system (1) when τ1 = 1.4348 > τ10 .
˜ = (1.3732, 13.1438, 0) of system Figure 6: Hopf-bifurcation behavior of system (1) around the equilibrium E ˜ (1) at τ1 = τ10 = 0.3587 and existence of subcritical bifurcating periodic solution around E.
From the above discussion when τ1 = 0 and τ2 > 0, we may determine the direction of Hopf-bifurcation and the direction of bifurcating periodic solution. We assume that: r = 4.1, β = 0.2, γ = 4.01, αb = 0.2, αg = 0.015, δ = 0.328, a = 0.1, b = 0.2, c1 = 0.9, c2 = 0.03, c3 = 0.15, eb = 0.15, eg = 0.3, db = 0.15, dg = 0.3, hb = 0.5, hg = 0.4. For these parameter values, conditions in Theorem 8 are satisfied and system (1) becomes 23
conditionally stable around the boundary equilibrium point E¯ for τ2 ∈ [0, τ20 ) (see, Figure 7) and unstable for τ2 > τ20 (see, Figure 8). d When τ2 passes through the critical value τ2 = τ20 = 0.8183 and (Reξ(τ2 ))|τ2 =τ20 = dτ2 0.0225 > 0, the equilibrium E¯ losses its stability and the system (1) experiences Hopf bifurcation.
Figure 7: Dynamics of system (1) when τ2 = 0.40915 < τ20 .
Figure 8: Dynamics of system (1) when τ2 = 1.6366 > τ20 .
24
¯ = (18.1235, 0, 199.6832) Figure 9: Hopf-bifurcation behavior of system (1) around the interior equilibrium E 0 of system (1) at τ2 = τ2 = 0.8183 and existence of subcritical bifurcating periodic solution around the ¯ interior equilibrium E.
7. Conclusion-discussion In this paper, we construct and analyze a delayed competition model between herbivores on a same resource in a given location in a Sahelian region. We consider three state variables namely, the resource biomass (R), the grazer biomass (Hg ) and the browser biomass (Hb ). As modelling option we assume that herbivory can stimulate resource biomass production. Consequently the herbivory optimization hypothesis (HOH) is considered for the resource growing term. Recall that the HOH was also used in others herbivory systems (e.g. Lebon et al. 2014 [19]). We also assume a Holling type III functional response for browsers while a Beddington-DeAngelis like-functional response is considered for grazers. Moreover, we also acknowledged two discrete delays as to account for times necessary to convert resource biomass into herbivore (grazers and browsers) biomasses. Finally, since for the Sahelian region drought can lead a low biomass production, we also consider the case where the resource intrinsic growth rate can be negative. To authors knowledge, it is the first time that a mathematical model designed to model resource-grazers-browsers dynamics in Sahelian region is presented and studied. We first analyzed the delays-free case and secondly study the case where delays are considered. In both cases, the model exhibited several equilibria including extinction equilibrium (both resource and herbivores biomasses are null), herbivores-free equilibrium, grazersfree equilibria, browsers-free equilibria and coexistence equilibria. Specifically, we found that, thanks to the HOH and the Holling type III functional response for browsers, up to two grazers-free equilibria may exist. Similarly, thanks to the HOH and the BeddingtonDeAngelis functional response for grazers, up to three browsers-free equilibria may also exist. Concerning the coexistence equilibria, we found that up to two coexistence equilibria may exist (see Proposition 1, page 6). We further carry out stability analysis of various equilibria in the delays-free case. Based on our mathematical analysis, we found that both resource and herbivores may go to extinction as a consequence of drought for example (i.e. the resource intrinsic growth rate r is such that r < 0). This phenomenon can motivate pastor to decide to migrate with the cattle 25
(transhumance). We also found several thresholds that shape stability of equilibria. Notably, we found that several multistability situations may exist involving aforementioned equilibria (see Theorem 1, page 8). Recall that a multistability situation occurs when several equilibria are simultaneously stable. In addition, we found that when varying some parameters, several bifurcations occurs including Hopf bifurcation (e.g. Theorem 4, page 9) and forward bifurcation. Indeed, with assumed parameter values, we obtained two bifurcation diagrams with the browser’s half saturation parameter b and the grazer’s catch rate c1 as bifurcation parameters respectively (see Figures 1 and 3). Notably, Figure 1 depicts monostability of herbivores-free equilibrium as well as two types of bistability situations when considering the browser’s half saturation parameter b as a bifurcation parameter. Specifically, we found bistability between the herbivores-free equilibrium and a limit cycle when b ∈ [0.2; 0.9]. When b ∈ [0.9; 1.5], we obtain a bistability situation involving the herbivores-free and the grazer-free equilibria. In Sahelian context, Hopf bifurcation denotes that herbivores and resource biomass exhibit periodic behavior. Similarly when the grazer-free (resp. browser-free) equilibrium is stable, its may be the consequence of migration of grazer (reps. browser) by pastor to another location. When discrete delays are considered, we obtain sufficient conditions on parameters for which the delay-induced system is asymptotically stable around an equilibrium (that was also stable in the delays-free case) for all values of the delay parameter. When the conditions are not satisfied, then there exists a critical value of the delay parameter below which the system is stable and above which the system is unstable. By applying the normal form theory and the center manifold theorem, the explicit formula which determine the stability and direction of the bifurcating periodic solutions are computed. Our analytical results show that when τ1 and τ2 pass through the critical values τ10 and τ20 , the equilibria E˜ and E¯ loss their stability and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate ¯ We further illustrate these results through numerical simulations (see Figures from E˜ and E. 4-9). A limitation of this study relies on the non-availability of parameter values that correspond to Sahelian regions. In this study we assume sets of parameter values in such a way that they fulfilled stability or bifurcation conditions. However, it is desirable that some field experiments may be carried out in order to assess at least the magnitude of (some) parameters. For example, assess the magnitude or the range of the browser’s half saturation parameter b and the grazer’s catch rate c1 who are two bifurcation parameters of the model. This study is a first step of a more general study that aims to provide through mathematical modelling some insights for the browsers, grazers and their resource interactions in Sahelian region. This study in a given location in Sahelian can be improved in several ways. First, one can take into account both temporal and spatial behaviors of browsers, grazers and their resource interactions. It can be done through only local spatial operator such as Laplace operator or with both local and non-local spatial behaviors with some kernels. Another improvement may consist on taking into account several spatial locations and look for the best strategy in order to reach both sustainable resource management in the different locations and a good herbivores nutrition. This last point will be the aim of a forthcoming paper. Acknowledgements We would like to thank the referees for their helpful comments and valuable suggestions. 26
References [1] U. African. African Union Policy Framework for Pastoralism in Africa: Securing, Protecting and Improving the Lives, Livelihoods and Rights of Pastoralist Communities. African Union, 2010. [2] R. Anguelov, Y. Dumont, and J.M. Lubuma. On nonstandard finite difference schemes in biosciences. AIP Conf. Proc., 1487:212–223, 2012. [3] M.L. Bates, R.A. Cropp, D.W. Hawker, and J. Norbury. Which functional responses preclude extinction in ecological population-dynamic models? Ecological Complexity, 26:57–67, 2016. [4] J. Bell. A grazing ecosystem in the serengeti. Sci Am., 225:86–89, 1971. [5] R.S. Cantrell and C. Cosner. On the dynamics of predator-prey models with the beddington-deangelis functional response. Journal of Mathematical Analysis and Applications, 257:206–222, 2001. [6] M. Carriere. Impact des syst`emes d’´elevage pastoraux sur l’environnement en afrique et en asie tropicale et sub-tropicale aride et sub-aride. livestock and the environment finding a balance. Scientific Environmental Monitoring Group, 1996. [7] G. Caughley. Plant-herbivore systems. Theoretical Ecology, Blackwell Scientific Publications, 320:530, 1976. [8] J.V. De Koppel and H.H.T. Prins. The importance of herbivore interactions for the dynamics of african savanna woodlands: an hypothesis. Journal of Tropical Ecology, 14:565–576, 1998. [9] D.L. DeAngelis, R.A. Goldstein, and R.V. O’Neill. A model for tropic interaction. Ecology, 56:881–892, 1975. [10] B. Dedieu, V. Ancey, P. Bommel, S. Cournut, L. Dobremez, B. Faye, A. Gibon, E. Josien, A. Ickowicz, B. Lemery, C. Macombe, S. Madelrieux, C.R. Poccard, H. Rapey, G. Serviere, and J.F.d Tourran. Transformations de l’´elevage et dynamiques des espaces. Projet TRANS, Rapport scientifique final, 2009. [11] Z. Feng, Z. Qiu, R. Liu, and D.L. DeAngelis. Dynamics of a plant herbivore-predator system with plant-toxicity. Mathematical Biosciences, 229:190–204, 2011. [12] H.I. Freedman and R.V.S. Hari. The trade-off between mutual interference and time lags in predator-prey systems. Bull. Math. Biol., 45:991–1004, 1983. [13] S. Gratton. Analyse matricielle et optimisation, edition ressources pedagogiques. Ouv. INPT, 0727, 2014. [14] K. Gu, S. Niculescu, and J. Chen. On stability crossing curves for general systems with two delays. J. Math. Anal. Appl., 311:231–252, 2005. [15] J.K. Hale. Asymptotic behavior of dissipative systems. American Mathematical Society, Providence, 1988. 27
[16] B.D. Hassard, N.D. Kazarinoff, and Y.H. Wan. Theory and applications of hopf bifurcation. Cambridge University, 1981. [17] D.W. Hilbert, D.M. Swift, J.K. Detling, and M.I. Dyer. Relative growth rates and the grazing optimization hypothesis. Oecologia, 51:14–18, 1981. [18] Y. Kuang. Delay differential equation with application in population dynamics. Academic Press New York, 1993. [19] A. Lebon, L. Mailleret, Y. Dumont, and F. Grognard. Direct and apparent compensation in plant-herbivore interactions. Ecological Modelling, 290:192–203, 2014. [20] Y. Li, Z. Feng, R. Swihart, J. Bryant, and N. Huntly. Modeling the impact of plant toxicity on plant-herbivore dynamics. Journal of Dynamics and Differential Equations, 18(4), 2006. [21] X. Lin and H. Wang. Stability analysis of delay differential equation with two discrete delays. Canadian Applied Mathematics Quarterly, 20, 2012. [22] A. Martin and S. Ruan. Predator-prey models with delay and prey harvesting. J. Math. Biol., 43:247–267, 2001. [23] S.J. McNaughton. Grazing lawns: animals in herds, plants form and coevolution. The American Naturalist, 124(6):863–886, 1984. [24] N. Nori, M. Taylor, and A. Sensi. Droits pastoraux, modes de vie et adaptation au changement climatique. Irish Aid, 2008. [25] I. Noy-Meir. Stability of grazing systems: An application of the predator-prey graphs. Journal of Ecology, 63:459–481, 1975. [26] Rapport. L’´elevage au sahel et en afrique de l’ouest. 26`eme runion annuelle du r´eseau de pr´evention des crises alimentaires (rpca). Accra (Ghana), 2010. [27] Rapport. Atlas des ´evolutions des syst`emes pastoraux au sahel, atlas of trends in pastoral systems in sahel. FAO et CIRAD, 2012. [28] C.V. Rebecca and S. Ruan. A delay-differential equation model of hiv infection of cd4 t-cells. Mathematical Biosciences (Elsevier), 165:27–39, 2000. [29] Research reports and studies. Humanitarian policy group, pastoralism policies and practice in the horn and east africa a review of current trends. 2009. [30] L.A. Shipley. Grazers and browsers: How digestived morphology affects diet selection. Journal of Tropical Ecology, 14:565–576, 1999. [31] A.R.E. Sinclair. Does interspecific competition or predation shape the african ungulate community? Anim Ecol, 54:899–918, 1985. [32] A.R.E Sinclair and M. Northon-Griffiths. Does competition or facilitation regulate migrant ungulate populations in the serengeti? a test of hypotheses. Oecologia, 53:364– 369, 1982. 28
[33] A.B. Smith. Origins and spread of pastoralism in africa. Nomadic Peoples, 32:91–105, 1993. [34] H. Smith. An introduction to delay differential equations with applications to the life sciences. Springer, New York, 2010. [35] J.J. Tewa, A. Bah, and S.C.O. Noutchie. Dynamical models of interactions between herds forage and water resources in sahelian region. Hindawi Publishing Corporation Abstract and Applied Analysis, 2014(13), 2014. [36] J.J. Tewa, S. Bowong, and N.S.C. Oukouomi. Mathematical analysis of a two–patch model of tuberculosis disease with staged progression. Applied Mathematical Modelling, 36(12):5792–5807, 2012. [37] J.J. Tewa, V.Y. Djeumen, and S. Bowong. Predator-prey model with holling response function of type ii and sis infectious disease. Applied Mathematical Modelling, 37(7):4825–4841, 2013. [38] S.C. Williamson, J.K. Detling, J.L. Dodd, and M.I. Dyer. Experimental evaluation of the grazing optimization hypothesis. Journal of Range Management ptimum gmrin9 Intensity, 42, 1989. [39] V. Yatat, Y. Dumont, J.J. Tewa, P. Couteron, and S. Bowong. Mathematical analysis of a size structured tree-grass competition model for savanna ecosystems. Biomath, (18), 2014.
29
Appendix A. Some results for the browser-free case Consider the following System when Hb = 0 i.e. dR α H α H c1 Hg g g g g = r 1+ 1− R − δR2 − R, dt β γ 1 + c2 Hg and let
dHg c1 Hg = eg R − (dg + hg )Hg , dt 1 + c2 Hg
αg Hg αg Hg c1 Hg r(1 + )(1 − )− β γ 1 + c2 Hg eg c1 RHg , f2 (R, Hg ) = µg (1 + c2 Hg ) F (R, Hg ) = (f1 , f2 )(R, Hg ). 1 f1 (R, Hg ) = δ
An equilibrium of system (A.1) satisfies dR(t) = 0 dt
(A.1)
,
f1 (R, Hg ) = R,
⇔ dH (t) f2 (R, Hg ) = Hg . g = 0 dt r , 0 as an equilibrium of (A.1). The computations give z0 = δ r Lemma A1. If L0 > 1, then ρ(DF (z0 )) < 1 where z0 = , 0 and ρ(A) denotes the δ spectral radius of matrix A. Proof: Let us consider the jacobian matrix at z0 . ∂f1 ∂f1 rαg (γ − β) 0 ∂R ∂Hg δβγ DF (z0 ) = = e g c1 r ∂f2 ∂f2 0 δµg ∂R ∂Hg (z=z0 ) Then ρ(DF (z0 )) =
.
eg c1 r 1 = . Since L0 > 1, one has δµg L0 ρ(DF (z0 )) < 1.
Lemma A2. If ρ(DF (z0 )) < 1, then there exists a matrix norm such that (see Gratton 2014 [13]): kDF (z0 )k < 1. 30
Proof: The following inequality holds ∀ε > 0, kDF (z0 )k ≤ ρ(DF (z0 )) + ε. For ε =
1 − ρ(DF (z0 )) , we have 2
kDF (z0 )k < 1.
Lemma A3. If kDF (z0 )k < 1, then F is strictly contracting. Proof: By continuity of DF at z0 we have: ∀ε > 0, ∃ρ > 0/ z ∈ B(z0 , ρ) ⇒ kDF (z)k < k(DF (z0 ))k + ε; ∀ε > 0, ∃ρ > 0/
sup z∈B(z0 ,ρ)
∀ε > 0, ∃ρ > 0/ kDF k = For ε =
kDF (z)k < kDF (z0 )k + ε;
sup z∈B(z0 ,ρ)
k(DF (z))k < kDF (z0 )k + ε.
1 − kDF (z0 )k , we have 2 ∃ρ > 0/
sup z∈B(z0 ,ρ)
kDF (z)k <
1 + kDF (z0 )k < 1, 2
i.e. kDF (z0 )k < 1. Lemma A4. If L0 > 1, then F is strictly contracting, and System (A.1) has a unique r equilibrium z0 = ,0 . δ
Proof: Since L0 > 1, Lemma r A1, Lemma A2 and Lemma A3 imply that F is strictly contracting. Moreover, z0 = , 0 is a fixed point of F . Then, from the fixed point theorem, δ System (A.1) has unique equilibrium z0 . Consequently, the following results hold. Lemma A5. If L0 > 1, then System (A.1) does not admit a coexistence equilibria. Lemma A6. If L0 < 1, then System (A.1) has one or three coexistence equilibria.
Appendix B. Proofs Appendix B.1. Proof of Lemma 1 (1) The resource’s equation is given by: dR(t) αb Hb + αg Hg αb Hb + αg Hg aR2 c1 Hg =r 1+ 1− R−δR2 − Hb − R, 2 dt β γ b+R 1 + c2 Hg + c3 Hb 31
Since αb Hb + αg Hg < γ, then dR γ ≤r 1+ R. dt β Therefore, r 1+
R(t) ≤ R(0)e
γ! β .
Since r ≤ 0, then R(t) → 0. If R(t) → 0, then Hb → 0 and Hg → 0. (2) (i) Let X(t) = (R, Hb , Hg ) (t) be a solution of System (1). Let suppose that the solution X starts in the positive orthant R3+ . To prove that system (1) is positively invariant, we only need to show that at any border of R3+ , either the resulting vector field stays on the border or is pointing inside R3+ . Case 1: {R = 0}. At this border and following the first equation of System (1), dR one has = 0. Consequently, the solution X cannot crosses R3+ throught the dt border {R = 0}. dHb Case 2: {Hb = 0}. Following the second equation of system (1), one has = 0. dt 3 Therefore, the solution X cannot crosses R+ throught the border {Hb = 0}. dHg Case 3: {Hg = 0}. Following the third equation of system (1), one has = 0. dt 3 Therefore, the solution X cannot crosses R+ through the border {Hg = 0}. Consequently the part 1 of Lemma 1 holds. (ii) Let us now prove that the set Ωε , is a positively invariant region and absorbing set for the system (1). The resource’s equation is given by: dR(t) αb Hb + αg Hg αb Hb + αg Hg aR2 = r 1+ 1− R − δR2 − Hb dt β γ b + R2 c1 Hg − R, 1 + c2 Hg + c3 Hb which implies dR(t) αb Hb + αg Hg αb Hb + αg Hg ≤r 1+ 1− R − δR2 . dt β γ Since αb Hb + αg Hg < γ, then dR(t) αb Hb + αg Hg ≤r 1+ R − δR2 . dt β which implies
dR(t) γ ≤r 1+ R − δR2 dt β
Then, lim sup R(t) ≤ t−→+∞
32
r γ 1+ , δ β
and hence for ε > 0,
r γ R(t) ≤ 1+ + ε. δ β
Since 0 < αb Hb + αg Hg < γ, consequently the part 2 of Lemma 1 holds. This completes the proof of Lemma 1. Appendix B.2. Proof of Proposition 1 Equilibria in System (1) are obtained by solving the following system: αb Hb + αg Hg αb Hb + αg Hg c1 Hg R aR2 2 r 1 + 1 − R − δR − − H = 0, 2 b β γ 1 + c H + c H b + R 2 g 3 b aR2 Hb − µb Hb = 0, e b b + R2 c1 RHg eg − µg Hg = 0. 1 + c2 Hg + c3 Hb (B.1) r 1. The trivial equilibria of system (1) are given by E0 = (0, 0, 0) and E1 = , 0, 0 . δ 2. When R 6= 0, Hb 6= 0, and Hg = 0 αb Hb αb Hb aR Hb = 0, (i) r 1+ 1− − δR − β γ b + R2 (B.2) 2 aR eb − µb = 0. (ii) b + R2 r bµb eb a ˜= From (ii ) one has R with S01 > 1 where S01 = . eb a − µb µb ˜ in (i ) gives: Substituting R ! ˜ 1 βγ βγ δ R ˜ b2 − ˜b + H (γ − β) − σ˜ H − 1 = 0. (B.3) αb rαb R αb2 r ˜ β δR 1− and ηa = , then equation (B.3) leads to: γ r ˜ 2 − 1 (γ − β)(1 − 1 )H ˜ b + βγ (ηa − 1) = 0. H b αb ηs αb2
rαb Let us set ηs = βσR˜
The discriminant of equation (B.4) is given by: 2 1 1 4βγ 2 ∆ = 2 (γ − β) 1 − − 2 (ηa − 1) . αb ηs αb Let us set N0 =
4βγηa . 1 2 2 (γ − β) (1 − ) + 4βγ ηs 33
(B.4)
• If N0 > 1 i.e. ∆ < 0, then (B.4) does not admit a real solution. ˜ 0 defined by: • If N0 = 1 i.e. ∆ = 0, then (B.4) admits a unique solution H b 1 (γ − β) 1 − ηs ˜0 = H . b 2αb ˜ H ˜ 0 , 0) is ecologically meaningful whenever N0 = 1, S01 > 1 Equilibrium E2 = (R, b and ηs > 1. ˜ 1 and H ˜ 2 given by: • If N0 < 1 i.e. ∆ > 0, then (B.4) admits two solutions H b b ˜1 = H b
˜2 = H b
r
1 (γ − β)(1 − ) − ηs 1 (γ − β)(1 − ) + ηs
r
(γ − β)2 (1 −
1 2 ) − 4βγ(ηa − 1) ηs
2αb
(γ − β)2 (1 −
1 2 ) − 4βγ(ηa − 1) ηs
2αb
,
.
˜ H ˜ 1 , 0) is ecologically meaningful whenever N0 < 1, ηs > 1 Equilibrium E3 = (R, b ˜ H ˜ 2 , 0) is ecologically meaningful and ηa > 1. Similarly, equilibrium E4 = (R, b 1 whenever (N0 < 1, S0 > 1 and ηa < 1) or (N0 < 1, S01 > 1, ηs > 1 and ηa > 1). Table 1: Existence of solutions of equation (B.4) N0 >1 =1 <1 <1 <1 <1
ηs >1 >1 >1 <1 <1
ηa
E2 √
>1 <1 <1 >1
E3
E4
√
√ √ √ √
Number of equilibrium 0 1 2 1 1 1
3. When R 6= 0, Hb = 0 and Hg 6= 0 αg Hg αg Hg c1 eg Hg r 1+ 1− − δR − = 0; (a) β γ 1 + c2 Hg
(B.5)
c1 eg R − µg = 0. (b) 1 + c2 Hg ¯ = µ g 1 + c2 H ¯g . Equation (b) implies R eg c1 ¯ Substituting R in (a) gives:
¯ 2 + βγc2 + αg (γ − β) (1 − L1 ) H ¯ g + βγ (L0 − 1) = 0, (B.6) ¯ 3 + γ − β (1 − L2 ) H H g g αg L2 αg c2 L1 αg2 c2 where L0 =
δµg βγc2 + αg (γ − β) ; L1 = 2βγδµg c2 βγc1 and L2 = reg c1 + reg c1
r
αg c2
γ−β
+
βγδµg c2 rαg eg c1
.
We summarize the problem of existence of solutions of (B.6) in the following table. 34
Table 2: Existence of solutions of equation (B.6) L2 <1 <1 >1 <1 >1 >1 >1 <1
L1 <1 >1 <1 >1 <1 >1 >1 <1
L0 <1 >1 >1 <1 <1 >1 <1 >1
Number of sign changes Number of positive solutions 1 1 2 2 or 0 2 2 or 0 1 1 3 3 or 1 2 2 or 0 1 1 0 0
Using Lemmas A5 and A6, Table 2 becomes Table 3: Existence of solutions of equation (B.6) L2 <1 <1 >1 >1
L1 <1 >1 <1 >1
L0 <1 <1 <1 <1
Number of sign changes Number of positive solutions 1 1 1 1 3 3 or 1 1 1
4. When R 6= 0, Hb 6= 0 and Hg 6= 0 α aR α c1 Hg b Hb + αg Hg b Hb + αg Hg r 1+ − 1− − δR − Hb = 0; (i) β γ 1 + c2 Hg + c3 Hb b + R2 aR2 e − µb = 0; (ii) b b + R2 c1 eg R − µg = 0. (iii) 1 + c2 Hg + c3 Hb (B.7) r ∗ e c R − µ (1 + c2 Hg∗ ) bµ g 1 g b Equations (ii) and (iii) give respectively R∗ = and Hb∗ = . eb a − µb c3 µ g e g c1 R ∗ . Hb∗ > 0 ⇔ eg c1 R∗ − µg (1 + c2 Hg∗ > 0 i.e. N ∗ > 1 with N ∗ = µg (1 + c2 Hg∗ ) Substituting R∗ and Hb∗ in equation (i) gives: D2 Hg∗2 + D1 Hg∗ + D0 = 0, rαg2 (1 − m1 )2 ; βγ rαg (1 − m1 )(γ − β) 2rαb αg (1 − m1 )(1 − m0 ) σR∗ c2 µg D1 = − + − ; βγ βγc3 c3 eg R ∗ rαb (m0 − 1)(γ − β) rαb2 (m0 − 1)2 σR∗ eg c1 R∗ σR∗ D0 = r + − − δR∗ − + , βγc3 βγc3 µ g c3 c3 α b c2 eg c1 R∗ m0 = and m1 = . µg αg c3 where D2 = −
35
(B.8)
We summarize the problem of existence of solutions of equation (B.8) in the following Table. Table 3: Existence of solutions of equation (B.8) D2 <0 <0 <0 <0
D1 <0 <0 >0 >0
D0 <0 >0 >0 <0
Number of sign changes Number of positive solutions 0 0 1 1 1 1 2 2 or 0
Appendix B.3. Proof of Theorem 1 1. It is easy to verify that the eigenvalues of the Jacobian matrix at E0 = (0, 0, 0) are r, −µb and −µg . Therefore, it is an unstable hyperbolic saddle critical point. • If r < 0, then E0 is locally asymptotically stable.
• If r > 0, then E0 is a saddle-point and therefore unstable. r eb ar2 2. The eigenvalues of the Jacobian matrix at E1 = ( , 0, 0) are −r, 2 − µb = δ bδ + r2 e g c1 r 1 µb S0 − 1 and − µg = µg − 1 . Hence, equilibrium E1 is locally asymptotδ L0 ically stable if S0 < 1 and L0 > 1. The instability of E1 follows if one of the previous conditions is not satisfied. 3. The jacobian matrix at the point E˜ is given by: m ˜ 11 m ˜ 12 m ˜ 13 ˜ = m ˜ 21 0 0 , J(E) 0 0 m ˜ 33 where,
˜H ˜ b (R ˜ 2 − b) ˜2 aR rαb ˜ m ˜b R ˜ − aR ; − δ R; ˜ 12 = γ − β − 2αb H ˜ 2 )2 ˜2 βγ (b + R b+R ˜˜ ˜H ˜b ˜ 2abR eg c1 R m ˜ 21 = ;m ˜ 33 = − µg = µg S0RHb − 1 . ˜ 2 )2 ˜b (b + R 1 + c3 H The characteristic polynomial is given by: P (λ) = (m ˜ 33 − λ) λ2 − m ˜ 11 λ − m ˜ 21 m ˜ 12 . m ˜ 11 =
The characteristic polynomial admits negative roots if
m ˜ 33 < 0, m ˜ 11 < 0 and m ˜ 21 m ˜ 12 < 0. ˜ ˜
m ˜ 33 < 0 ⇔ S0RHb < 1, ˜H ˜ b (R ˜ 2 − b) ˜2H ˜b aR aR ˜ < 0 ⇔ S1R˜ H˜b < 1 with S1R˜ H˜b = m ˜ 11 = − δR . ˜ 2 )2 ˜ b + δ(b + R ˜ 2 )2 (b + R abH m ˜ 21 m ˜ 12 < 0 ⇔ m ˜ 12 < 0 because m ˜ 21 > 0, ˙ ∂R m ˜ 12 = = ψ 0 (Hb ) − G0 (Hb ) with ∂Hb αb Hb αb Hb ˜ + σ ˜ Hb . 1− and G(Hb ) = δ R ψ(Hb ) = r 1 + R β γ 36
Moreover,
dR = 0 is equivalent to dt ψ(Hb ) − G(Hb ) = 0.
therefore
(B.9)
dR = 0 is equivalent to dt ψ(Hb ) − G(Hb ) = 0.
(B.10)
Furthermore, by using relation (B.10) we deduce graphically the stability of equilib˜ see Figure B.10. rium E,
Figure B.10: We have three different configurations. Case (i): There is no equilibria. Case (ii): The unique equilibrium is locally asymptotically stable. Case (iii): There exist two equilibria and one is locally asymptotically stable.
4. The jacobian matrix at the point E¯ is given by: m ¯ 11 m ¯ 12 m ¯ 13 ¯ = 0 m ¯ 22 0 , J(E) m ¯ 31 0 m ¯ 33 where,
¯ ¯ rαb ¯g R ¯ − aR + c1 c3 Hg R; ¯ ¯ m γ − β − 2αg H m ¯ 11 = −δ R; ¯ 12 = ¯ 2 (1 + c2 H ¯ g )2 βγ b+R m ¯ 13 =
¯ ¯ rαg ¯g R ¯ − c1 R + c1 c2 Hg R; ¯ γ − β − 2αg H ¯ g (1 + c2 H ¯ g )2 βγ 1+H
¯2 ¯g ¯g aeb R eg c1 H µ g c2 H ¯ R − µ = µ S − 1 ; m ¯ = ; m ¯ = − b b 31 33 0 ¯2 ¯g ¯g . b+R 1 + c2 H (1 + c2 H The characteristic polynomial, given by: m ¯ 22 =
P (λ) = (m ¯ 22 − λ) [λ2 − (m ¯ 11 + m ¯ 33 )λ + m ¯ 11 m ¯ 33 − m ¯ 13 m ¯ 31 ], admits negative roots if m ¯ 22 < 0, m ¯ 11 + m ¯ 33 < 0, 37
and m ¯ 11 m ¯ 33 − m ¯ 13 m ¯ 31 > 0.
m ¯ 11 + m ¯ 33 ¯
¯g µ g c2 H ¯ = − δR + ¯ g < 0. (1 + c2 H
Observe that m ˜ 22 < 0 ⇔ S0R < 1, m ¯ 11 m ¯ 33 − m ¯ 13 m ¯ 31 > 0 if m ¯ 13 < 0, m ¯ 13 =
∂ R˙ = ψ 0 (Hg ) − T 0 (Hg ), with ∂Hg αg Hg c1 Hg αg Hg 1− and T (Hg ) = δR + . ψ(Hg ) = r 1 + β γ 1 + c2 Hg
Moreover,
dR = 0 is equivalent to dt ψ 0 (Hg ) − T 0 (Hg ) = 0.
(B.11)
Therefore, m ¯ 13 < 0 is equivalent to ψ 0 (Hg )−T 0 (Hg ) < 0. Furthermore, by using relation ¯ see Figure B.11. (B.11), we deduce graphically the stability of equilibrium E,
¯ Case (a): The unique equilibrium is Figure B.11: We have different configurations for the stability of E. locally asymptotically stable. Case (b): there exist three equilibria and two are locally asymptotically stable.
38
Appendix B.4. Proof of Theorem 2 The characteristic polynomial of the jacobian matrix at the point E˜ is given by: P (λ) = (m ˜ 33 − λ) λ2 − m ˜ 11 λ − m ˜ 21 m ˜ 12 ,
˜ H˜b (R ˜ 2 − b) aR ˜ Observe that where m ˜ 11 = − δ R. ˜ 2 )2 (b + R P (λ) = 0 ⇔ λ = m ˜ 33 or λ2 − m ˜ 11 λ − m ˜ 21 m ˜ 12 = 0. (B.12) From expression of m ˜ 11 , E˜ is the only candidate for a Hopf bifurcation. It follows from expression of m ˜ 11 that if a Hopf bifurcation occurs, there exists δh that m ˜ 11 = 0 or equivalently ˜ H˜b (R ˜ 2 − b) aR ˜ = 0 i.e. − δR 2 2 ˜ (b + R ) ˜ 2 − b) aH˜b (R . δh = ˜ 2 )2 (b + R The discriminant of equation (B.12) is given by ∆ = m ˜ 211 + 4m ˜ 21 m ˜ 12 . 1. If δ 6= δh and ∆ < 0, then the eigenvalues of P (λ) are given by: λ1,2 (δ) = α1,2 (δ) + iβ1,2 (δ). 2. If δ = δh and ∆ < 0, then λ1,2 (δh ) = iβ1,2 (δh ). 3. Moreover,
∂α1,2 ˜ 6= 0. (δh ) = −R ∂δ
∂α1,2 (δh ) < 0, and the conclusion follows : ∂δ • δ = δh is a bifurcation value.
Thus,
• there exists δh1 > δh such that for all δ ∈ [δh , δh1 ], the equilibrium point E˜ is a stable focus. ˜ there exists δh2 < δh such that for every δ ∈ • for every neighborhood U of E, [δh2 , δh [, the equilibrium E˜ is an unstable focus surrounded by a√stable limit cycle contained in U , whose amplitude increases and is of the order δh − δ.
Therefore, in Theorem 2, we have a supercritical Hopf bifurcation with the appearance of a stable limit cycle.
Appendix B.5. Proof of Theorem 3 The jacobian matrix at the point E ∗ is given by: ∗ m11 m∗12 m∗13 0 , J(E ∗ ) = m∗21 0 ∗ ∗ m31 m32 m∗33 39
where,
αb Hb∗ + αg Hg∗ c1 Hg∗ 2abR∗ ∗ 1− −2δR − − Hb∗ ; ∗ ∗ ∗2 2 γ 1 + c2 Hg + c3 Hb (b + R )
m∗11
αb Hb∗ + αg Hg∗ =r 1+ β
m∗12
∗ c1 c3 Hg∗ rαb aR∗2 ∗ ∗ R∗ ; = γ − β − 2(αb Hb + αg Hg ) R − + ∗ 2 ∗2 ∗ βγ b+R (1 + c2 Hg + c3 Hb )
m∗13 =
m∗21 =
rαb 2eb abR∗ ∗ c1 R∗ (1 + c3 Hb∗ ) ∗ = γ − β − 2(αb Hb∗ + αg Hg∗ ) R∗ − ; m H 31 βγ (1 + c2 Hg∗ + c3 Hb∗ )2 (b + R∗2 )2 b eg c1 c3 Hg∗ R∗ µg c2 Hg∗ eg c1 Hg∗ ∗ ∗ = − = − ; m ; m . 32 33 1 + c2 Hg∗ + c3 Hb∗ (1 + c2 Hg∗ + c3 Hb∗ )2 1 + c2 Hg∗ + c3 Hb∗
The characteristic polynomial is given by: P (λ) = −λ3 + (m∗11 + m∗33 )λ2 + (m∗13 m∗31 + m∗12 m∗21 − m∗11 m∗33 )λ + m∗12 m∗21 m∗33 − m∗13 m∗21 m∗32 . The Routh-Hurwitz conditions for stability of this equilibrium are C1 = m∗11 + m∗33 < 0, C2 = (m∗11 + m∗33 )(m∗11 m∗33 − m∗13 m∗31 − m∗12 m∗21 ) + (m∗12 m∗21 m∗33 − m∗13 m∗21 m∗32 ) < 0, C3 = m∗12 m∗21 m∗33 − m∗13 m∗21 m∗32 > 0. When these conditions are satisfied, a coexistence equilibrium is locally asymptotically stable. Appendix B.6. Proof of Theorem 4 Since ω0 is a solution of equation (11), the characteristic equation (8) has the pair of purely imaginary roots ±iω0 . From equations (9) and (10), we have cos ω0 τ1 =
B3 ω0 4 + (B1 B4 − B2 B3 )ω0 2 . B32 ω0 2 + B42
0 Hence, τ1,n for n = 0, 1, ... as a function of ω0 is given by 1 B3 ω0 4 + (B1 B4 − B2 B3 )ω0 2 2πn 0 τ1,n = arccos + . 2 2 ω0 B3 ω0 2 + B4 ω0
(B.13)
For τ1 = 0, Theorem 1 ensures that E˜ is locally asymptotically stable. Hence, by Butler’s 0 lemma Freedman et al. 1983 [12] E˜ remains stable up to the minimum value of τ1,n , obtained 0 here thus for n = 0 i.e. for τ1 < τ1,0 , so that 0 0 = τ1,0 . τ10 = min τ1,n n≥0
The theorem will be completely proven if we can show that d(Reξ(τ1 )) > 0. sign dτ1 ξ=iω0 40
Differentiating equation (18) with respect to τ1 yields dξ 2 3ξ + 2B1 ξ + B2 + B3 e−ξτ1 − τ1 (B3 ξ + B4 )e−ξτ1 = ξ(B3 ξ + B4 )e−ξτ1 , dτ1
which implies that −1 dξ 3ξ 2 + 2B1 ξ + B2 + B3 e−ξτ1 − τ1 (B3 ξ + B4 )e−ξτ1 . = dτ1 ξ(B3 ξ + B4 )e−ξτ1 Since (B3 ξ + B4 )e−ξτ1 = −(ξ 3 + B1 ξ 2 + B2 ξ), we obtain −1 dξ 3ξ 2 + 2B1 ξ + B2 B3 τ1 = + − , 3 2 dτ1 −ξ(ξ + B1 ξ + B2 ξ) ξ(B3 ξ + B4 ξ =
ξ 2 + B1 ξ + B2 + 2ξ 2 + B1 ξ B3 τ1 + − , 2 2 −ξ (ξ + B1 ξ + B2 ) ξ(B3 ξ + B4 ξ
=
2ξ 2 + B1 ξ −B4 τ1 + − . −ξ 2 (ξ 2 + B1 ξ + B2 ) ξ 2 (B3 ξ + B4 ξ
Thus, ( −1 ) d(Reξ(τ1 )) dξ sign = sign Re , dτ1 dτ1 ξ=iω0 ξ=iω0 2ξ 2 + B1 ξ −B4 = sign Re + Re 2 , −ξ 2 (ξ 2 + B1 ξ + B2 ) ξ=iω0 ξ (B3 ξ + B4 ) ξ=iω0 1 2B32 ω0 8 + (3B42 + B32 θ20 )ω0 6 + 2B42 θ20 ω0 4 + B42 B22 ω0 2 . = sign ω0 2 [B12 ω0 4 + (ω0 3 − B2 ω0 )2 ] [B42 + B32 ω0 2 ] Since θ20 > 0, we have
d (Reξ(τ1 )) > 0. dτ1 ξ=iω0 Hence, the transversality condition is satisfied and a Hopf bifurcation occurs at τ1 = τ10 . This achieves the proof. Appendix B.7. Proof of Theorem 5 ˜ x2 (t) = Hb (t) − H ˜ b and x3 (t) = Hg (t); then system (1) is equivalent Let x1 (t) = R(t) − R, to the following three dimensional system: dx1 (t) ˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b x1 (t) + ψH (H ˜ b , 0)R ˜ − P (R) ˜ x2 (t) = ψ( H b dt ˜ ˜ ˜ ˜ + ψ ( H , 0) R − g ( H , 0) R x (t) + f x (t), x (t), x (t) ; Hg b Hg b 3 1 1 2 3 dx2 (t) ˜ ˜ = e P ( R) H x (t − τ ) + f x (t), x (t), x (t), x (t − τ ) ; b b 1 1 2 1 2 3 1 1 dt dx3 (t) = e g (H ˜ b , 0) − µg x3 (t) + f3 x1 (t), x2 (t), x3 (t) , g Hg dt 41
(B.14)
where 2 ˜ ˜ ˜ f1 x1 (t), x2 (t), x3 (t) = ψ x2 (t) + Hb , x3 (t) x1 (t) + R − δ x1 (t) + R ˜ − g x2 (t) + H ˜ b , x3 (t) x1 (t) + R ˜ − P x1 (t) + R h ˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b x1 (t) + ψH (H ˜ b , 0)R ˜ − P (R) ˜ x2 (t) − ψ(H b +
i ˜ b , 0)R ˜ − gHg (H ˜ b , 0)R ˜ x3 (t) ; ψHg (H
f2 x1 (t), x2 (t), x3 (t), x1 (t − τ2 )
˜ ˜ b ) − µb (x2 (t) + H ˜ b) = eb P x1 (t) + R (x2 (t) + H
˜ H ˜ b x1 (t − τ1 ); − eb P (R)
and ˜ b , x3 (t) x1 (t) + R ˜ − µg x3 (t) f3 x1 (t), x2 (t), x3 (t) = eg g x2 (t) + H ˜ − eg gHg (Hb , 0) − µg x3 (t).
Let τ1 = τ10 + µ, then µ = 0 is the Hopf bifurcation value of system (1) at the positive ˜ H ˜ b , 0). Since system (1) is equivalent to system (B.14), in the following equilibrium E˜ = (R, discussion we shall consider mainly the system (B.14). In system (B.14), let ui (t) = xi (τ1 t). Then, system (B.14) can be rewritten as follows: h du1 (t) 0 ˜ b , 0)R ˜ − P (R) ˜ u2 (t) ˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b u1 (t) + ψH (H + µ) ψ( H = (τ 1 b dt i ˜ ˜ ˜ ˜ ( H , 0) R − g ( H , 0) R u (t) + f u (t), u (t), u (t) ; + ψ 1 2 3 b 3 1 Hg b Hg h i du2 (t) 0 ˜ ˜ + µ) e P ( R) H u (t − τ ) + f u (t), u (t), u (t), u (t − τ ) ; = (τ 1 2 3 1 1 b b 1 1 2 1 dt h i du3 (t) = (τ 0 + µ) e g (H ˜ , 0) − µ u (t) + f u (t), u (t), u (t) . g H b g 3 3 1 2 3 g 1 dt
(B.15) System (B.15) can be written as functional differential equation (FDE) in C([−1, 0], R3 ) of the form: u(t) ˙ = Lµ (ut ) + f (µ, ut ), (B.16)
where, u(t) = (u1 (t), u2 (t), u3 (t))T , Lµ : C −→ R is the linear operator and f : (R×C) −→ R is the nonlinear operator, given respectively by
s11 s12 s13
Lµ (φ) = (τ10 + µ) 0 0
0 0
0
0 0
φ1 (0) φ2 (0) + t21 0 0 0 φ3 (0) s33 0 0 0 42
φ1 (−1) φ2 (−1) (B.17) φ3 (−1)
˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b , s12 = ψH (H ˜ b , 0)R ˜ − P (R), ˜ where s11 = ψ(H b ˜ ˜ ˜ ˜ ˜ ˜ H ˜b s13 = ψHg (Hb , 0)R − gHg (Hb , 0)R, s33 = eg gHg (Hb , 0) − µg and t21 = eb P 0 (R) and f1 (φ1 (0), φ2 (0), φ3 (0)) , f (φ (0), φ (0), φ (0), φ (−1)) f (φ, ν) = (τ10 + µ) 2 1 2 3 1 f3 (φ1 (0), φ2 (0), φ3 (0))
(B.18)
where φ = (φ1 , φ2 , φ3 )T ∈ C. By the Riesz representation theorem, there exists a (3 × 3) matrix, η(θ, µ), (−1 ≤ θ ≤ 0) whose elements are bounded variation functions such that Z 0 dη(θ, µ)φ(θ) f or φ ∈ C([−1, 0], R3 ). (B.19) Lµ (φ) = −1
In fact, we can choose
η(θ, µ) =
(τ10
s11 s12 s13
+ µ) 0 0
0 0
0
0 0
t21 0 0 0 δ(θ) + s33 0 0 0
where δ is the Dirac function defined by 0, δ(θ) = 1,
θ 6= 0;
δ(θ + 1) ,
(B.20)
(B.21)
θ = 0.
For φ ∈ C([−1, 0], R3 ), we define
dφ(θ) , θ ∈ [−1, 0); dθ A(µ)φ(θ) = Z 0 dη(µ, s)φ(s), θ = 0.
(B.22)
−1
and
R(µ)φ(θ) = Then, the system (B.15) is equivalent to
0,
θ ∈ [−1, 0);
f (µ, θ),
θ = 0.
u(t) ˙ = A(µ)ut + R(µ)ut ,
43
(B.23)
(B.24)
where ut (θ) = u(t + θ) for θ ∈ [−1, 0]. For ψ ∈ C([−1, 0], R3 ), we define dψ(s) − , s ∈ (0, 1]; ds A∗ ψ(s) = Z 0 dη T (t, 0)ψ(−t), s = 0.
(B.25)
−1
and a bilinear inner product
¯ hψ(s), φ(θ)i = ψ(0)φ(0) −
Z
0
−1
Z
θ
ξ=0
¯ − θ)dη(θ)φ(ξ)dξ, ψ(ξ
(B.26)
where η(θ) = η(θ, 0). Then, A(0) and A∗ are adjoint operators. By Theorem 4, we know that ±iω0 τ10 are eigenvalues of A(0). Thus, they are also eigenvalues of A∗ . Let us consider q(θ) as the eigenvector of A(0) corresponding to +iω0 τ10 and q ∗ (θ) the eigenvector of A∗ corresponding to −iω0 τ10 . T T 0 0 Then, q(θ) = 1, β3 , γ3 eiθω0 τ1 and q ∗ (s) = D 1, β3∗ , γ3∗ eiω0 τ1 s . From the above discussion, it is easy to know that A(0)q(0) = iω0 τ10 q(0) and A∗ (0)q ∗ (0) = −iω0 τ10 q ∗ (0). That is
τ10
and
τ10
Then,
s11 s12 s13
0 0
s11 0
and
0
0 0
0 0 q(0) + τ1 t21 0 0 s33 0 0 0
0 0
0
s12 0 0 s13 0 s33
0 t21 0
∗ 0 q (0) + τ1 0 0
iω0 − s11
−t21 e−iω0 τ10 0
0 0
−s12
−s13
iω0
0
0
iω0 − s33
44
q(−1) = iω0 τ10 q(0),
∗ 0 ∗ 0 q (−1) = −iω0 τ1 q (0). 0
1
0
β3 = 0 γ3 0
(B.27)
(B.28)
(B.29)
0
−iω0 − s11 −t21 e−iω0 τ1
0
−s12
−iω0
0
−s13
0
−iω0 − s33
Thus, we can easily obtain
1
∗
D1 1, β3∗ , γ3∗
q (s) = where
T
0
∗ β3 = 0 . ∗ γ3 0
T 0 q(θ) = 1, β3 , 0 eiθω0 τ1 ,
and
(B.30)
(B.31)
0
eiω0 τ1 s ,
(B.32)
˜ b , 0)R ˜ ˜ b , 0)R ˜ + gHg (H ˜ b , 0)R ˜ + P (R) ˜ ˜ H ˜ b e−iω0 τ10 s −ψHg (H −ψHb (H eb P 0 (R) , β3∗ = , γ3∗ = β3 = − ¯ g )R ¯ − µg . iω0 iω0 iω0 + eg gHg (0, H In order to assure hq ∗ (s), q(θ)i = 1, we need to determine the value of D1 . From (B.26), we have Z 0Z θ ∗ ∗ q¯∗ (ξ − θ)dη(θ)φ(ξ)dξ, hq (s), q(θ)i = q¯ (0)q(0) − −1 0
Z
ξ=0 θ
Z
T ¯ 1, β¯∗ , γ¯ ∗ e−iω0 τ10 (ξ−θ) dη(θ) 1, β3 , 0 eiω0 τ10 ξ dξ, D 3 3 −1 ξ=0 0 0 0 ˜ H ˜ b 0 0 (−e−iω0 τ10 )q(0), = q¯∗ (0)q(0) − q¯∗ τ10 eb P 0 (R) 0 0 0 = q¯∗ (0)q(0) −
¯ 1 [1 + β3 β¯3∗ + τ10 β¯3∗ eb P 0 (R) ˜ H ˜ b e−iω0 τ10 ]. = D Therefore,
1 , ∗ 0 ¯∗ ¯ ˜ H ˜ b e−iω0 τ10 1 + β3 β3 + τ1 β3 eb P 0 (R) 1 . = ∗ 0 ∗ ¯ ˜ H ˜ b eiω0 τ10 1 + β3 β3 + τ1 β3 eb P 0 (R)
D¯1 =
(B.33)
D1
(B.34)
Using the same notations as in Hassard et al. 1981 [16], we first compute the coordinates to describe the center manifold C2 (0) at µ = 0. Let ut be the solution of Equation (B.24) when µ = 0. Define z0 (t) = hq ∗ , ut i, W (t, θ) = ut (θ) − 2Re(z0 (t)q(θ)) = ut (θ) − z0 (t)q(θ) + z¯0 (t)¯ q (θ) . (B.35) On the center manifold C2 (0), we have
W (t, θ) = W (z0 , z¯0 , θ), 45
(B.36)
where, W (z0 , z¯0 , θ) = W20 (θ)
z02 z¯2 z3 + W11 (θ)z0 z¯0 + W02 (θ) 0 + W30 (θ) 0 + ..., 2 2 6
(B.37)
where z0 and z¯0 are local coordinates for center manifold C2 (0) in the direction of q ∗ and q¯∗ . Note that W is real if ut is real. We only consider real solutions. For solution ut ∈ C2 (0) of (B.15), since µ = 0, we have def z˙0 (t) = iω0 τ10 z0 + q¯∗ (0)f 0, W (z0 , z¯0 , 0) + 2Re(z0 q(θ)) = iω0 τ10 z0 + q¯∗ (0)f0 (z0 , z¯0 ). (B.38)
We rewrite this equation as
z˙0 (t) = iω0 τ10 z0 + p(z0 , z¯0 ),
(B.39)
where, p(z0 , z¯0 ) = q¯∗ (0)f0 (z0 , z¯0 ) = p20 (θ)
z02 z¯2 z 2 z¯0 + p11 (θ)z0 z¯0 + p02 (θ) 0 + p21 (θ) 0 + ..., 2 2 2
(B.40)
Remark 7. In what follows, our objective is to compute coefficients p20 , p11 , p02 and p21 of p(z0 , z¯0 ). These coefficients will be use to find the direction of Hopf bifurcation. T 0 Let us recall, ut (θ) = u1t (θ), u2t (θ), u3t (θ) and q(θ) = 1, β3 , 0 eiω0 τ1 θ . So from (B.35) and (B.37), it follows that ut (θ) = W (t, θ) + 2Re(z0 (t)q(θ)), z¯2 z02 + W11 (θ)z0 z¯0 + W02 (θ) 0 2 2 T T 0 0 + 1, β3 , 0 eiω0 τ1 θ z0 + 1, β¯3 , 0 e−iω0 τ1 θ z¯0 + ... = W20 (θ)
Then,
(1)
z02 z¯2 (1) (1) + W11 (0)z0 z¯0 + W02 (0) 0 + ..., 2 2
u1t (0)
= z0 + z¯0 + W20 (0)
u2t (0)
z z¯ (2) (2) (2) = β3 z0 + β¯3 z¯0 + W20 (0) 0 + W11 (0)z0 z¯0 + W02 (0) 0 + ..., 2 2
u3t (0)
= W20 (0)
2
(3)
(B.41)
2
z02 z¯2 (3) (3) + W11 (0)z0 z¯0 + W02 (0) 0 + ..., 2 2
0
0
(1)
u1t (−1) = e−iω0 τ1 z0 + eiω0 τ1 z¯0 + W20 (−1)
z02 z¯2 (1) (1) + W11 (−1)z0 z¯0 + W02 (−1) 0 + ..., 2 2 2
2
z z¯ 0 0 (2) (2) (2) u2t (−1) = β3 e−iω0 τ1 z0 + β¯3 eiω0 τ1 z¯0 + W20 (−1) 0 + W11 (−1)z0 z¯0 + W02 (−1) 0 + ..., 2 2 (3)
u3t (−1) = W20 (−1)
z02 z¯2 (3) (3) + W11 (−1)z0 z¯0 + W02 (−1) 0 + .... 2 2 46
(B.42)
Also taking into account Eq.(B.18), it follows that
where,
p(z0 , z¯0 ) = q¯∗ f0 (z0 , z¯0 ), = q¯∗ f (0,ut ), ¯ 1 1, β¯3 , 0 f (0, ut ), = τ10 D m1 ¯ 1 1, β¯3 , 0 × m2 , = τ10 D m3
(B.43)
2 r h ˜ b αb u2t (0)u1t (0) + αg u3t (0)u1t (0) − αb u2t (0) + αg u3t (0) γ − β − 2αb H βγ i ˜ − δu2 (0) − P 0 (R)u ˜ 2t (0)u1t (0) − 1 P 00 (R)u ˜ 2t (0)u2 (0) − 1 P 00 (R) ˜ H ˜ b u2 (0) u1t (0) + R 1t 1t 1t 2 2 2 ˜ b , 0)u3t (0)u1t (0) − 1 ∂ g (H ˜ b , 0) u2 (0)u1t (0) + Ru ˜ 2 (0) − gHg (H 3t 3t 2 ∂Hg2 ∂ 2g ˜ ˜ + (Hb , 0) u3t (0)u2t (0)u1t (0) + Ru3t (0)u2t (0) ∂Hb Hg
m1 =
˜ 2t (0)u2 (−1) + 1 eb P 00 (R) ˜ H ˜ b u2 (−1) + ... ˜ 2t (0)u1t (−1) + 1 eb P 00 (R)u m2 = eb P 0 (R)u 1t 1t 2 2 2 ˜ b , 0)u2t (0)u1t (0) + 1 eg ∂ g (H ˜ b , 0) u23t (0)u1t (0) + Ru ˜ 23t (0) m3 = eg gHb (H 2 ∂Hg2 ∂ 2g ˜ ˜ (Hb , 0) u3t (0)u2t (0)u1t (0) + Ru3t (0)u2t (0) + .... + eg ∂Hb Hg
p(z0 , z¯0 ) is given by n h ¯1 r ˜ b αb u2t (0)u1t (0) + αg u3t (0)u1t (0) − αb u2t (0) p(z0 , z¯0 ) = τ10 D γ − β − 2αb H βγ 2 i ˜ − δu21t (0) − P 0 (R)u ˜ 2t (0)u1t (0) − 1 P 00 (R)u ˜ 2t (0)u21t (0) + αg u3t (0) u1t (0) + R 2 1 00 ˜ ˜ 2 1 ∂ 2g ˜ 2 2 ˜ ˜ − P (R)Hb u1t (0) − gHg (Hb , 0)u3t (0)u1t (0) − (Hb , 0) u3t (0)u1t (0) + Ru3t (0) 2 2 ∂Hg2 ∂ 2g ˜ b , 0) u3t (0)u2t (0)u1t (0) + Ru ˜ 3t (0)u2t (0) + (H ∂Hb Hg h io 1 1 0 ˜ 00 ˜ 2 00 ˜ ˜ 2 ¯ + β3 eb P (R)u2t (0)u1t (−1) + eb P (R)u2t (0)u1t (−1) + eb P (R)Hb u1t (−1) + ... 2 2
This expression can be extended to
P 00 (R) ˜ H ˜b z02 n 0 ¯ h rαb β3 2˜ 2 0 ˜ ˜ p(z0 , z¯0 ) = 2τ1 D1 γ − β − 2αb Hb − αb Rβ3 − δ + β3 P (R) − 2 βγ 2 47
io ˜ b P 00 (R)e ˜ −2iω0 τ10 ˜ 3 e−iω0 τ10 + 1 eb H + β¯3 eb P 0 (R)β 2
P 00 (R) ˜ H ˜b z¯02 n 0 ¯ h rαb β¯3 2 ˜ ¯2 0 ˜ ¯ ˜ + 2τ1 D1 γ − β − 2αb Hb − αb Rβ3 − δ + β3 P (R) − 2 βγ 2 io ˜ β¯3 eiω0 τ10 + 1 eb H ˜ b P 00 (R)e ˜ 2iω0 τ10 + β¯3 eb P 0 (R) 2
n h rα (β + β¯ ) b 3 3 0 ¯ 2˜ 0 ˜ ¯ ¯ ˜ + z0 z¯0 τ1 D1 γ − β − 2αb Hb − 2αb Rβ3 β3 − 2δ − β3 + β3 )P (R βγ io 0 0 ˜ b P 00 (R) ˜ + β¯3 eb P 0 (R)Re{( ˜ ˜ b P 00 (R) ˜ − H β¯3 e−iω0 τ1 + β3 eiω0 τ1 )} + eb H z02 z¯0 n 0 ¯ h r (2) (1) (2) ˜ b 2αb W11 τ1 D1 γ − β − 2αb H (0) + 2β3 αb W11 (0) + αb W20 (0) 2 βγ (1) (3) (3) + β3 αb W20 (0) + 2αg W11 (0) + αg W20 (0) − 2 αb2 β32 + 2β3 β¯3 αb2 +
+
(2) 2W11 (0)
+
(2) W20 (0)
(3) (3) ˜ + 2αb αg Rβ3 W11 (0) + αb αg W20 (0)
(1) (1) (2) (1) ˜ 2W11 − 2δ 2W11 (0) + W20 (0) − P 0 (R) (0) + 2β3 W11 (0)
(2) (1) (1) (1) ˜ β¯3 + 2β3 + 2H ˜ b W11 ˜ b W20 + W20 (0) + β¯3 W11 (0) − P 00 (R) (0) + H (0) ∂ 2g (3) (3) (3) (3) ˜ b , 0) 2W11 ˜ b , 0) 2β3 W11 − gHg (H (0) + W20 (0) − (H (0) + W20 (0) ∂Hb Hg 0 0 (2) (1) (2) (1) ˜ W11 + β¯3 eb P 0 (R) (0)e−iω0 τ1 + β3 W11 (−1) + W20 (0)eiω0 τ1 + β¯3 W11 (−1)
io 1 0 0 0 (1) (1) ˜ 2H ˜ b W11 ˜ b W20 eb P 00 (R) (−1)e−iω0 τ1 + H (−1)eiω0 τ1 + 2β¯3 e−2iω0 τ1 + 4β3 2 Comparing the coefficients with (B.40), gives h 00 ˜ ˜ ¯ 1 rαb β3 γ − β − 2αb H ˜ b − α2 Rβ ˜ 2 − δ + β3 P 0 (R) ˜ − P (R)Hb p20 = 2τ10 D 3 b βγ 2 i ˜ b P 00 (R)e ˜ −2iω0 τ10 , ˜ 3 e−iω0 τ10 + 1 eb H + β¯3 eb P 0 (R)β 2 +
p02
h rα β¯ P 00 (R) ˜ H ˜b b 3 2 ˜ ¯2 0 ˜ ¯ ˜ γ − β − 2αb Hb − αb Rβ3 − δ + β3 P (R) − = βγ 2 i ˜ β¯3 eiω0 τ10 + 1 eb H ˜ b P 00 (R)e ˜ 2iω0 τ10 , + β¯3 eb P 0 (R) 2
¯1 p11 = τ10 D
¯1 2τ10 D
h rα (β + β¯ ) b 3 3 ˜ b − 2α2 Rβ ˜ 3 β¯3 − 2δ − β3 + β¯3 )P 0 (R ˜ −H ˜ b P 00 (R) ˜ γ − β − 2αb H b βγ 48
p21
i 0 0 ˜ b P 00 (R) ˜ ˜ , + β¯3 eb P 0 (R)Re{ β¯3 e−iω0 τ1 + β3 eiω0 τ1 } + eb H h r (2) (1) (2) (1) ¯1 ˜ b 2αb W11 = τ10 D γ − β − 2αb H (0) + 2β3 αb W11 (0) + αb W20 (0) + β3 αb W20 (0) βγ (2) (2) (3) (3) + 2αg W11 (0) + αg W20 (0) − 2 αb2 β32 + 2β3 β¯3 αb2 + 2W11 (0) + W20 (0) (3) (3) (1) (1) ˜ 3 W11 + 2αb αg Rβ (0) + αb αg W20 (0) − 2δ 2W11 (0) + W20 (0) (2) (1) (2) (1) ˜ 2W11 − P 0 (R) (0) + 2β3 W11 (0) + W20 (0) + β¯3 W11 (0)
(1) (1) (3) (3) ˜ β¯3 + 2β3 + 2H ˜ b W11 ˜ b W20 ˜ b , 0) 2W11 − P 00 (R) (0) + H (0) − gHg (H (0) + W20 (0) −
∂ 2g (3) (3) ˜ b , 0) 2β3 W11 (H (0) + W20 (0) ∂Hb Hg
0 0 (1) (1) (2) (2) ˜ W11 + β¯3 eb P 0 (R) (0)e−iω0 τ1 + β3 W11 (−1) + W20 (0)eiω0 τ1 + β¯3 W11 (−1) +
1 ˜ e P 00 (R) 2 b
i 0 0 0 (1) (1) ˜ b W20 ˜ b W11 . (−1)eiω0 τ1 + 2β¯3 e−2iω0 τ1 + 4β3 2H (−1)e−iω0 τ1 + H
(B.44)
Since there are W20 and W11 in p21 , we still need to compute them. From (B.24) and (B.35), we have: ˙ W
=
=
def
=
u˙ t − z˙0 q − z¯˙0 q¯,
AW − 2Re(q¯∗ (0)f0 q(θ)),
θ ∈ [−1, 0),
AW − 2Re(q¯∗ (0)f0 q(θ)) + f0 ,
(B.45)
θ = 0,
AW + H(z0 , z¯0 , θ)
where
z02 z¯2 + H11 (θ)z0 z¯0 + H02 (θ) 0 + .... (B.46) 2 2 Substituting the corresponding series into Eq. (B.45), and comparing the coefficients, give H(z0 , z¯0 , θ) = H20 (θ)
(A − 2iω0 τ10 )W20 (θ) = −H20 (θ),
(B.47)
= −H11 (θ).
AW11 (θ) From Eq. (B.45) we know that for θ ∈ [−1, 0),
H(z0 , z¯0 , θ) = −q ∗ f0 q(θ) − q ∗ f¯0 q¯(θ) = −p(z0 , z¯0 )q(θ) − p¯(z0 , z¯0 )¯ q (θ).
(B.48)
Comparing the coefficients with Eq. (B.46), gives H20 (θ) = −p20 q(θ) − p¯02 q¯(θ), 49
(B.49)
H11 (θ) = −p11 q(θ) − p¯11 q¯(θ).
(B.50)
From Eqs. (B.46) and (B.49) and the definition of A, it follows that ˙ 20 (θ) = 2iω0 τ 0 W20 (θ) + p20 q(θ) + p¯02 q¯(θ). W 1
(B.51)
T 0 Notice that q(θ) = 1, β3 , 0 eiω0 τ1 θ . Hence,
i¯ p02 ip20 0 0 iω0 τ10 θ q(0)e + q¯(0)e−iω0 τ1 θ + E3 e2iω0 τ1 θ , (B.52) 0 0 ω 0 τ1 3ω0 τ1 (1) (2) (3) where, E3 = E3 , E3 , E3 ∈ R3 is a constant vector. Similarly, from (B.46) and (B.50)), we obtain W20 (θ) =
W11 (θ) = −
i¯ p11 ip11 0 0 q(0)eiω0 τ1 θ + q¯(0)e−iω0 τ1 θ + E4 , 0 0 ω 0 τ1 ω 0 τ1
(B.53)
(1) (2) (3) where, E4 = E4 , E4 , E4 ∈ R3 is also a constant vector. In what follows, we will seek appropriate E3 and E4 . From the definition of A and (B.46), we obtain Z 0
−1
and
dη(θ)W20 (θ) = 2iω0 τ10 W20 (0) − H20 (0), Z
(B.54)
0
−1
dη(θ)W11 (θ) = −H11 (0),
(B.55)
where, η(θ) = η(0, θ). By (B.45), we have H20 (0) = −p20 q(0) − p¯02 q¯(0) + 2τ10 ×
00 ˜ ˜ rαg β3 ˜ b ) − αb2 Rβ ˜ 32 − (δ + β3 P 0 (R)) ˜ − P (R)Hb (γ − β − 2α H b βγ 2 ˜ 3 e−iω0 τ10 + 1 eb H ˜ b P 00 (R)e ˜ −2iω0 τ10 eb P 0 (R)β 2 0
H11 (0) = −p11 q(0) − p¯11 q¯(0) + 2τ10 ×
,
(B.56)
rαb (β3 + β¯3 ) 2˜ 0 ˜ 00 ˜ ¯ ¯ ˜ ˜ (γ − β − 2αb Hb ) − 2αb Rβ3 β3 − 2δ − (β3 + β3 )P (R) − Hb P (R) βγ . 0 ˜ −iω0 τ10 iω0 τ10 00 ˜ ¯ ˜ eb P (R)Re{(β3 e + β3 e )} + eb Hb P (R) 0 (B.57)
50
Substituting Eqs. (B.52) and (B.56) into Eq. (B.54) with the fact that Z 0 0 0 iω0 τ1 I − eiω0 τ1 dη(θ) q(0) = 0, −1
− iω0 τ10 I −
Z
0
−iω0 τ10
e
−1
(B.58)
dη(θ) q¯(0) = 0.
gives,
h1 0 e2iω0 τ1 dη(θ) E3 = 2τ20 × h2 , 2iω0 τ10 I − −1 0 Z
0
00 ˜ ˜ rαg β3 ˜ b ) − αb2 Rβ ˜ 32 − (δ + β3 P 0 (R)) ˜ − P (R)Hb (γ − β − 2αb H βγ 2 1 00 ˜ −2iω0 τ10 0 ˜ −iω0 τ10 ˜ . + 2 eb Hb P (R)e and h2 = eb P (R)β3 e This leads to 2iω0 − h3 h4 h5 h 1 ˜ H ˜ b e−iω0 τ10 2iω0 eb P 0 (R) E3 = 2 h2 , 0 0 ˜ ˜ 0 0 2iω0 − eg gHg (Hb , 0)R + µg
where h1 =
˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b , h4 = −ψH (H ˜ b , 0)R ˜ + P (R) ˜ where, h3 = ψ(H b ˜ b , 0)R ˜ + gHg (H ˜ b , 0)R. ˜ and h5 = −ψHg (H Solving this system for E3 gives ˜ b , 0)R ˜ + P (R) ˜ −ψHg (H ˜ b , 0)R ˜ + gHg (H ˜ b , 0)R ˜ h1 −ψHb (H 2 (1) h2 2iω0 0 E3 = ∆3 0 ˜ b , 0)R ˜ + µg 0 2iω0 − eg gH (H g
(2)
E3
and
˜ b , 0) + 2δ R ˜ + P 0 (R) ˜ H ˜ b h1 −ψHg (H ˜ b , 0)R ˜ + gHg (H ˜ b , 0)R ˜ 2iω0 − ψ(H 2 ˜ H ˜ b e−iω0 τ10 = eb P 0 (R) h2 0 ∆3 ˜ b , 0)R ˜ + µg 0 0 2iω0 − eg gHg (H
(3)
E3
,
,
˜ b , 0) + 2δ R ˜ + P 0 (R) ˜ H ˜ b −ψH (H ˜ b , 0)R ˜ + P (R) ˜ h1 2iω0 − ψ(H b 2 0 0 −iω τ ˜ H ˜ be 0 1 = eb P (R) 2iω0 h2 = 0, ∆3 0 0 0 51
where, 2iω0 − h3 h4 h5 ˜ H ˜ b e−iω0 τ10 0 2iω0 0 ∆3 = eb P 0 (R) ˜ b , 0)R ˜ + µg 0 0 2iω0 − eg gHg (H
.
Similarly, substituting (B.53) and (B.57) into (B.55) gives h3 −h4 −h5 l 1 ˜ H ˜b 0 −eb P 0 (R) E4 = 2 l2 , 0 0 ˜ ˜ 0 0 eg gHg (Hb , 0)R − µg
rαb (β3 + β¯3 ) ˜ b ) − 2αb2 Rβ ˜ 3 β¯3 − 2δ − (β3 + β¯3 )P 0 (R) ˜ −H ˜ b P 00 (R) ˜ and (γ − β − 2αb H βγ 0 0 ˜ ˜ b P 00 (R). ˜ l2 = eb P 0 (R)Re{( β¯3 e−iω0 τ1 + β3 eiω0 τ1 )} + eb H Hence, ˜ b , 0)R ˜ − P (R) ˜ ψHg (H ˜ b , 0)R ˜ + gHg (H ˜ b , 0)R ˜ l1 ψHb (H 2 (1) , l 0 0 E4 = 2 ∆4 0 ˜ b , 0)R ˜ − µg 0 eg gH (H where l1 =
g
(2)
E4
and
˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b l1 ψHg (0, H ¯ g )R ¯ − gHg (0, H ¯ g )R ¯ ψ(H 2 ˜ H ˜b = −eb P 0 (R) l2 0 ∆4 ˜ b , 0)R ˜ − µg 0 0 eg gHg (H
(3)
E4
where
,
˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b ψH (H ˜ b , 0)R ˜ − P (R) ˜ l1 ψ(H b 2 0 ˜ ˜ = 0, = −e P ( R) H 0 l b b 2 ∆4 0 0 0 h3 −h4 −h5 ˜ H ˜b 0 0 ∆4 = −eb P 0 (R) ˜ b , 0)R ˜ − µg 0 0 eg gHg (H 52
.
Thus, we can determine W20 and W11 from (B.52) and (B.53). Furthermore, g21 in (B.44) can be expressed using the parameters and delay. Thus, we can compute the following values: C2 (0) = µ2
|p02 |2 p21 i 2 p p − 2|p | − + , 20 11 11 2ω0 τ10 3 2
= −
Re{C2 (0)} , Re{ξ 0 (τ10 )}
β2
= 2Re{C2 (0)},
T2
= −
(B.59)
Im(C2 (0)) + µ2 Im(ξ 0 (τ10 )) . ω0 τ10
Appendix B.8. Proof of Lemma 3 It suffices to prove that there exists ω ¯ ∈ R+ such that |P1 |2 − |P0 |2 − |P2 |2 = 0.
(B.60)
Note that equation (B.60) is equivalent to ω ¯6 + ω ¯ 4 (R12 − 2R2 ) + ω ¯ 2 (T12 + R22 − S12 ) − S22 = 0.
(B.61)
By the intermediate value theorem the above equation has at least one positive root for each sign of (T12 + R22 − S12 ). This achieves the proof. Appendix B.9. Proof of Lemma 4 Let 2 2 2 2 F (ω) = |P1 | − |P0 | − |P2 | − 4|P0 |2 |P¯2 |2 .
There exists ω ˆ such that F (ˆ ω ) = 0. In fact, we have F (0) = S22 > 0. By lemma 3, there 2 exists ω ¯ such that |P1 | − |P0 |2 − |P2 |2 = 0. Hence, F (¯ ω ) = −4|P0 |2 |P¯2 |2 < 0 and ω ˆ ∈]0, ω ¯ [. Moreover, F (ω) is a polynomial and then, has a finite number of roots. This concludes the proof.
53
PII: DOI: Reference:
S0378-4754(19)30054-0 https://doi.org/10.1016/j.matcom.2019.02.006 MATCOM 4740
To appear in:
Mathematics and Computers in Simulation
Received date : 30 November 2017 Revised date : 26 December 2018 Accepted date : 11 February 2019 Please cite this article as: A. Mendy, J.J. Tewa, M. Lam et al., Hopf bifurcation in a grazing system with two delays, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.02.006 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript Click here to view linked References
Hopf bifurcation in a grazing system with two delays A. Mendya , J. J. Tewab , M. Lama , P. Tchinda Mouofob,1 a
University Cheikh Anta Diop, Department of Mathematics and Computer Science, Faculty of Science and Technic, P.O. Box 5005 Dakar b Department of Mathematics and Physics, National Advanced School of Engineering, University of Yaound´e I, P.O. Box 8390 Yaound´e, Cameroon
Abstract Life in the Sahelian zone is very difficult due to lack of resources on a permanent basis. This scarcity of resources leads the pastor, with their flock to organize themselves to better resist the difficult conditions. Pastoralism is of paramount importance for this region. There is an urgent need to understand this way of life in order to help pastors to better organize themselves. In this paper, we propose mathematical models that consist of three food trophic chains. We consider and formalize, in a given location, interactions between browsers, grazers and forage resources. We first analyze the resulting system without delays and in a second time take into account two positive and discrete time-delays. In our model, delays denote an average time required in order that food consumed by browsers and grazers are beneficial for them. Qualitative analysis of the delays-free model shows several equilibria as well as various situations of multistability. Considering delay as parameter, we investigate the effect of delay on the stability of equilibria. It is found that the delays can lead the system dynamic behavior to exhibit stability switches, and Hopf bifurcation occurs when the delay crosses some critical value. By applying the normal form theory and the center manifold theorem, the explicit formula which determine the stability and direction of the bifurcating periodic solutions are determined. Finally based on a nonstandard numerical scheme, we provide some numerical simulations in order to illustrate our qualitative results and support discussions. Keywords: Pastoralism, Browsers, Grazers, Hopf Bifurcation, Multistability. 1. Introduction Pastoralism describes the interdependent relationship between pastoralists, their ruminant herds and their habitat. The exploitation of livestock remains strongly dominated by extensive methods of herd management. The pastoral system is dominant along the Sahelian belt with relatively large herds of cattle and small ruminants (Report 2010 [26]). It is based on extensive exploitation of natural resources without recourse to zootechnical inputs, with the exception of years with critical forage deficit (see Report 2010 [26]). This type of livestock relies on the mobility or transhumance of herders and herds in search of water, pastures (during the long period of dry season, 8 to 10 months), and areas of salt cures (see also Report 2010 [26]). The dynamics of change in livestock activities are at the heart of Email addresses: [email protected] (A. Mendy), [email protected] (J. J. Tewa), [email protected] (M. Lam), [email protected] (P. Tchinda Mouofo) 1 Corresponding author: +(237) 679 44 12 29
Preprint submitted to Mathematics and Computers in Simulation
December 26, 2018
the challenges of sustainable development of many territories in the grassland and pastoral areas (Dedieu et al. 2009 [10]). The pastoral farming systems are those in which more than 90% of the dry matter consumed by livestock comes from grazing(Carriere 1996 [6]). Extensive pastoral production is practiced in 25% of the world’s land, from the arid zones of Africa (66% of the continent’s land) and the Arabian Peninsula to the tableland of Asia and Latin America (Nori et al. 2008 [24]). It provides 10% of the world’s meat production and supports some 200 million pastoralist households that raise almost 1 billion head of camels, cattle and small ruminants, about one third of them in sub-Saharan Africa, where pastoral farming represents about 20% of the Gross Domestic Product (GDP) of the different states (Nori et al. 2008 [24]). In the Sahel, the environmental consequences are less obvious (Smith 1996 [33]). The great majority of food consumed by ruminants in the sahelian region is still formed of natural pasture (Carriere 1996 [6], Research reports and studies 2009 [29]). Despite this important role, pastoral systems face many challenges related notably to the marginalization of pastoralists and to the increasing competition with other user groups especially farmers for accessing to natural resources (pasture, water points etc.) (see Report 2012 [27]). Recent studies indicate that the characteristics of the food and animal’s body size, rather than anatomical adaptations for grazing and browsing, have greater significance for the nutritional ecology of the herbivores (Shipley 1996 [30]). To better manage range livestock and habitats of free-ranging herbivores, ecologists and range managers often wish to estimate their forage intake rate. Morphological characteristics of grasses makes it easier to estimate the intake rate of grazers than that of browsers (Shipley 1996 [30]). The relative role of processes such as competition, facilitation and predation in structuring herbivore assemblages are inherently difficult to ascertain, and this has caused considerable controversy (McNaughton 1984 [23], Sinclair 1985 [31], Sinclair and Northon-Griffiths 1982 [32]). Inter-specific competition arises whenever one species reduces shared food resources below the level that can be exploited efficiently by another species (Bell 1971 [4]). However, by reducing plant biomass one species may also benefit another by facilitating access to forage of a suitable height or quality (Bell 1971 [4], McNaughton 1984 [23]). Long-term sustainable pastoral development requires a good knowledge of the dynamics of multiple factors underlining pastoralism, and here research has a crucial role to play. The links between research and policy also need to be strengthened, so that policy responds to the dynamics of pastoral livelihoods in an efficient manner (African Union 2010 [1]). The issue of livestock impacts on the dryland environment generally refers to a number of major themes, the main ones being desertification, drought, overgrazing, sedentarization of livestock breeders, rights of access to livestock resources (Carriere 1996 [6]). These highly mediatized subjects are the matter of intense controversy between specialists because, depending on whether an economic, sociological or ecological approach is adopted, the impacts are assessed in a different, even contradictory way (Carriere 1996 [6]). But it was also interesting to observe the evolution of pastoral areas according to the modalities of their exploitation by the animals. The demand for resources for herds depends on the number and type of animals and their nutritional requirements. The question of the resource is very important throughout the world and particularly in the Sahel, which has a wet season (june to october) and a dry season (november to may) (Tewa et al. 2014 [35]). An important issue in this part of Africa is the management of herds by pastoralists in order to maintain diversity between browsers and grazers in the sahel as well as manage the coexistence between resource and herbivores. 2
Mathematical models can help to better understand the action of pastoralists and to help them to better organize their work for an optimal yield. To authors knowledge, there is no mathematical model that addresses interactions between resource, browsers and grazers in Sahelian regions. In order to contribute to these problems related to Sahelian pastoralism, we therefore develop in this paper a model of interactions between the resource and a herd of herbivores in a given location. In our model, we make two important hypothesis. The first one is that as herbivory reduces the plant biomass of its maximum crop, the productivity of plants will increase. In fact, some empirical evidences show that when herbivory reduces plant biomass, for maintaining the ecosystem and avoid extinction, plants become resilient to herbivory and therefore increase their productivity (McNaughton 1984 [23], Lebon et al. 2014 [19]). The second hypothesis come from the fact that after consumption of resources, browsers and grazers need some time so that the consumed resource is beneficial to them. Taking the second hypothesis leads to delay differential systems and a rich dynamic of our model. This work is the first one which models interaction of herbivory and resource and taking time delays. The paper is organized as follows. We formulate a delayed mathematical model for herbivores and their resource interactions in section 2. In section 3, basics mathematical results are obtained and stability of equilibria are established when there is no delay in section 4. In section 5, the existence of Hopf bifurcation is established when delays are considered. Moreover, direction and stability of the Hopf bifurcation are also presented in section 5. Some numerical simulations that illustrate our analytical results are presented in section 6 and finally, a summary is presented in section 7. 2. Model construction 2.1. Dynamic of resource The resource includes grass, shrub and tree biomasses. The resource-herbivore interactions can be considered as a particular case of the well-mastered prey-predator system in which the logistic growth curve is most appropriate for writing primary productivity (NoyMeir 1975 [25], Caughley 1976 [7], Samuel et al. 1989 [38]). To understand the impact of herbivore pasture on the dynamics of the resource, several theories have been proposed including the Herbivores Optimization Hypothesis (HOH) (McNaughton 1984 [23], Hilbert et al. 1981 [17]). The HOH indicates that as the intensity of pasture increases, net primary productivity increases to optimum grazing intensity, then decreases to a lower level than the non-fattened resource (Samuel et al. 1989 [38]). This hypothesis of optimization of the grazing of herbivores was also used by Lebon et al. 2014 [19]. The resource growth function that models the HOH is defined by (see Lebon et al. 2014 [19]): ψ(B, H) = rB (H)B − δB 2 , with rB (H) =
H 1+ β
H 1− , γ
where B represents plant biomass and H herbivore (or pest) population. In this paper, we use this approach to model the growth of the resource in the presence of two populations of herbivores. We consider the following HOH-like function ψ(R, Hb , Hg ) = rφR (Hb , Hg )R − δR2 , 3
with φR (Hb , Hg ) =
αb Hb + αg Hg 1+ β
αb Hb + αg Hg 1− , γ
where αb , αg ∈ {0; 1}, γ > β and αb Hb + αg Hg < γ. The term δR2 is the density-dependent death rate of the resource as a consequence of local competition for water and nutrients. We note that in the herbivore-free case (i.e. ψ = 1), the resource follows a logistic growth (e.g. Tewa et al. 2014 [35], Zhilan et al. 2011 [11], Li et al. 2006 [20]). Primary production is heavily dependent on precipitation (Carriere 1996 [6]), which means that the intrinsic growth rate may be negative. The intrinsic growth rate has been defined as (see Tewa et al. 2014 [35]): r = r01 p + r02 − r2 , where p denotes the mean annual precipitation, r01 is the growth rate of the resource with precipitation, r02 is the growth rate of the resource without precipitation, r2 is the natural death rate. The particularity of this model, which is the difference with classic predator-prey models in the literature, is the fact that the intrinsic production rate can be negative (r < 0) (Tewa et al. 2014 [35]). This situation is frequently observed in Sub-Saharan Africa, particularly in Sahelian region when precipitations are scarce (Tewa et al. 2014 [35]). Intensification of this situation can conduct pastoralists to migrate with their herds: this is the transhumance phenomenon. 2.2. Dynamic of browsers Browsers include camels, beefs, giraffes, horses etc. (Johan et al. 1998 [8]). They feed on grass, shrubs and stress biomass. In the study of interactions between herbivores, and resource, it is important to determine which specific form of functional response, describing the amount of resource consumed by an animal in the herd per unit of time, is ecologically meaningful and provides a solid basis for theoretical development (Tewa et al. 2013 [37], Tewa et al. 2012 [36]). In the simplest case, such a function is a linear function of the forage resources density, which is integrated into the classic Lotka-Volterra model. The linear functional response is a limited case and can only be seen over a short period. However, the Holling type II or Holling type III functional response can be used as an improvement of the linear-like function. With the Holling type II function, at low resource biomass, resource experiences herbivory while with the Holling type III function herbivory is very low at low resource biomass. The Holling type III function is typical of generalist consumers moving from one resource to another and concentrating their activities in a region where the resource is abundant (Bates et al. 2016 [3]). This function is defined by: P (R) =
aR2 , b + R2
where a is the maximum consumption rate of browsers on the resource, b is the half-saturation parameter. As modelling hypothesis, we assume that there is no alternative food for herbivores. Consequently, in absence of resource biomass, both browsers and grazers go to extinction.
4
2.3. Dynamic of grazers Grazers include sheep, goats, etc. (Koppel et al. 1998 [8]). Grazers feed on grass and can not reach shrubs and trees. The interactions between the resource and grazers follow the functional response that takes into account the interspecific competition (competition between browsers and grazers). An interspecific competition is a competition between two different species for the same resource. The competition for the resource arises from the fact that the resource is limited and that the competitors interfere to obtain it. BeddingtonDeAngelis 1975 [9] proposed a function defined by: g(R, Hg , Hb ) =
c1 Hg R, 1 + c2 Hg + c3 Hb
where c1 , c2 and c3 are positive parameters. c1 is the catch rate, c2 is the manipulation time, and c3 is the interference amplitude between browsers and grazers. The term c1 Hg corresponds to the consumption of the resource by the small herbivores 1 + c2 Hg + c3 Hb taking into account its competition with the large herbivores (DeAngelis et al. 1975 [9]). The functional response of Beddington-DeAngelis suggests that the effects of interferences between herbivores become negligible when the density of the resource biomass is high. This function characterizes the fact that in addition to encountering the resource, herbivores likely collide with each other because the resource is limited. This function also assumes that herbivores spend their time in three activities: searching for the resource, consuming the resource and manipulating other herbivores. After consumption of the resource, browsers need a time so that the consumed resource is beneficial to them. This period is represented by τ1 . Similarly, for the consumption of the resource by grazers, this period is represented by τ2 . The following system of ordinary differential equations models interactions between forage resources, browsers and grazers : dR(t) dt
αb Hb (t) + αg Hg (t) αb Hb (t) + αg Hg (t) = r 1+ 1− R(t) − δR2 (t) β γ −
c1 Hg (t) aR2 (t) R(t) − Hb (t), 1 + c2 Hg (t) + c3 Hb (t) b + R2 (t)
dHb (t) aR2 (t − τ1 ) = e Hb (t) − (db + hb )Hb (t), b dt b + R2 (t − τ1 ) dHg (t) c1 Hg (t) = eg R(t − τ2 ) − (dg + hg )Hg (t), dt 1 + c2 Hg (t) + c3 Hb (t)
(1)
where, db and dg are the natural mortality rates of large herbivores and small herbivores respectively, hb and hg are the harvest rates of large herbivores and small herbivores respectively, eb and eg are the conversion rates of the resource consumed into large and small herbivorous biomass, respectively. System (1) is considered with the following initial values R(θ) = ϕ1 (θ) ≥ 0, Hb (θ) = ϕ2 (θ) ≥ 0, Hb (θ) = ϕ3 (θ) ≥ 0, θ ∈ [−τ, 0], where τ = max{τ1 , τ2 }, ϕ1 , ϕ2 , ϕ3 ∈ C [−τ, 0], R3+ , the Banach space of continuous functions mapping the interval [−τ, 0] into R3+ , with R3+ = {(R, Hb , Hg )/ R ≥ 0, Hb ≥ 0, Hg ≥ 0}. 5
3. Basic results 3.1. Boundedness of solutions We start by showing that solutions of System (1) that start in R3+ will remain there and are uniformly bounded. The following Lemma is valid. Lemma 1. 1) If r ≤ 0, then all the trajectories of the System (1) go to the zero solution. 2) If r > 0, (i) The nonnegative orthant R3+ is a positively invariant region for the System (1). (ii) Let ε > 0, the set Ωε defined by r γ 3 Ωε = (R, Hb , Hg ) ∈ R+ : 0 ≤ R ≤ 1+ + ε; 0 ≤ αb Hb + αg Hg ≤ γ δ β is a positively invariant and absorbing set for the System (1). Proof: See Appendix B.1. From Lemma 1, one deduces that System (1) is biologically well-posed, that is, solutions are non negative for non negative initial data and remain bounded. Moreover, we also deduce from Lemma 1 that System (1) is a dissipative system (Hale 1988 [15]). In the sequel we give the results concerning existence of equilibria for System (1). 3.2. Equilibria of the model We reach the step of characterizing equilibria of System (1). Let us set, ˜ ˜ rαb δR eb a β aR 4βγηa , ηs = , ηa = 1− , S01 = , , σR˜ = N0 = 2 2 ˜ βσR˜ γ r µb b+R 1 2 (γ − β) 1 − + 4βγ ηs δµg βγc2 + αg (γ − β) γ−β L0 = , L1 = and L2 = . 2βγδµg c2 βγc1 αg βγδµg c2 reg c1 + + reg c1 r c2 rαg eg c1 Proposition 1 (Existence of equilibria). The following results hold for System (1). r (1) Equilibrium E0 = (0, 0, 0) always exists while equilibrium E1 = , 0, 0 exists whenδ ever r > 0. ˜ H ˜ b , 0 of System (1), there exist two cases. (2) For the grazer-free equilibrium E˜ = R, ˜ H ˜ 0 , 0) with (a) When N0 = 1, we have a unique equilibrium E˜0 = (R, b s bµb ˜0 = 1 γ − β 1 − 1 . ˜= and H R b e b a − µb 2αb ηs E˜0 is ecologically meaningful whenever S01 > 1 and ηs > 1.
6
˜ H ˜ 1 , 0) and E˜2 = (R, ˜ H ˜ 2 , 0) with (b) When N0 < 1, we have two equilibria: E˜1 = (R, b b s 2 ˜ b1 = 1 (γ − β) 1 − 1 − (γ − β)2 1 − 1 H − 4βγ(ηa − 1) 2αb ηs ηs and
˜ b2 = 1 (γ − β) 1 − 1 H 2αb ηs
+
2 1 (γ − β)2 1 − − 4βγ(ηa − 1) . ηs
s
E˜1 is ecologically meaningful whenever N0 < 1, S01 > 1, ηa > 1 and ηs > 1 while E˜2 is ecologically meaningful whenever (S01 > 1 and ηa < 1) or (S01 > 1, ηa < 1 and ηs < 1). ¯ 0, H ¯ g ) of System (1) exists with (3) The browser-free equilibrium E¯ = (R, µ g ¯= ¯ g ) and H ¯ g is a positive solution of equation R (1 + c2 H eg c1 ¯ 3 + γ − β (1 − L2 ) H ¯ 2 + βγc2 + αg (γ − β) (1 − L1 ) H ¯ g + βγ (L0 − 1) = 0. H g g αg L 2 αg c2 L1 αg2 c2
(2)
(4) A coexistence equilibrium E ∗ = (R∗ , Hb∗ , Hg∗ ) exists if S01 > 1, N ∗ > 1, and Hg∗ being a positive solution of equation D2 Hg∗2 + D1 Hg∗ + D0 = 0,
(3)
r
eg c1 R∗ − µg (1 + c2 Hg∗ ) bµb , Hb∗ = and the expression of D2 , D1 e b a − µb c3 µ g and D0 are given in the proof. ∗
where, R =
Proof: See Appendix B.2. The complexity of results given in Proposition 1 calls for some remarks: first, concerning existence of browsers-free equilibrium, one may obtain up to two equilibria of that type thanks to the modelling option that consists on considering the HOH and BeddingtonDeAngelis’s function. A similar observation can be done for the grazer-free equilibrium as well as the coexistence equilibrium. Indeed, thank to the modelling choice of HOH and Holling type III functional response, on may obtain up to three grazer-free equilibria and up to two coexistence equilibria. This is not a common feature of predator-prey models. However in our case, it is desirable since it may lead to several cases of multistability i.e. more than one equilibrium is stable in the same time (see below). In the sequel, we deal with the problem of stability characterization of various equilibria of model (1). For that, we first consider the case where there is no delay and, secondly focus on the case where delays are considered. This process is motivated by the fact that if an equilibrium is unstable without delay, it will remain unstable when delays are considered (see Rebecca et al. 2000 [28] and Martin et al. 2001 [22]).
7
4. Stability analysis when τ1 = τ2 = 0 Let us set, ˜ eb a R˜ H˜b e g c1 R eb ar2 δµg 1 = = S0 = ; L = ; S ; S ; 0 0 0 µb (bδ 2 + r2 ) e g c1 r µb µg (1 + c3 H˜b ) ˜2H ˜b ¯2 aR eb aR ˜ H˜b ¯ R = ; S S0R = . ¯2) 1 ˜ b + δ(b + R ˜ 2 )2 µb (b + R abH Theorem 1 (Local stability results). The following results hold. (1) If r < 0, then E0 = 0, 0, 0 is locally asymptotically stable. r , 0, 0 is locally asymptotically stable. (2) If S0 < 1 and L0 > 1, then E1 = δ ˜ ˜ ˜ ˜ (3) Suppose that S0RHb < 1 and S1RHb < 1. We have two cases. ˜ ˜ ˜ (a) Grazer-free equilibrium E = R, Hb , 0 is locally asymptotically stable whenever it is unique. (b) If there exist two grazer-free equilibria, then one is locally asymptotically stable and the other is unstable. ¯ R ¯ (4) Suppose that S0 < 1, then we have two cases for the browser-free equilibria E = ¯ 0, H ¯g . R, (a) If we have one equilibrium, then it is locally asymptotically stable.
(b) If we have three equilibria, then two are locally asymptotically stable and the other is unstable. Proof: See Appendix B.3. Before going further, let us stress the fact that in Theorem 1 the case r = 0 is not considered. However, when r = 0, using a comparison argument (e.g. Yatat et al. 2014 [39]), solution of System (1) go to extinction. Furthermore, from Theorem 1, one deduces that multistability involving several equilibria of the same type or of different types may exist. This means that depending on initial conditions, the system may converge to different equilibria. From an ecological point of view, these multistabilities mean that transitions between herd, browsers and grazers amount may be more complex and confirm that these systems in Sahelian are very complex systems. In the sequel we provide some interpretations of Theorem 1 in terms of pastoral mechanisms. Pastoral interpretation
• The stability of equilibrium E0 = 0, 0, 0 means that there will be no more forage resources and no animals of the herd in the considered area. There are many explanations of this situation. Droughts can cause animals mortalities through starvation, emergency slaughtering, and sales, or definitive herds migrations (transhumance), which can create severe drops in the herd sizes. • When this removal concerns only herbivores, depending on some climatic changes, the forage resources growth towards the maximal quantity needed for livestock and the r equilibrium E1 = , 0, 0 is stable. δ 8
• If the disappearance concerns only grazers, we have a coexistence between the resource and the browsers. A situation that can be explained by the extinction of the grass resource and therefore the small herbivores may disappear or the pastor may decide to migrate with the grazers (transhumance). Consequently, the grazers-free equilibrium ˜ H ˜ b , 0 is stable. E˜ = R, • If the disappearance concerns only the browsers, we have a coexistence between the resource and the grazers. A situation that can be explained by the insufficiency of the available resource because browsers need a large availability of resources. Consequently, this situation can cause pastoralists to migrate with browsers and therefore, ¯ 0, H ¯ g is sable. the browsers-free equilibrium E¯ = R,
Theorem 2 (Hopf bifurcation). ˜ H ˜ b , 0) when δ The System (1) has a Hopf bifurcation around grazer-free equilibria E˜ = (R, passes through δh , where ˜ 2 − b) aH˜b (R . (4) δh = ˜ 2 )2 (b + R Proof: See Appendix B.4. The pastoral interpretation of the Hopf bifurcation is that herbivores will coexist with the resource, exhibiting oscillatory balance behavior. Theorem 3 (Local stability of coexistence equilibria). The coexistence equilibrium E ∗ = (R∗ , Hb∗ , Hg∗ ) is locally asymptotically stable if the conditions of the Routh-Hurwitz criterion are verified, that is: C1 < 0,
C2 < 0
and
C3 > 0,
(5)
where C1 = m∗11 + m∗33 , C2 = (m∗11 + m∗33 )(m∗11 m∗33 − m∗13 m∗31 − m∗12 m∗21 ) + (m∗12 m∗21 m∗33 − m∗13 m∗21 m∗32 ), C3 = m∗12 m∗21 m∗33 − m∗13 m∗21 m∗32 , with m∗11 , m∗33 , m∗12 , m∗21 , m∗13 and m∗31 as given in the proof. Proof: See Appendix B.5. 5. Stability analysis with positive delays In this section we study the stability behavior of equilibria in the presence of discrete delays τ1 > 0 and τ2 > 0. Let E ∗ = (R∗ , Hb∗ , Hg∗ ) be an equilibrium of the model, and let us define u(t) = R(t) − R∗ , v(t) = Hb (t) − Hb∗ and w(t) = Hg (t) − Hg∗ . Then the linearized system of (1) at E ∗ is given by
9
du(t) = a1 u(t) + a2 v(t) + a3 w(t); dt dv(t) (6) = m1 u(t − τ1 ) + b2 v(t); dt dw(t) = m u(t − τ ) + c v(t) + c w(t); 2 2 2 3 dt where αb Hb∗ + αg Hg∗ c1 Hg∗ αb Hb∗ + αg Hg∗ 2abR∗ 1− −2δR∗ − − a1 = r 1 + Hb∗ ; ∗ ∗ ∗2 2 β γ 1 + c2 Hg + c3 Hb (b + R ) ∗2 c1 c3 Hg∗ rαb aR a2 = γ − β − 2(αb Hb∗ + αg Hg∗ ) R∗ − R∗ ; + βγ b + R∗2 (1 + c2 Hg∗ + c3 Hb∗ )2 c1 R∗ (1 + c3 Hb∗ ) rαb 2eb abR∗ ∗ ∗ a3 = γ − β − 2(αb Hb∗ + αg Hg∗ ) R∗ − = ; m H ; 31 βγ (1 + c2 Hg∗ + c3 Hb∗ )2 (b + R∗2 )2 b eg c1 Hg∗ 2eb abR∗ ∗ eb aR2 m1 = H ; b2 = − µb ; m2 = ; (b + R∗2 )2 b b + R2 1 + c2 Hg∗ + c3 Hb∗ eg c1 c3 Hg∗ R∗ µg c2 Hg∗ c2 = − ; c3 = − . (1 + c2 Hg∗ + c3 Hb∗ )2 1 + c2 Hg∗ + c3 Hb∗ We then express system (6) in matrix form as follow u(t) u(t) u(t − τ1 ) u(t − τ2 ) d v(t) = A v(t) + A1 v(t − τ1 ) + A2 v(t − τ2 ) , dt w(t) w(t) w(t − τ1 ) w(t − τ2 ) where
a1 a2 a3 0 0 0 A = 0 b2 0 , A1 = m1 0 0 , 0 c2 c3 0 0 0
0 0 0 A2 = 0 0 0 . m2 0 0
The characteristic equation of system (6) at equilibrium E ∗ is given by ∆(ξ, τ1 , τ2 ) = |ξI − A − A1 e−ξτ1 − A2 e−ξτ2 | = 0, that is, P0 (ξ) + P1 (ξ)e−ξτ1 + P2 (ξ)e−ξτ2 = 0,
(7)
where, P0 (ξ) = ξ 3 − (a1 + c3 + b2 )ξ 2 + (a1 c3 + a1 b2 + c3 b2 )ξ − a1 c3 b2 , P1 (ξ) = −m1 a2 ξ + m1 a2 c3 − m1 a3 c2 and P2 (ξ) = −a3 m2 ξ + a3 m2 b2 . 5.1. Stability analysis when τ1 > 0 and τ2 = 0
0, 0, 0 ,
In the case where τ1 > 0 and τ2 = 0, the characteristic equations at E0 = r ¯ 0, H ¯ g are independent of the delay τ1 . Then equilibria E0 , E1 , 0, 0 and E¯ = R, E1 = δ and E¯ have the same stability as in the delays-free case.
10
˜ existence of a switch threshold 5.1.1. Stability analysis of E: ˜ When τ1 > 0 and τ2 = 0, the characteristic equation In the sequel, we only focus on E. ˜ at E is given by (8) ξ 3 + B1 ξ 2 + B2 ξ + (B3 ξ + B4 )e−ξτ1 = 0. where B1 = −(˜ a1 + c˜3 ), B2 = a ˜1 c˜3 , B3 = −m ˜ 1a ˜2 and B4 = a ˜2 c˜3 m ˜ 1. Substituting ξ = iω in the characteristic equation (8), the real and the imaginary part satisfy ω 3 − B2 ω = B3 ω cos ωτ1 − B4 sin ωτ1 ,
(9)
B1 ω 2 = B3 ω sin ωτ1 + B4 cos ωτ1 .
(10)
Eliminating τ1 by squaring and adding equations (9) and (10), we get the algebraic equation ω 6 + (B12 − 2B2 )ω 4 + (B22 − B32 )ω 2 − B42 = 0.
(11)
Substituting ω 2 = η in the above equation, we obtain the following cubic equation η 3 + θ20 η 2 + θ10 η + θ00 = 0,
(12)
where, θ20 = B12 − 2B2 , θ10 = B22 − B32 and θ00 = −B42 . Note that θ20 = B12 − 2B2 = (˜ a1 + c˜3 )2 − 2˜ a1 c˜3 = a ˜21 + c˜23 > 0 and θ00 = −B42 < 0. By Descartes’rule of sign, the equation (12) has at least one positive root. Remark 1. When τ1 = τ2 = 0, by Theorem 1 (3), E˜ is locally asymptotically stable if ˜ ˜ ˜ ˜ ˜1 < 0) and a ˜2 < 0. S0RHb < 1 (i.e. c˜3 < 0), S1RHb < 1 (i.e. a The following theorem gives a criterion for the switching in the stability behavior of E˜ in terms of the delay parameter τ1 . Theorem 4. Suppose that E˜ exists and is locally asymptotically stable for (1) with τ1 = τ2 = 0. Also let, ω02 = η0 be a positive root of (12). Then, there exists τ1 = τ10 such that E˜ is locally asymptotically stable for τ1 ∈ (0, τ10 ] and unstable for τ1 > τ10 . Furthermore, the system undergoes a Hopf bifurcation at E˜ when τ1 = τ10 . Proof: See Appendix B.6. ˜ direction and stability of the Hopf Bifurcation 5.1.2. Stability analysis of E: ˜ H ˜ b , 0) of system In the previous section, the stability of the positive equilibrium E˜ = (R, (1) has been analized as well as the existence of Hopf bifurcation at this equilibrium. In this section, we shall study the properties of the Hopf bifurcation obtained in Theorem 4 and the stability of bifurcated periodic solutions occurring through Hopf bifurcation by using the normal form theory and the center manifold reduction for retarded functional differential equations (RFDEs) due to Hassard et al. 1981 [16]. Throughout this section, we always assume that system (1) undergoes Hopf bifurcation at the positive equilibrium 0 ˜ H ˜ b , 0) for τ1 = τ1j E˜ = (R, , (j = 0, 1, 2....) and then ±iω0 is the corresponding purely imaginary roots of the characteristic equation (8). Theorem 5. There exist three real numbers µ2 , β2 and T2 which determine the qualities of bifurcating periodic solution in the center manifold at the critical value τ10 such that: 11
(a) When µ2 > 0, then the Hopf bifurcation is supercritical and the bifurcating periodic solution exists for τ1 > τ10 . When µ2 < 0, then the Hopf bifurcation is subcritical and the bifurcating periodic solution exists for τ1 < τ10 . (b) When β2 < 0, the bifurcating periodic solutions are stable and when β2 > 0, the bifurcating periodic solutions are unstable. (c) The period increases when T2 > 0, and decreases when T2 < 0. Proof: See Appendix B.7. 5.1.3. Stability analysis of E ∗ when τ1 > 0 and τ2 = 0 When τ1 > 0 and τ2 = 0, the characteristic equation at E ∗ is given by ∆(ξ, τ1 ) = M (ξ) + N (ξ)e−ξτ1 = 0,
(13)
where M (ξ) = ξ 3 + A2 ξ 2 + A1 ξ and N (ξ) = S1 ξ + S0 with A2 = −(a∗1 + c∗3 ), A1 = a∗1 c∗3 − a∗3 m∗2 , S1 = m∗1 a∗2 and S0 = m∗1 a∗2 c∗3 − m∗1 a∗3 c∗2 . Substituting ξ = iω in the characteristic equation (13), the real and imaginary parts satisfy ω 3 − A1 ω = S1 ω cos ωτ1 − S0 sin ωτ1 ,
(14)
A2 ω 2 = S1 ω sin ωτ1 + S0 cos ωτ1 .
(15)
Eliminating τ1 by squaring and adding equations (14) and (15), gives the algebraic equation ω 6 + (A22 − 2A1 )ω 4 + (A21 − S12 )ω 2 − S02 = 0.
(16)
Substituting ω 2 = η in the above equation, gives the following cubic equation η 3 + (A22 − 2A1 )η 2 + (A21 − S12 )η − S02 = 0.
(17)
By the intermediate value theorem, the equation (17) has at least one positive root. The following theorem gives a criterion for the switching in the stability behavior of E ∗ in terms of the delay parameter τ1 . Theorem 6. Suppose that E ∗ exists and is locally asymptotically stable for (1) with τ1 = τ2 = 0. Also suppose that ω12 = η1 is a positive root of (17). Then, there exists τ1 = τ100 such that E ∗ is locally asymptotically stable for τ1 ∈ (0, τ100 ] and unstable for τ1 > τ100 . Furthermore, the system undergoes a Hopf bifurcation at E ∗ when τ1 = τ100 . The proof of Theorem 6 is done in the same way as that of Theorem 4. 5.1.4. Stability analysis of E ∗ : direction and stability of the Hopf Bifurcation Similar to section 5.1.2, the direction of Hopf of E ∗ is given on the sequel Theorem 7. There exist three real numbers µ3 , β3 and T3 which determine the quality of the bifurcating periodic solution in the center manifold at the critical value τ100 such that: (a) When µ3 > 0, then the Hopf bifurcation is supercritical and the bifurcating periodic solution exists for τ1 > τ100 . When µ3 < 0, then the Hopf bifurcation is subcritical and the bifurcating periodic solution exists for τ1 < τ100 . (b) The bifurcating periodic solutions are stable when β3 < 0, and unstable if β3 > 0 (c) The period increases when T3 > 0 and decreases when T3 < 0. The proof of the Theorem 7 is done in the same way as that of Theorem 5. 12
5.2. Stability analysis when τ1 = 0 and τ2 > 0
0, 0, 0 ,
In the case where τ1 = 0 and τ2 > 0, the characteristic equations at E0 = r ˜ 0, H ˜ g are independent of the delay τ2 . Then, equilibria E0 , E1 E1 = , 0, 0 and E˜ = R, δ and E˜ have the same stability properties as in the delays-free cases. 5.2.1. Analysis of E¯ when τ1 = 0 and τ2 > 0 When τ1 = 0 and τ2 > 0, the characteristic equation at E¯ is given by ξ 3 + C1 ξ 2 + C2 ξ + C3 + (C4 ξ + C5 )e−ξτ2 = 0.
(18)
where C1 = −(¯ a1 + c¯3 + ¯b2 ); C2 = a ¯1 c¯3 + ¯b2 a ¯1 + ¯b2 c¯3 ; C3 = −¯ a1 c¯3¯b2 ; C4 = −¯ a3 m ¯ 2 and C5 = a ¯3 m ¯ 2¯b2 . Substituting ξ = iω in the characteristic equation (18), the real and imaginary parts satisfy ω 3 − C2 ω = C4 ω cos ωτ2 − C5 sin ωτ2 ,
(19)
C1 ω 2 − C3 = C4 ω sin ωτ2 + C5 cos ωτ2 .
(20)
Eliminating τ2 by squaring and adding equations (19) and (20), gives the algebraic equation ω 6 + (C12 − 2C2 )ω 4 + (C22 − C42 − 2C3 C1 )ω 2 + C32 − C52 = 0.
(21)
Substituting ω 2 = η in the above equation, gives the following cubic equation η 3 + θ2 η 2 + θ1 η + θ0 = 0,
(22)
where, θ2 = C12 − 2C2 ; θ1 = C22 − C42 − 2C3 C1 and θ0 = C32 − C52 . We have θ2 = C12 − 2C2 = (¯ a1 + c¯3 + ¯b2 )2 − 2(¯ a1 c¯3 + ¯b2 a ¯1 + ¯b2 c¯3 ) = a ¯21 + c¯23 + ¯b22 > 0 and θ0 = C32 − C52 = (C3 − C5 )(C3 + C5 ) = (−¯ a1 c¯3¯b2 − a ¯3 m ¯ 2¯b2 )(−¯ a1 c¯3¯b2 + a ¯3 m ¯ 2¯b2 ) = ¯b22 (¯ a1 c¯3 + a ¯3 m ¯ 2 )(¯ a1 c¯3 − a ¯3 m ¯ 2 ). We consider only the case when (22) possesses exactly one positive root. This occurs if ¯1 c¯3 + a ¯3 m ¯ 2 < 0. Under this condition, the equation (22) possesses a θ0 = C32 − C52 < 0 i.e. a 2 positive root η2 = ω2 and thus (18) has the pair of purely imaginary roots ±iω2 . Remark 2. When τ1 = τ2 = 0, by Theorem 1 (4) E¯ is locally asymptotically stable, if ¯ ¯3 < 1. S0R < 1 and a Using a similar approach as in Theorem 4, one deduces the following result. Theorem 8. Suppose that E¯ exists and is locally asymptotically stable for (1) when τ1 = τ2 = 0. Also suppose that, ω22 = η2 is a positive root of (22). Then, there exists a τ2 = τ20 such that E¯ is locally asymptotically stable for τ2 ∈ (0, τ20 ] and unstable for τ2 > τ20 . Furthermore, the system undergoes a Hopf bifurcation at E¯ when τ2 = τ20 . The proof of Theorem 8 is done in the same way as that of Theorem 4.
13
5.2.2. Nonexistence of switching stability around E¯ When τ1 = 0 and τ2 > 0, the characteristic equation at E¯ is given by ξ 3 + C1 ξ 2 + C2 ξ + C3 + (C4 ξ + C5 )e−ξτ2 = 0.
(23)
Equation (23) can be written as P (ξ) + Q(ξ)e−ξτ2 = 0,
(24)
where, P (ξ) = ξ 3 + C1 ξ 2 + C2 ξ + C3 and Q(ξ) = C4 ξ + C5 . Clearly P (ξ) and Q(ξ) are both analytic functions in Re (ξ) > 0. The following results hold: i) P (0) + Q(0) = C3 + C5 > 0, ¯ ii) P (−iy) = P¯ (iy) and Q(−iy) = Q(iy), Q(ξ) iii) lim sup = 0. |ξ|→+∞ P (ξ)
In addition
F (y) = |P (−iy)|2 − |Q(iy)|2 = y 6 + (C12 − 2C2 )y 4 + (C22 − C42 − 2C3 C1 )y 2 + C32 − C52 , is a cubic expression in y 2 . By using simple algebraic calculation we have, C12 − 2C2 = a ¯21 + c¯23 + ¯b22 > 0, C22 − C42 − 2C3 C1 = a ¯21¯b22 + ¯b22 c¯23 + (¯ a1 c¯3 + a ¯3 m ¯ 2 )(¯ a1 c¯3 − a ¯3 m ¯ 2 ), 2 2 2 C3 − C5 = ¯b2 (¯ a1 c¯3 + a ¯3 m ¯ 2 )(¯ a1 c¯3 − a ¯3 m ¯ 2 ), ¯1 c¯3 + a ¯3 m ¯ 2 > 0. C32 − C52 > 0 and C22 − C42 − 2C3 C1 > 0 if a Therefore F (y) = 0, has no positive real root if a ¯1 c¯3 +¯ a3 m ¯ 2 > 0. Hence by applying Theorem 4.1 in Kuang 1993 [18], the system (1) does not possess any stability switches around E¯ when a ¯1 c¯3 + a ¯3 m ¯ 2 > 0. Remark 3. Since the system (1) does not always undergo a Hopf bifurcation at the positive ¯ we are not interested by the direction of the Hopf bifurcation at equilibrium equilibrium E, E¯ 5.2.3. Stability analysis of E ∗ when τ1 = 0 and τ2 > 0 When τ1 = 0 and τ2 > 0, the characteristic equation at E ∗ is given by P0 (ξ) + Q0 (ξ)e−ξτ2 = 0,
(25)
where, P0 (ξ) = ξ 3 + R2 ξ 2 + R1 ξ + R0 and Q0 (ξ) = T ξ with R2 = −(a∗1 + b∗2 + c∗3 ), R1 = a∗1 b∗2 + c∗3 a∗1 + b∗2 c∗3 − a∗2 m∗1 , R0 = c∗3 a∗2 m∗1 − m∗1 a∗3 c∗2 and T = −a∗3 m∗2 . Clearly P0 (ξ) and Q0 (ξ) are both analytic functions in Re (ξ) > 0. The following results hold: i) P0 (0) + Q0 (0) = R0 6= 0, ¯ 0 (ix), ii) P0 (−ix) = P¯0 (ix) and Q0 (−ix) = Q 14
Q (ξ) 0 iii) lim sup = 0. |ξ|→+∞ P0 (ξ)
Then,
F (x) = |P0 (−ix)|2 − |Q0 (ix)|2 = x6 + (R22 − 2R1 )x4 + (R12 − T 2 − 2R0 R2 )x2 + R02 = 0, is a cubic expression in x2 . Therefore, by Descartes’rule of sign, the equation F (w) = 0 has no positive real root if R22 − 2R1 > 0 and R12 − T 2 − 2R0 R2 > 0.
(26)
By applying Theorem 4.1 in Kuang 1993 [18], when τ1 = 0, τ2 > 0 and conditions (26) hold, the system does not possess any stability switches around E ∗ . 5.3. Stability analysis when τ1 > 0 and τ2 > 0 Before going further, let us recall the following results. (a) When τ1 > 0 and τ2 > 0, the characteristic equation at E¯ can be written as P0 (ξ) + P1 (ξ) + P2 (ξ)e−ξτ2 = 0, where, P0 , P1 and P2 are some polynomials. Therefore, the delay τ1 has no influence ¯ on equilibrium E. (b) When τ1 > 0 and τ2 > 0, the characteristic equation at E˜ can be written as P0 (ξ) + P2 (ξ) + P1 (ξ)e−ξτ1 = 0, where, P0 , P1 and P2 are some polynomials. Therefore, the delay τ2 has no influence ˜ on equilibrium E. In what follows, we use the algebraic method developed in Lin et al. 2012 [21] to derive a closed form for stability switching curves of our delayed system when τ1 > 0 and τ2 > 0. We provide some properties of these curves and stability switching direction. To applied this method, let us write the following characteristic equation at equilibrium E ∗ , D(ξ, τ1 , τ2 ) = P0 (ξ) + P1 (ξ)e−ξτ1 + P2 (ξ)e−ξτ2 = 0,
(27)
where, P0 = ξ 3 + R1 ξ 2 + R2 ξ, P1 = S1 ξ + S2 and P2 = T1 ξ, with R1 , R2 , S1 , S2 and T1 given by equation (7). A- Stability crossing curves Lemma 2. As (τ1 , τ2 ) varies continuously in R2+ , the number of zeros (with multiplicity counted) of D(ξ, τ1 , τ2 ) on C+ can change only if a zero appears on or cross the imaginary axis. The proof of Lemma 2 can be found in any book on functional differential equations, for example in Kuang 1993 [18] and Smith 2010 [34]. Remark 4. From the above lemma, to study stability switch, we need to find the roots on the imaginary axis. 15
In the following, we investigate the existence of purely imaginary roots ξ = iω (w > 0) of equation (27). Substituting this into (27), gives −iωτ2 + P1 (iω)e−iωτ1 = 0. D(iω, τ1 , τ2 ) = P0 (iω) + P2 (iω)e Since |e−iωτ1 | = 1, we have
|P0 (iω) + P2 (iω)e−iωτ2 | = |P1 (iω)|, which is equivalent to
P 0 + P2 e
−iωτ2
(28)
iωτ2 ¯ ¯ = |P1 |2 P0 + P2 e
where P0 = −R1 ω 2 + i(−ω 3 + R2 ω), P1 = iS1 ω + S2 and P2 = iT1 ω. After simplification, (28) becomes |P0 |2 + |P2 |2 + 2Re (P0 P¯2 ) cos(ωτ2 ) − 2Im (P0 P¯2 ) sin(ωτ2 ) = |P1 |2 . Thus,
|P1 |2 − |P0 |2 − |P2 |2 = 2Re (P0 P¯2 ) cos(ωτ2 ) − 2Im (P0 P¯2 ) sin(ωτ2 ). (29) Since P0 6= 0 and P2 6= 0, we have [Re (P0 P¯2 )]2 + [Im (P0 P¯2 )]2 = |P0 P¯2 |2 > 0. Therefore, equation (29) can be written as |P1 |2 − |P0 |2 − |P2 |2 Re (P0 P¯2 ) Im (P0 P¯2 ) = cos(ωτ ) − sin(ωτ2 ). 2 2|P0 P¯2 | |P0 P¯2 | |P0 P¯2 |
(30)
Then, there exists a continuous function φ2 (ω) such that cos(φ2 (ω)) =
Re (P0 P¯2 ) Im (P0 P¯2 ) and sin(φ (ω)) = . 2 |P0 P¯2 | |P0 P¯2 |
Therefore, equation (30) becomes |P1 |2 − |P0 |2 − |P2 |2 = cos(φ2 (ω) + ωτ2 ). 2|P0 P¯2 |
(31)
Obviously, a sufficient and necessary condition for the existence of τ2 ∈ R+ satisfying equation (31) is 2 2 2 (32) |P1 | − |P0 | − |P2 | ≤ 2|P0 P¯2 |. Let us consider
n o Ω2 = ω ∈ R+ : |P1 |2 − |P0 |2 − |P2 |2 ≤ 2|P0 P¯2 | .
We have the following result.
Lemma 3. The set Ω2 is non empty.
16
Proof: See Appendix B.8. Let |P1 |2 − |P0 |2 − |P2 |2 ˜ , ψ2 ∈ [0, π]. 2|P0 P¯2 |
cos(ψ˜2 (ω)) = We obtain
ψ˜2 (ω) − φ2 (ω) + 2πn2 , n2 ∈ Z. (33) ω Substituting (33) into (27), gives the expression of τ1 (ω) unconditionally with each ω ∈ Ω2 i.e. P1 1 + 2πn1 , n1 ∈ Z. (34) τ1,n1 (ω) = arg − ω P0 + P2 e−iωτ2 Thus, the stability crossing curves are I = (τ1 (ω), τ2 (ω)) ∈ R2+ , ω ∈ Ω2 . (35) τ2,n2 (ω) =
Another way to find τ1 is to analyze τ1 similarly to the analysis of τ2 , which gives τ1,n1 (ω) =
ψ˜1 (ω) − φ1 (ω) + 2πn1 , n1 ∈ Z, ω
(36)
where, cos(ψ1 (ω)) =
|P2 |2 − |P0 |2 − |P1 |2 , ψ1 ∈ [0, π], 2|P0 P¯1 |
with the following condition on ω: 2 2 2 |P2 | − |P0 | − |P1 | ≤ 2|P0 P¯1 |. Let us consider
(37)
n o Ω1 = ω ∈ R+ : |P2 |2 − |P0 |2 − |P1 |2 ≤ 2|P0 P¯1 | .
By squaring both sides of the two conditions (32) and (37), one can show that (32) and (37) are equivalent. Thus, Ω1 = Ω2 ≡ Ω. We call Ω the crossing set.
Lemma 4. The set Ω consists of a finite number of intervals of finite length. Proof: See Appendix B.9. Remark 5. Since F (0) > 0, then F (ω) has roots 0 < a1 < b1 ≤ a2 < b2 < ... ≤ aN < bN < ∞, (N ≤ 6) and Ω=
N [
Ωk , where Ωk = [ak , bk ].
k=1
We allow bk−1 = ak with concern that ak may be a root with even multiplicity. 17
It is obvious that the stability crossing curves obtained in this way should be the same as I in (35). Since τ1 , τ2 obtained in both ways depend only on ω ∈ Ω, they should be the same respectively and can be defined as in (33) and (36) as well to get the crossing curves, i.e., n o I= τ1,n1 (ω), τ2,n2 (ω) ∈ R2+ , ω ∈ Ω, n1 , n2 = −1, 0, 1, 2... . Since Ω =
N [
Ωk , we have
k=1
I=
[
k=1,2,...,N
Ink1 ,n2 ,
(38)
n1 =−1,0,1,2... n2 =−1,0,1,2...
where Ink1 ,n2
n o 2 = τ1 (ω), τ2 (ω) ∈ R+ , ω ∈ Ωk , n1 , n2 = −1, 0, 1, 2... .
(39)
One can show that each Ink1 ,n2 is continuous in R2+ . In addition Ink1 ,n2 may connect with each other. Similar analysis can be found in Gu et al. 2005 [14]. B- Crossing direction Next, we will discuss the direction in which the solutions of the characteristic equation given by D(ξ, τ1 , τ2 ) = 0 cross the imaginary axis as (τ1 , τ2 ) deviate from a curve in Ink1 ,n2 . We will call the direction of the curve that corresponds to increasing ω, the positive direction. Observe that, the curve passes through the points corresponding to the end points of Ωk , the positive direction is reversed. For the purpose of discussing the direction of crossing, we need to consider τ1 and τ2 as functions of ξ = σ + iω, i.e., functions of two real variables σ and ω, and partial derivative notation needs to be adopted instead. Since the tangent of Ink1 ,n2 along the ∂τ ∂τ 2 1 , , the normal to Ink1 ,n2 pointing to the right region of the positive direction is ∂ω ∂ω ∂τ ∂τ1 2 positive direction is ,− . Also, as a pair of complex conjugate solutions of ∂ω ∂ω D(ξ, τ1 , τ2 ) = 0 cross the imaginary axis to the C+ , (τ1 , τ2 ) moves along the direction ∂τ ∂τ 1 2 , . We can therefore conclude that the inner product of these two vectors are ∂σ ∂σ positive, i.e., ∂τ ∂τ ∂τ ∂τ1 1 2 2 Γ(ω) = , . ,− = det ∆(ω) > 0, ∂σ ∂σ ∂ω ∂ω where ∂τ1 ∂τ1 ∂σ ∂ω ∆(ω) = ∂τ . ∂τ2 2 ∂σ ∂ω σ=0,ω∈Ω
The region at the right of Ink1 ,n2 has two more characteristic roots with positive real parts. On the other hand, if the previous inequality is reversed, then the region at the left of Ink1 ,n2 has two more characteristic roots with positive real parts. By the implicit function theory, we have ˆ1 R ˆ 2 −1 R0 −I0 R ∆(ω) = , I0 R0 Iˆ1 Iˆ2 18
where, R0 = I0 = ˆ1 = R ˆ2 = R Iˆ1 = Iˆ2 = Since det
∂Re D(ξ, τ1 , τ2 ) ∂σ ξ=iω ∂Im D(ξ, τ1 , τ2 ) ∂σ ξ=iω ∂Re D(ξ, τ1 , τ2 ) ∂τ1 ξ=iω ∂Re D(ξ, τ1 , τ2 ) ∂τ2 ξ=iω ∂Im D(ξ, τ1 , τ2 ) ∂τ1 ξ=iω ∂Im D(ξ, τ1 , τ2 ) ∂τ2 ξ=iω R0 −I0 I0 R0
) 2 X , = Re P00 (iω) + Pk0 (iω) − τk Pk (iω)e−iωτk k=1 ) ( 2 X , = Im P00 (iω) + Pk0 (iω) − τk Pk (iω)e−iωτk (
k=1
−iωτ1
= ω 2 S1 cos ωτ1 − ωS2 sin ωτ1 , = Re − iωP1 (iω)e −iωτ2 = ω 2 T1 cos ωτ2 , = Re − iωP2 (iω)e −iωτ1 = −ω 2 S1 sin ωτ1 − ωS2 cos ωτ1 , = Im − iωP1 (iω)e = Im − iωP2 (iω)e−iωτ2 = −ω 2 T1 sin ωτ2 .
ˆ 1 Iˆ2 − R ˆ 2 Iˆ1 ). = R02 + I02 > 0, we have sign(Γ(ω)) = sign(R
ˆ 1 Iˆ2 − R ˆ 2 Iˆ1 = −ω 2 T1 sin ωτ2 ω 2 S1 cos ωτ1 − ωS2 sin ωτ1 R 2 2 + ω T1 cos ωτ2 ω S1 sin ωτ1 + ωS2 cos ωτ1 , h i = ω 2 T1 ωS2 cos(ω(τ1 − τ2 )) + ω 2 S1 sin(ω(τ1 − τ2 )) , 1 −τ2 ) = ω 2 Im P¯1 P2 eiω(τ , o n = −ω 2 |P0 P¯1 |Im eiφ1 +iωτ1 , = −ω 2 |P0 P¯1 | sin φ1 + ωτ1 .
Hence, we have the following result.
Theorem 9. For any k = 1, 2, ..., N , (i) If sin φ1 + ωτ1 < 0, then the region at the right of Ink1 ,n2 has two more characteristic roots with positive real parts. (ii) If sin φ1 + ωτ1 > 0, then the region at the left of Ink1 ,n2 has two more characteristic roots with positive real parts. Remark 6. If we know the number of characteristic roots with positive real parts when τ1 = τ2 = 0, we can use the criterion in Theorem 9 to find the number of characteristic roots with positive real parts for any (τ1 , τ2 ) in R2+ . Hence, the stability of the characteristic equation is completely known. 6. Numerical Simulations Using a nonstandard numerical scheme as in Anguelov et al. 2012 [2], we provide the following numerical simulations where parameter values have been chosen in such a way that they obey the conditions for stability or bifurcation. 19
6.1. Numerical Simulations when τ1 = τ2 = 0 In this section, we present some numerical simulations to illustrate the analytical results observed in the previous sections. We first consider the following set of parameter values: r = 0.01, β = 0.001, γ = 2, αb = 0.1, αg = 0.15, δ = 0.73, a = 0.6, c1 = 0.32, c2 = 0.3, c3 = 0.15, eb = 0.45, db = 0.015, hb = 0.15, eg = 0.65, dg = 0.75, hg = 0.8. Figure 1 depicts monostability as well as two types of bistabilities situations when considering the browser’s half saturation parameter b as a bifurcation parameter. Namely, monostability of the equilibrium E1 = (K, 0, 0), bistability between limit cycle and equilibrium E1 for b ∈ [0.2, 0.9] and the bistability between equilibrium E1 = (K, 0, 0) and equilibrium E˜ = (R, Hb , 0) for b ∈ [0.9, 1.5]. Specifically, in Figure 2, we illustrate these two bistabilities situation for two particular values of the browser’s half saturation parameter b. L0 >1, S0<1
14
12
Browsers (Hb)
10
Stable equilibrium (R, 0, 0)
(1)
Stable equilibrium (R, Hb , 0)
(1)
Unstable equilibrium (R, Hb , 0)
8
Stable limit cycle in the 2D R−Hb plan
6
(2)
Unstable equilibrium (R, Hb , 0)
4
2
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
Bifurcation parameter: b
Figure 1: A bifurcation diagram of system (1) without delays and when considering the browser’s half saturation parameter b as a bifurcation parameter.
20
Bistability of E1=(R, 0, 0) and E2=(R, Hb, 0) with b=1.4
15
• Browsers
0.1
Browsers
10
0.08 0.06 0.04 0.02
Zoom
0 0
•0.02
0.04 0.06 0.08
0.1
Resource
5
0 0
5
10
15
Resource
Bistability of E =(R, 0, 0) and the limit cycle with b=0.4 1
15
Browser
0.8
Browsers
10
0.6 0.4 0.2 0
5
0 0
0
•
0.05
0.1
Zoom
1
0.15
0.2
0.25
0.3
0.35
Resource
2
3
4
5
6
Resource
Figure 2: Phase portrait of resource and browsers without delays. Illustration of bistability between equi˜ = (R, Hb , 0) (top panel). Phase portrait of resource and browsers librium E1 = (K, 0, 0) and equilibrium E without delays. Illustration of bistability between equilibrium E1 = (K, 0, 0) and the limit cycle (bottom panel).
Let us now assume that r = 1.4, β = 0.1, γ = 2, αb = 0.1, αg = 0.015, δ = 0.33, a = 0.1, b = 0.2, c1 = 0.9, c2 = 0.03, c3 = 0.15, eb = 0.15, db = 0.015, hb = 0.15, eg = 0.3, dg = 0.3, hg = 0.4. In Figure 3 we illustrate a (forward) bifurcation phenomenon involving equilibria E¯ and E1 when we consider the grazer’s catch rate c1 as a bifurcation parameter. Indeed, for L0 > 1 i.e. c1 < 0.55, equilibrium E1 is stable while for L0 < 1 i.e. c1 > 0.55, equilibrium E1 becomes unstable and equilibrium E¯ is stable (see Figure 3-(i)).
21
(i)
(ii)
Forward bifurcation in system (1)
Phase portrait in the Resource−Grazers 2D plan
0.8
4.5 4
0.7
Stable equilibrium (R, 0, Hg)
Stable equilibrium (R, 0, 0)
Unstable equilibrium (R, 0, 0)
3
0.5
Grazers
g
Grazers (H )
0.6
3.5
0.4 0.3
L0>1
2.5 2 1.5
L0<1
0.2
1
0.1 0 0
•
0.5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0 0
1
1
2
3
4
5
6
7
8
9
Resource
Bifurcation parameter: c1
(iii)
(iv)
Phase portrait in the Resource−Browsers 2D plan
Phase portrait in the Browsers−grazers 2D plan
2
6
1.8 5
1.6 1.4 Grazers
Browsers
4 1.2 1 0.8
2
0.6 0.4
1
•
0.2 0 0
3
1
•
2
3
4
5
6
7
8
0 0
9
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Browsers
Resource
Figure 3: (i) Illustration, without delays, of the stability of equilibrium E1 = (K, 0, 0) for the bifurcation ¯ = (R, 0, Hg ) for c1 > 0.55. (ii) Phase portrait of resource parameter c1 < 0.55 and stability of equilibrium E ¯ = (R, 0, Hg ). (iii) Phase portrait of and grazers without delays. Illustration of stability of equilibrium E ¯ = (R, 0, Hg ). (iv) Phase resource and browsers without delays. Illustration of stability of equilibrium E ¯ = (R, 0, Hg ). portrait of browsers and grazers without delays. Illustration of stability of equilibrium E
6.2. Numerical Simulations with positive delays In what follows, the numerical experiments are performed on system (1) to illustrate our theoretical findings when delays are considered. From the above discussion when τ1 > 0 and τ2 = 0, we may determine the direction of the Hopf bifurcation and the stability of the bifurcating periodic solution. We consider the following set of parameter values: r = 0.01, β = 0.001, γ = 2, αb = 0.1, αg = 0.015, δ = 0.728, a = 0.6, b = 1.2, c1 = 0.322, c2 = 0.3, c3 = 0.15, eb = 0.45, eg = 0.65, db = 0.015, dg = 0.75, hb = 0.15, hg = 0.8. For these parameter values, conditions textcolorredin Theorem 4, page 11 are satisfied and system (1) becomes conditionally stable around the boundary equilibrium point E˜ for τ1 ∈ [0, τ10 ) (see, Figure 4) and unstable for τ1 > τ10 (e.g. see Figure 5). Indeed, when τ1 passes through d the critical value τ1 = τ10 = 0.3587 and (Reξ(τ1 ))|τ1 =τ10 = 0.1570 > 0, the equilibrium dτ1 E˜ losses its stability and the system (1) experiences a Hopf bifurcation. From Sections 5 and 5.1.2, we can determine the nature of the stability and direction of the periodic solution bifurcating from the boundary equilibrium E˜ at the critical point τ10 . We can compute C2 (0) = 4.2289 + 1.5425, µ2 = −26.9288 < 0, β2 = 8.4578 > 0 and T2 = −31.5892 < 0. Hence, the bifurcating periodic solution exists and the corresponding periodic solution is subcritical as illustrate in Figure 6. The period of the periodic solution is 31.5892 and its period decreases with τ1 since T2 < 0. 22
Figure 4: Dynamics of system (1) when τ1 = 0.1793 < τ10 .
Figure 5: Dynamics of system (1) when τ1 = 1.4348 > τ10 .
˜ = (1.3732, 13.1438, 0) of system Figure 6: Hopf-bifurcation behavior of system (1) around the equilibrium E ˜ (1) at τ1 = τ10 = 0.3587 and existence of subcritical bifurcating periodic solution around E.
From the above discussion when τ1 = 0 and τ2 > 0, we may determine the direction of Hopf-bifurcation and the direction of bifurcating periodic solution. We assume that: r = 4.1, β = 0.2, γ = 4.01, αb = 0.2, αg = 0.015, δ = 0.328, a = 0.1, b = 0.2, c1 = 0.9, c2 = 0.03, c3 = 0.15, eb = 0.15, eg = 0.3, db = 0.15, dg = 0.3, hb = 0.5, hg = 0.4. For these parameter values, conditions in Theorem 8 are satisfied and system (1) becomes 23
conditionally stable around the boundary equilibrium point E¯ for τ2 ∈ [0, τ20 ) (see, Figure 7) and unstable for τ2 > τ20 (see, Figure 8). d When τ2 passes through the critical value τ2 = τ20 = 0.8183 and (Reξ(τ2 ))|τ2 =τ20 = dτ2 0.0225 > 0, the equilibrium E¯ losses its stability and the system (1) experiences Hopf bifurcation.
Figure 7: Dynamics of system (1) when τ2 = 0.40915 < τ20 .
Figure 8: Dynamics of system (1) when τ2 = 1.6366 > τ20 .
24
¯ = (18.1235, 0, 199.6832) Figure 9: Hopf-bifurcation behavior of system (1) around the interior equilibrium E 0 of system (1) at τ2 = τ2 = 0.8183 and existence of subcritical bifurcating periodic solution around the ¯ interior equilibrium E.
7. Conclusion-discussion In this paper, we construct and analyze a delayed competition model between herbivores on a same resource in a given location in a Sahelian region. We consider three state variables namely, the resource biomass (R), the grazer biomass (Hg ) and the browser biomass (Hb ). As modelling option we assume that herbivory can stimulate resource biomass production. Consequently the herbivory optimization hypothesis (HOH) is considered for the resource growing term. Recall that the HOH was also used in others herbivory systems (e.g. Lebon et al. 2014 [19]). We also assume a Holling type III functional response for browsers while a Beddington-DeAngelis like-functional response is considered for grazers. Moreover, we also acknowledged two discrete delays as to account for times necessary to convert resource biomass into herbivore (grazers and browsers) biomasses. Finally, since for the Sahelian region drought can lead a low biomass production, we also consider the case where the resource intrinsic growth rate can be negative. To authors knowledge, it is the first time that a mathematical model designed to model resource-grazers-browsers dynamics in Sahelian region is presented and studied. We first analyzed the delays-free case and secondly study the case where delays are considered. In both cases, the model exhibited several equilibria including extinction equilibrium (both resource and herbivores biomasses are null), herbivores-free equilibrium, grazersfree equilibria, browsers-free equilibria and coexistence equilibria. Specifically, we found that, thanks to the HOH and the Holling type III functional response for browsers, up to two grazers-free equilibria may exist. Similarly, thanks to the HOH and the BeddingtonDeAngelis functional response for grazers, up to three browsers-free equilibria may also exist. Concerning the coexistence equilibria, we found that up to two coexistence equilibria may exist (see Proposition 1, page 6). We further carry out stability analysis of various equilibria in the delays-free case. Based on our mathematical analysis, we found that both resource and herbivores may go to extinction as a consequence of drought for example (i.e. the resource intrinsic growth rate r is such that r < 0). This phenomenon can motivate pastor to decide to migrate with the cattle 25
(transhumance). We also found several thresholds that shape stability of equilibria. Notably, we found that several multistability situations may exist involving aforementioned equilibria (see Theorem 1, page 8). Recall that a multistability situation occurs when several equilibria are simultaneously stable. In addition, we found that when varying some parameters, several bifurcations occurs including Hopf bifurcation (e.g. Theorem 4, page 9) and forward bifurcation. Indeed, with assumed parameter values, we obtained two bifurcation diagrams with the browser’s half saturation parameter b and the grazer’s catch rate c1 as bifurcation parameters respectively (see Figures 1 and 3). Notably, Figure 1 depicts monostability of herbivores-free equilibrium as well as two types of bistability situations when considering the browser’s half saturation parameter b as a bifurcation parameter. Specifically, we found bistability between the herbivores-free equilibrium and a limit cycle when b ∈ [0.2; 0.9]. When b ∈ [0.9; 1.5], we obtain a bistability situation involving the herbivores-free and the grazer-free equilibria. In Sahelian context, Hopf bifurcation denotes that herbivores and resource biomass exhibit periodic behavior. Similarly when the grazer-free (resp. browser-free) equilibrium is stable, its may be the consequence of migration of grazer (reps. browser) by pastor to another location. When discrete delays are considered, we obtain sufficient conditions on parameters for which the delay-induced system is asymptotically stable around an equilibrium (that was also stable in the delays-free case) for all values of the delay parameter. When the conditions are not satisfied, then there exists a critical value of the delay parameter below which the system is stable and above which the system is unstable. By applying the normal form theory and the center manifold theorem, the explicit formula which determine the stability and direction of the bifurcating periodic solutions are computed. Our analytical results show that when τ1 and τ2 pass through the critical values τ10 and τ20 , the equilibria E˜ and E¯ loss their stability and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate ¯ We further illustrate these results through numerical simulations (see Figures from E˜ and E. 4-9). A limitation of this study relies on the non-availability of parameter values that correspond to Sahelian regions. In this study we assume sets of parameter values in such a way that they fulfilled stability or bifurcation conditions. However, it is desirable that some field experiments may be carried out in order to assess at least the magnitude of (some) parameters. For example, assess the magnitude or the range of the browser’s half saturation parameter b and the grazer’s catch rate c1 who are two bifurcation parameters of the model. This study is a first step of a more general study that aims to provide through mathematical modelling some insights for the browsers, grazers and their resource interactions in Sahelian region. This study in a given location in Sahelian can be improved in several ways. First, one can take into account both temporal and spatial behaviors of browsers, grazers and their resource interactions. It can be done through only local spatial operator such as Laplace operator or with both local and non-local spatial behaviors with some kernels. Another improvement may consist on taking into account several spatial locations and look for the best strategy in order to reach both sustainable resource management in the different locations and a good herbivores nutrition. This last point will be the aim of a forthcoming paper. Acknowledgements We would like to thank the referees for their helpful comments and valuable suggestions. 26
References [1] U. African. African Union Policy Framework for Pastoralism in Africa: Securing, Protecting and Improving the Lives, Livelihoods and Rights of Pastoralist Communities. African Union, 2010. [2] R. Anguelov, Y. Dumont, and J.M. Lubuma. On nonstandard finite difference schemes in biosciences. AIP Conf. Proc., 1487:212–223, 2012. [3] M.L. Bates, R.A. Cropp, D.W. Hawker, and J. Norbury. Which functional responses preclude extinction in ecological population-dynamic models? Ecological Complexity, 26:57–67, 2016. [4] J. Bell. A grazing ecosystem in the serengeti. Sci Am., 225:86–89, 1971. [5] R.S. Cantrell and C. Cosner. On the dynamics of predator-prey models with the beddington-deangelis functional response. Journal of Mathematical Analysis and Applications, 257:206–222, 2001. [6] M. Carriere. Impact des syst`emes d’´elevage pastoraux sur l’environnement en afrique et en asie tropicale et sub-tropicale aride et sub-aride. livestock and the environment finding a balance. Scientific Environmental Monitoring Group, 1996. [7] G. Caughley. Plant-herbivore systems. Theoretical Ecology, Blackwell Scientific Publications, 320:530, 1976. [8] J.V. De Koppel and H.H.T. Prins. The importance of herbivore interactions for the dynamics of african savanna woodlands: an hypothesis. Journal of Tropical Ecology, 14:565–576, 1998. [9] D.L. DeAngelis, R.A. Goldstein, and R.V. O’Neill. A model for tropic interaction. Ecology, 56:881–892, 1975. [10] B. Dedieu, V. Ancey, P. Bommel, S. Cournut, L. Dobremez, B. Faye, A. Gibon, E. Josien, A. Ickowicz, B. Lemery, C. Macombe, S. Madelrieux, C.R. Poccard, H. Rapey, G. Serviere, and J.F.d Tourran. Transformations de l’´elevage et dynamiques des espaces. Projet TRANS, Rapport scientifique final, 2009. [11] Z. Feng, Z. Qiu, R. Liu, and D.L. DeAngelis. Dynamics of a plant herbivore-predator system with plant-toxicity. Mathematical Biosciences, 229:190–204, 2011. [12] H.I. Freedman and R.V.S. Hari. The trade-off between mutual interference and time lags in predator-prey systems. Bull. Math. Biol., 45:991–1004, 1983. [13] S. Gratton. Analyse matricielle et optimisation, edition ressources pedagogiques. Ouv. INPT, 0727, 2014. [14] K. Gu, S. Niculescu, and J. Chen. On stability crossing curves for general systems with two delays. J. Math. Anal. Appl., 311:231–252, 2005. [15] J.K. Hale. Asymptotic behavior of dissipative systems. American Mathematical Society, Providence, 1988. 27
[16] B.D. Hassard, N.D. Kazarinoff, and Y.H. Wan. Theory and applications of hopf bifurcation. Cambridge University, 1981. [17] D.W. Hilbert, D.M. Swift, J.K. Detling, and M.I. Dyer. Relative growth rates and the grazing optimization hypothesis. Oecologia, 51:14–18, 1981. [18] Y. Kuang. Delay differential equation with application in population dynamics. Academic Press New York, 1993. [19] A. Lebon, L. Mailleret, Y. Dumont, and F. Grognard. Direct and apparent compensation in plant-herbivore interactions. Ecological Modelling, 290:192–203, 2014. [20] Y. Li, Z. Feng, R. Swihart, J. Bryant, and N. Huntly. Modeling the impact of plant toxicity on plant-herbivore dynamics. Journal of Dynamics and Differential Equations, 18(4), 2006. [21] X. Lin and H. Wang. Stability analysis of delay differential equation with two discrete delays. Canadian Applied Mathematics Quarterly, 20, 2012. [22] A. Martin and S. Ruan. Predator-prey models with delay and prey harvesting. J. Math. Biol., 43:247–267, 2001. [23] S.J. McNaughton. Grazing lawns: animals in herds, plants form and coevolution. The American Naturalist, 124(6):863–886, 1984. [24] N. Nori, M. Taylor, and A. Sensi. Droits pastoraux, modes de vie et adaptation au changement climatique. Irish Aid, 2008. [25] I. Noy-Meir. Stability of grazing systems: An application of the predator-prey graphs. Journal of Ecology, 63:459–481, 1975. [26] Rapport. L’´elevage au sahel et en afrique de l’ouest. 26`eme runion annuelle du r´eseau de pr´evention des crises alimentaires (rpca). Accra (Ghana), 2010. [27] Rapport. Atlas des ´evolutions des syst`emes pastoraux au sahel, atlas of trends in pastoral systems in sahel. FAO et CIRAD, 2012. [28] C.V. Rebecca and S. Ruan. A delay-differential equation model of hiv infection of cd4 t-cells. Mathematical Biosciences (Elsevier), 165:27–39, 2000. [29] Research reports and studies. Humanitarian policy group, pastoralism policies and practice in the horn and east africa a review of current trends. 2009. [30] L.A. Shipley. Grazers and browsers: How digestived morphology affects diet selection. Journal of Tropical Ecology, 14:565–576, 1999. [31] A.R.E. Sinclair. Does interspecific competition or predation shape the african ungulate community? Anim Ecol, 54:899–918, 1985. [32] A.R.E Sinclair and M. Northon-Griffiths. Does competition or facilitation regulate migrant ungulate populations in the serengeti? a test of hypotheses. Oecologia, 53:364– 369, 1982. 28
[33] A.B. Smith. Origins and spread of pastoralism in africa. Nomadic Peoples, 32:91–105, 1993. [34] H. Smith. An introduction to delay differential equations with applications to the life sciences. Springer, New York, 2010. [35] J.J. Tewa, A. Bah, and S.C.O. Noutchie. Dynamical models of interactions between herds forage and water resources in sahelian region. Hindawi Publishing Corporation Abstract and Applied Analysis, 2014(13), 2014. [36] J.J. Tewa, S. Bowong, and N.S.C. Oukouomi. Mathematical analysis of a two–patch model of tuberculosis disease with staged progression. Applied Mathematical Modelling, 36(12):5792–5807, 2012. [37] J.J. Tewa, V.Y. Djeumen, and S. Bowong. Predator-prey model with holling response function of type ii and sis infectious disease. Applied Mathematical Modelling, 37(7):4825–4841, 2013. [38] S.C. Williamson, J.K. Detling, J.L. Dodd, and M.I. Dyer. Experimental evaluation of the grazing optimization hypothesis. Journal of Range Management ptimum gmrin9 Intensity, 42, 1989. [39] V. Yatat, Y. Dumont, J.J. Tewa, P. Couteron, and S. Bowong. Mathematical analysis of a size structured tree-grass competition model for savanna ecosystems. Biomath, (18), 2014.
29
Appendix A. Some results for the browser-free case Consider the following System when Hb = 0 i.e. dR α H α H c1 Hg g g g g = r 1+ 1− R − δR2 − R, dt β γ 1 + c2 Hg and let
dHg c1 Hg = eg R − (dg + hg )Hg , dt 1 + c2 Hg
αg Hg αg Hg c1 Hg r(1 + )(1 − )− β γ 1 + c2 Hg eg c1 RHg , f2 (R, Hg ) = µg (1 + c2 Hg ) F (R, Hg ) = (f1 , f2 )(R, Hg ). 1 f1 (R, Hg ) = δ
An equilibrium of system (A.1) satisfies dR(t) = 0 dt
(A.1)
,
f1 (R, Hg ) = R,
⇔ dH (t) f2 (R, Hg ) = Hg . g = 0 dt r , 0 as an equilibrium of (A.1). The computations give z0 = δ r Lemma A1. If L0 > 1, then ρ(DF (z0 )) < 1 where z0 = , 0 and ρ(A) denotes the δ spectral radius of matrix A. Proof: Let us consider the jacobian matrix at z0 . ∂f1 ∂f1 rαg (γ − β) 0 ∂R ∂Hg δβγ DF (z0 ) = = e g c1 r ∂f2 ∂f2 0 δµg ∂R ∂Hg (z=z0 ) Then ρ(DF (z0 )) =
.
eg c1 r 1 = . Since L0 > 1, one has δµg L0 ρ(DF (z0 )) < 1.
Lemma A2. If ρ(DF (z0 )) < 1, then there exists a matrix norm such that (see Gratton 2014 [13]): kDF (z0 )k < 1. 30
Proof: The following inequality holds ∀ε > 0, kDF (z0 )k ≤ ρ(DF (z0 )) + ε. For ε =
1 − ρ(DF (z0 )) , we have 2
kDF (z0 )k < 1.
Lemma A3. If kDF (z0 )k < 1, then F is strictly contracting. Proof: By continuity of DF at z0 we have: ∀ε > 0, ∃ρ > 0/ z ∈ B(z0 , ρ) ⇒ kDF (z)k < k(DF (z0 ))k + ε; ∀ε > 0, ∃ρ > 0/
sup z∈B(z0 ,ρ)
∀ε > 0, ∃ρ > 0/ kDF k = For ε =
kDF (z)k < kDF (z0 )k + ε;
sup z∈B(z0 ,ρ)
k(DF (z))k < kDF (z0 )k + ε.
1 − kDF (z0 )k , we have 2 ∃ρ > 0/
sup z∈B(z0 ,ρ)
kDF (z)k <
1 + kDF (z0 )k < 1, 2
i.e. kDF (z0 )k < 1. Lemma A4. If L0 > 1, then F is strictly contracting, and System (A.1) has a unique r equilibrium z0 = ,0 . δ
Proof: Since L0 > 1, Lemma r A1, Lemma A2 and Lemma A3 imply that F is strictly contracting. Moreover, z0 = , 0 is a fixed point of F . Then, from the fixed point theorem, δ System (A.1) has unique equilibrium z0 . Consequently, the following results hold. Lemma A5. If L0 > 1, then System (A.1) does not admit a coexistence equilibria. Lemma A6. If L0 < 1, then System (A.1) has one or three coexistence equilibria.
Appendix B. Proofs Appendix B.1. Proof of Lemma 1 (1) The resource’s equation is given by: dR(t) αb Hb + αg Hg αb Hb + αg Hg aR2 c1 Hg =r 1+ 1− R−δR2 − Hb − R, 2 dt β γ b+R 1 + c2 Hg + c3 Hb 31
Since αb Hb + αg Hg < γ, then dR γ ≤r 1+ R. dt β Therefore, r 1+
R(t) ≤ R(0)e
γ! β .
Since r ≤ 0, then R(t) → 0. If R(t) → 0, then Hb → 0 and Hg → 0. (2) (i) Let X(t) = (R, Hb , Hg ) (t) be a solution of System (1). Let suppose that the solution X starts in the positive orthant R3+ . To prove that system (1) is positively invariant, we only need to show that at any border of R3+ , either the resulting vector field stays on the border or is pointing inside R3+ . Case 1: {R = 0}. At this border and following the first equation of System (1), dR one has = 0. Consequently, the solution X cannot crosses R3+ throught the dt border {R = 0}. dHb Case 2: {Hb = 0}. Following the second equation of system (1), one has = 0. dt 3 Therefore, the solution X cannot crosses R+ throught the border {Hb = 0}. dHg Case 3: {Hg = 0}. Following the third equation of system (1), one has = 0. dt 3 Therefore, the solution X cannot crosses R+ through the border {Hg = 0}. Consequently the part 1 of Lemma 1 holds. (ii) Let us now prove that the set Ωε , is a positively invariant region and absorbing set for the system (1). The resource’s equation is given by: dR(t) αb Hb + αg Hg αb Hb + αg Hg aR2 = r 1+ 1− R − δR2 − Hb dt β γ b + R2 c1 Hg − R, 1 + c2 Hg + c3 Hb which implies dR(t) αb Hb + αg Hg αb Hb + αg Hg ≤r 1+ 1− R − δR2 . dt β γ Since αb Hb + αg Hg < γ, then dR(t) αb Hb + αg Hg ≤r 1+ R − δR2 . dt β which implies
dR(t) γ ≤r 1+ R − δR2 dt β
Then, lim sup R(t) ≤ t−→+∞
32
r γ 1+ , δ β
and hence for ε > 0,
r γ R(t) ≤ 1+ + ε. δ β
Since 0 < αb Hb + αg Hg < γ, consequently the part 2 of Lemma 1 holds. This completes the proof of Lemma 1. Appendix B.2. Proof of Proposition 1 Equilibria in System (1) are obtained by solving the following system: αb Hb + αg Hg αb Hb + αg Hg c1 Hg R aR2 2 r 1 + 1 − R − δR − − H = 0, 2 b β γ 1 + c H + c H b + R 2 g 3 b aR2 Hb − µb Hb = 0, e b b + R2 c1 RHg eg − µg Hg = 0. 1 + c2 Hg + c3 Hb (B.1) r 1. The trivial equilibria of system (1) are given by E0 = (0, 0, 0) and E1 = , 0, 0 . δ 2. When R 6= 0, Hb 6= 0, and Hg = 0 αb Hb αb Hb aR Hb = 0, (i) r 1+ 1− − δR − β γ b + R2 (B.2) 2 aR eb − µb = 0. (ii) b + R2 r bµb eb a ˜= From (ii ) one has R with S01 > 1 where S01 = . eb a − µb µb ˜ in (i ) gives: Substituting R ! ˜ 1 βγ βγ δ R ˜ b2 − ˜b + H (γ − β) − σ˜ H − 1 = 0. (B.3) αb rαb R αb2 r ˜ β δR 1− and ηa = , then equation (B.3) leads to: γ r ˜ 2 − 1 (γ − β)(1 − 1 )H ˜ b + βγ (ηa − 1) = 0. H b αb ηs αb2
rαb Let us set ηs = βσR˜
The discriminant of equation (B.4) is given by: 2 1 1 4βγ 2 ∆ = 2 (γ − β) 1 − − 2 (ηa − 1) . αb ηs αb Let us set N0 =
4βγηa . 1 2 2 (γ − β) (1 − ) + 4βγ ηs 33
(B.4)
• If N0 > 1 i.e. ∆ < 0, then (B.4) does not admit a real solution. ˜ 0 defined by: • If N0 = 1 i.e. ∆ = 0, then (B.4) admits a unique solution H b 1 (γ − β) 1 − ηs ˜0 = H . b 2αb ˜ H ˜ 0 , 0) is ecologically meaningful whenever N0 = 1, S01 > 1 Equilibrium E2 = (R, b and ηs > 1. ˜ 1 and H ˜ 2 given by: • If N0 < 1 i.e. ∆ > 0, then (B.4) admits two solutions H b b ˜1 = H b
˜2 = H b
r
1 (γ − β)(1 − ) − ηs 1 (γ − β)(1 − ) + ηs
r
(γ − β)2 (1 −
1 2 ) − 4βγ(ηa − 1) ηs
2αb
(γ − β)2 (1 −
1 2 ) − 4βγ(ηa − 1) ηs
2αb
,
.
˜ H ˜ 1 , 0) is ecologically meaningful whenever N0 < 1, ηs > 1 Equilibrium E3 = (R, b ˜ H ˜ 2 , 0) is ecologically meaningful and ηa > 1. Similarly, equilibrium E4 = (R, b 1 whenever (N0 < 1, S0 > 1 and ηa < 1) or (N0 < 1, S01 > 1, ηs > 1 and ηa > 1). Table 1: Existence of solutions of equation (B.4) N0 >1 =1 <1 <1 <1 <1
ηs >1 >1 >1 <1 <1
ηa
E2 √
>1 <1 <1 >1
E3
E4
√
√ √ √ √
Number of equilibrium 0 1 2 1 1 1
3. When R 6= 0, Hb = 0 and Hg 6= 0 αg Hg αg Hg c1 eg Hg r 1+ 1− − δR − = 0; (a) β γ 1 + c2 Hg
(B.5)
c1 eg R − µg = 0. (b) 1 + c2 Hg ¯ = µ g 1 + c2 H ¯g . Equation (b) implies R eg c1 ¯ Substituting R in (a) gives:
¯ 2 + βγc2 + αg (γ − β) (1 − L1 ) H ¯ g + βγ (L0 − 1) = 0, (B.6) ¯ 3 + γ − β (1 − L2 ) H H g g αg L2 αg c2 L1 αg2 c2 where L0 =
δµg βγc2 + αg (γ − β) ; L1 = 2βγδµg c2 βγc1 and L2 = reg c1 + reg c1
r
αg c2
γ−β
+
βγδµg c2 rαg eg c1
.
We summarize the problem of existence of solutions of (B.6) in the following table. 34
Table 2: Existence of solutions of equation (B.6) L2 <1 <1 >1 <1 >1 >1 >1 <1
L1 <1 >1 <1 >1 <1 >1 >1 <1
L0 <1 >1 >1 <1 <1 >1 <1 >1
Number of sign changes Number of positive solutions 1 1 2 2 or 0 2 2 or 0 1 1 3 3 or 1 2 2 or 0 1 1 0 0
Using Lemmas A5 and A6, Table 2 becomes Table 3: Existence of solutions of equation (B.6) L2 <1 <1 >1 >1
L1 <1 >1 <1 >1
L0 <1 <1 <1 <1
Number of sign changes Number of positive solutions 1 1 1 1 3 3 or 1 1 1
4. When R 6= 0, Hb 6= 0 and Hg 6= 0 α aR α c1 Hg b Hb + αg Hg b Hb + αg Hg r 1+ − 1− − δR − Hb = 0; (i) β γ 1 + c2 Hg + c3 Hb b + R2 aR2 e − µb = 0; (ii) b b + R2 c1 eg R − µg = 0. (iii) 1 + c2 Hg + c3 Hb (B.7) r ∗ e c R − µ (1 + c2 Hg∗ ) bµ g 1 g b Equations (ii) and (iii) give respectively R∗ = and Hb∗ = . eb a − µb c3 µ g e g c1 R ∗ . Hb∗ > 0 ⇔ eg c1 R∗ − µg (1 + c2 Hg∗ > 0 i.e. N ∗ > 1 with N ∗ = µg (1 + c2 Hg∗ ) Substituting R∗ and Hb∗ in equation (i) gives: D2 Hg∗2 + D1 Hg∗ + D0 = 0, rαg2 (1 − m1 )2 ; βγ rαg (1 − m1 )(γ − β) 2rαb αg (1 − m1 )(1 − m0 ) σR∗ c2 µg D1 = − + − ; βγ βγc3 c3 eg R ∗ rαb (m0 − 1)(γ − β) rαb2 (m0 − 1)2 σR∗ eg c1 R∗ σR∗ D0 = r + − − δR∗ − + , βγc3 βγc3 µ g c3 c3 α b c2 eg c1 R∗ m0 = and m1 = . µg αg c3 where D2 = −
35
(B.8)
We summarize the problem of existence of solutions of equation (B.8) in the following Table. Table 3: Existence of solutions of equation (B.8) D2 <0 <0 <0 <0
D1 <0 <0 >0 >0
D0 <0 >0 >0 <0
Number of sign changes Number of positive solutions 0 0 1 1 1 1 2 2 or 0
Appendix B.3. Proof of Theorem 1 1. It is easy to verify that the eigenvalues of the Jacobian matrix at E0 = (0, 0, 0) are r, −µb and −µg . Therefore, it is an unstable hyperbolic saddle critical point. • If r < 0, then E0 is locally asymptotically stable.
• If r > 0, then E0 is a saddle-point and therefore unstable. r eb ar2 2. The eigenvalues of the Jacobian matrix at E1 = ( , 0, 0) are −r, 2 − µb = δ bδ + r2 e g c1 r 1 µb S0 − 1 and − µg = µg − 1 . Hence, equilibrium E1 is locally asymptotδ L0 ically stable if S0 < 1 and L0 > 1. The instability of E1 follows if one of the previous conditions is not satisfied. 3. The jacobian matrix at the point E˜ is given by: m ˜ 11 m ˜ 12 m ˜ 13 ˜ = m ˜ 21 0 0 , J(E) 0 0 m ˜ 33 where,
˜H ˜ b (R ˜ 2 − b) ˜2 aR rαb ˜ m ˜b R ˜ − aR ; − δ R; ˜ 12 = γ − β − 2αb H ˜ 2 )2 ˜2 βγ (b + R b+R ˜˜ ˜H ˜b ˜ 2abR eg c1 R m ˜ 21 = ;m ˜ 33 = − µg = µg S0RHb − 1 . ˜ 2 )2 ˜b (b + R 1 + c3 H The characteristic polynomial is given by: P (λ) = (m ˜ 33 − λ) λ2 − m ˜ 11 λ − m ˜ 21 m ˜ 12 . m ˜ 11 =
The characteristic polynomial admits negative roots if
m ˜ 33 < 0, m ˜ 11 < 0 and m ˜ 21 m ˜ 12 < 0. ˜ ˜
m ˜ 33 < 0 ⇔ S0RHb < 1, ˜H ˜ b (R ˜ 2 − b) ˜2H ˜b aR aR ˜ < 0 ⇔ S1R˜ H˜b < 1 with S1R˜ H˜b = m ˜ 11 = − δR . ˜ 2 )2 ˜ b + δ(b + R ˜ 2 )2 (b + R abH m ˜ 21 m ˜ 12 < 0 ⇔ m ˜ 12 < 0 because m ˜ 21 > 0, ˙ ∂R m ˜ 12 = = ψ 0 (Hb ) − G0 (Hb ) with ∂Hb αb Hb αb Hb ˜ + σ ˜ Hb . 1− and G(Hb ) = δ R ψ(Hb ) = r 1 + R β γ 36
Moreover,
dR = 0 is equivalent to dt ψ(Hb ) − G(Hb ) = 0.
therefore
(B.9)
dR = 0 is equivalent to dt ψ(Hb ) − G(Hb ) = 0.
(B.10)
Furthermore, by using relation (B.10) we deduce graphically the stability of equilib˜ see Figure B.10. rium E,
Figure B.10: We have three different configurations. Case (i): There is no equilibria. Case (ii): The unique equilibrium is locally asymptotically stable. Case (iii): There exist two equilibria and one is locally asymptotically stable.
4. The jacobian matrix at the point E¯ is given by: m ¯ 11 m ¯ 12 m ¯ 13 ¯ = 0 m ¯ 22 0 , J(E) m ¯ 31 0 m ¯ 33 where,
¯ ¯ rαb ¯g R ¯ − aR + c1 c3 Hg R; ¯ ¯ m γ − β − 2αg H m ¯ 11 = −δ R; ¯ 12 = ¯ 2 (1 + c2 H ¯ g )2 βγ b+R m ¯ 13 =
¯ ¯ rαg ¯g R ¯ − c1 R + c1 c2 Hg R; ¯ γ − β − 2αg H ¯ g (1 + c2 H ¯ g )2 βγ 1+H
¯2 ¯g ¯g aeb R eg c1 H µ g c2 H ¯ R − µ = µ S − 1 ; m ¯ = ; m ¯ = − b b 31 33 0 ¯2 ¯g ¯g . b+R 1 + c2 H (1 + c2 H The characteristic polynomial, given by: m ¯ 22 =
P (λ) = (m ¯ 22 − λ) [λ2 − (m ¯ 11 + m ¯ 33 )λ + m ¯ 11 m ¯ 33 − m ¯ 13 m ¯ 31 ], admits negative roots if m ¯ 22 < 0, m ¯ 11 + m ¯ 33 < 0, 37
and m ¯ 11 m ¯ 33 − m ¯ 13 m ¯ 31 > 0.
m ¯ 11 + m ¯ 33 ¯
¯g µ g c2 H ¯ = − δR + ¯ g < 0. (1 + c2 H
Observe that m ˜ 22 < 0 ⇔ S0R < 1, m ¯ 11 m ¯ 33 − m ¯ 13 m ¯ 31 > 0 if m ¯ 13 < 0, m ¯ 13 =
∂ R˙ = ψ 0 (Hg ) − T 0 (Hg ), with ∂Hg αg Hg c1 Hg αg Hg 1− and T (Hg ) = δR + . ψ(Hg ) = r 1 + β γ 1 + c2 Hg
Moreover,
dR = 0 is equivalent to dt ψ 0 (Hg ) − T 0 (Hg ) = 0.
(B.11)
Therefore, m ¯ 13 < 0 is equivalent to ψ 0 (Hg )−T 0 (Hg ) < 0. Furthermore, by using relation ¯ see Figure B.11. (B.11), we deduce graphically the stability of equilibrium E,
¯ Case (a): The unique equilibrium is Figure B.11: We have different configurations for the stability of E. locally asymptotically stable. Case (b): there exist three equilibria and two are locally asymptotically stable.
38
Appendix B.4. Proof of Theorem 2 The characteristic polynomial of the jacobian matrix at the point E˜ is given by: P (λ) = (m ˜ 33 − λ) λ2 − m ˜ 11 λ − m ˜ 21 m ˜ 12 ,
˜ H˜b (R ˜ 2 − b) aR ˜ Observe that where m ˜ 11 = − δ R. ˜ 2 )2 (b + R P (λ) = 0 ⇔ λ = m ˜ 33 or λ2 − m ˜ 11 λ − m ˜ 21 m ˜ 12 = 0. (B.12) From expression of m ˜ 11 , E˜ is the only candidate for a Hopf bifurcation. It follows from expression of m ˜ 11 that if a Hopf bifurcation occurs, there exists δh that m ˜ 11 = 0 or equivalently ˜ H˜b (R ˜ 2 − b) aR ˜ = 0 i.e. − δR 2 2 ˜ (b + R ) ˜ 2 − b) aH˜b (R . δh = ˜ 2 )2 (b + R The discriminant of equation (B.12) is given by ∆ = m ˜ 211 + 4m ˜ 21 m ˜ 12 . 1. If δ 6= δh and ∆ < 0, then the eigenvalues of P (λ) are given by: λ1,2 (δ) = α1,2 (δ) + iβ1,2 (δ). 2. If δ = δh and ∆ < 0, then λ1,2 (δh ) = iβ1,2 (δh ). 3. Moreover,
∂α1,2 ˜ 6= 0. (δh ) = −R ∂δ
∂α1,2 (δh ) < 0, and the conclusion follows : ∂δ • δ = δh is a bifurcation value.
Thus,
• there exists δh1 > δh such that for all δ ∈ [δh , δh1 ], the equilibrium point E˜ is a stable focus. ˜ there exists δh2 < δh such that for every δ ∈ • for every neighborhood U of E, [δh2 , δh [, the equilibrium E˜ is an unstable focus surrounded by a√stable limit cycle contained in U , whose amplitude increases and is of the order δh − δ.
Therefore, in Theorem 2, we have a supercritical Hopf bifurcation with the appearance of a stable limit cycle.
Appendix B.5. Proof of Theorem 3 The jacobian matrix at the point E ∗ is given by: ∗ m11 m∗12 m∗13 0 , J(E ∗ ) = m∗21 0 ∗ ∗ m31 m32 m∗33 39
where,
αb Hb∗ + αg Hg∗ c1 Hg∗ 2abR∗ ∗ 1− −2δR − − Hb∗ ; ∗ ∗ ∗2 2 γ 1 + c2 Hg + c3 Hb (b + R )
m∗11
αb Hb∗ + αg Hg∗ =r 1+ β
m∗12
∗ c1 c3 Hg∗ rαb aR∗2 ∗ ∗ R∗ ; = γ − β − 2(αb Hb + αg Hg ) R − + ∗ 2 ∗2 ∗ βγ b+R (1 + c2 Hg + c3 Hb )
m∗13 =
m∗21 =
rαb 2eb abR∗ ∗ c1 R∗ (1 + c3 Hb∗ ) ∗ = γ − β − 2(αb Hb∗ + αg Hg∗ ) R∗ − ; m H 31 βγ (1 + c2 Hg∗ + c3 Hb∗ )2 (b + R∗2 )2 b eg c1 c3 Hg∗ R∗ µg c2 Hg∗ eg c1 Hg∗ ∗ ∗ = − = − ; m ; m . 32 33 1 + c2 Hg∗ + c3 Hb∗ (1 + c2 Hg∗ + c3 Hb∗ )2 1 + c2 Hg∗ + c3 Hb∗
The characteristic polynomial is given by: P (λ) = −λ3 + (m∗11 + m∗33 )λ2 + (m∗13 m∗31 + m∗12 m∗21 − m∗11 m∗33 )λ + m∗12 m∗21 m∗33 − m∗13 m∗21 m∗32 . The Routh-Hurwitz conditions for stability of this equilibrium are C1 = m∗11 + m∗33 < 0, C2 = (m∗11 + m∗33 )(m∗11 m∗33 − m∗13 m∗31 − m∗12 m∗21 ) + (m∗12 m∗21 m∗33 − m∗13 m∗21 m∗32 ) < 0, C3 = m∗12 m∗21 m∗33 − m∗13 m∗21 m∗32 > 0. When these conditions are satisfied, a coexistence equilibrium is locally asymptotically stable. Appendix B.6. Proof of Theorem 4 Since ω0 is a solution of equation (11), the characteristic equation (8) has the pair of purely imaginary roots ±iω0 . From equations (9) and (10), we have cos ω0 τ1 =
B3 ω0 4 + (B1 B4 − B2 B3 )ω0 2 . B32 ω0 2 + B42
0 Hence, τ1,n for n = 0, 1, ... as a function of ω0 is given by 1 B3 ω0 4 + (B1 B4 − B2 B3 )ω0 2 2πn 0 τ1,n = arccos + . 2 2 ω0 B3 ω0 2 + B4 ω0
(B.13)
For τ1 = 0, Theorem 1 ensures that E˜ is locally asymptotically stable. Hence, by Butler’s 0 lemma Freedman et al. 1983 [12] E˜ remains stable up to the minimum value of τ1,n , obtained 0 here thus for n = 0 i.e. for τ1 < τ1,0 , so that 0 0 = τ1,0 . τ10 = min τ1,n n≥0
The theorem will be completely proven if we can show that d(Reξ(τ1 )) > 0. sign dτ1 ξ=iω0 40
Differentiating equation (18) with respect to τ1 yields dξ 2 3ξ + 2B1 ξ + B2 + B3 e−ξτ1 − τ1 (B3 ξ + B4 )e−ξτ1 = ξ(B3 ξ + B4 )e−ξτ1 , dτ1
which implies that −1 dξ 3ξ 2 + 2B1 ξ + B2 + B3 e−ξτ1 − τ1 (B3 ξ + B4 )e−ξτ1 . = dτ1 ξ(B3 ξ + B4 )e−ξτ1 Since (B3 ξ + B4 )e−ξτ1 = −(ξ 3 + B1 ξ 2 + B2 ξ), we obtain −1 dξ 3ξ 2 + 2B1 ξ + B2 B3 τ1 = + − , 3 2 dτ1 −ξ(ξ + B1 ξ + B2 ξ) ξ(B3 ξ + B4 ξ =
ξ 2 + B1 ξ + B2 + 2ξ 2 + B1 ξ B3 τ1 + − , 2 2 −ξ (ξ + B1 ξ + B2 ) ξ(B3 ξ + B4 ξ
=
2ξ 2 + B1 ξ −B4 τ1 + − . −ξ 2 (ξ 2 + B1 ξ + B2 ) ξ 2 (B3 ξ + B4 ξ
Thus, ( −1 ) d(Reξ(τ1 )) dξ sign = sign Re , dτ1 dτ1 ξ=iω0 ξ=iω0 2ξ 2 + B1 ξ −B4 = sign Re + Re 2 , −ξ 2 (ξ 2 + B1 ξ + B2 ) ξ=iω0 ξ (B3 ξ + B4 ) ξ=iω0 1 2B32 ω0 8 + (3B42 + B32 θ20 )ω0 6 + 2B42 θ20 ω0 4 + B42 B22 ω0 2 . = sign ω0 2 [B12 ω0 4 + (ω0 3 − B2 ω0 )2 ] [B42 + B32 ω0 2 ] Since θ20 > 0, we have
d (Reξ(τ1 )) > 0. dτ1 ξ=iω0 Hence, the transversality condition is satisfied and a Hopf bifurcation occurs at τ1 = τ10 . This achieves the proof. Appendix B.7. Proof of Theorem 5 ˜ x2 (t) = Hb (t) − H ˜ b and x3 (t) = Hg (t); then system (1) is equivalent Let x1 (t) = R(t) − R, to the following three dimensional system: dx1 (t) ˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b x1 (t) + ψH (H ˜ b , 0)R ˜ − P (R) ˜ x2 (t) = ψ( H b dt ˜ ˜ ˜ ˜ + ψ ( H , 0) R − g ( H , 0) R x (t) + f x (t), x (t), x (t) ; Hg b Hg b 3 1 1 2 3 dx2 (t) ˜ ˜ = e P ( R) H x (t − τ ) + f x (t), x (t), x (t), x (t − τ ) ; b b 1 1 2 1 2 3 1 1 dt dx3 (t) = e g (H ˜ b , 0) − µg x3 (t) + f3 x1 (t), x2 (t), x3 (t) , g Hg dt 41
(B.14)
where 2 ˜ ˜ ˜ f1 x1 (t), x2 (t), x3 (t) = ψ x2 (t) + Hb , x3 (t) x1 (t) + R − δ x1 (t) + R ˜ − g x2 (t) + H ˜ b , x3 (t) x1 (t) + R ˜ − P x1 (t) + R h ˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b x1 (t) + ψH (H ˜ b , 0)R ˜ − P (R) ˜ x2 (t) − ψ(H b +
i ˜ b , 0)R ˜ − gHg (H ˜ b , 0)R ˜ x3 (t) ; ψHg (H
f2 x1 (t), x2 (t), x3 (t), x1 (t − τ2 )
˜ ˜ b ) − µb (x2 (t) + H ˜ b) = eb P x1 (t) + R (x2 (t) + H
˜ H ˜ b x1 (t − τ1 ); − eb P (R)
and ˜ b , x3 (t) x1 (t) + R ˜ − µg x3 (t) f3 x1 (t), x2 (t), x3 (t) = eg g x2 (t) + H ˜ − eg gHg (Hb , 0) − µg x3 (t).
Let τ1 = τ10 + µ, then µ = 0 is the Hopf bifurcation value of system (1) at the positive ˜ H ˜ b , 0). Since system (1) is equivalent to system (B.14), in the following equilibrium E˜ = (R, discussion we shall consider mainly the system (B.14). In system (B.14), let ui (t) = xi (τ1 t). Then, system (B.14) can be rewritten as follows: h du1 (t) 0 ˜ b , 0)R ˜ − P (R) ˜ u2 (t) ˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b u1 (t) + ψH (H + µ) ψ( H = (τ 1 b dt i ˜ ˜ ˜ ˜ ( H , 0) R − g ( H , 0) R u (t) + f u (t), u (t), u (t) ; + ψ 1 2 3 b 3 1 Hg b Hg h i du2 (t) 0 ˜ ˜ + µ) e P ( R) H u (t − τ ) + f u (t), u (t), u (t), u (t − τ ) ; = (τ 1 2 3 1 1 b b 1 1 2 1 dt h i du3 (t) = (τ 0 + µ) e g (H ˜ , 0) − µ u (t) + f u (t), u (t), u (t) . g H b g 3 3 1 2 3 g 1 dt
(B.15) System (B.15) can be written as functional differential equation (FDE) in C([−1, 0], R3 ) of the form: u(t) ˙ = Lµ (ut ) + f (µ, ut ), (B.16)
where, u(t) = (u1 (t), u2 (t), u3 (t))T , Lµ : C −→ R is the linear operator and f : (R×C) −→ R is the nonlinear operator, given respectively by
s11 s12 s13
Lµ (φ) = (τ10 + µ) 0 0
0 0
0
0 0
φ1 (0) φ2 (0) + t21 0 0 0 φ3 (0) s33 0 0 0 42
φ1 (−1) φ2 (−1) (B.17) φ3 (−1)
˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b , s12 = ψH (H ˜ b , 0)R ˜ − P (R), ˜ where s11 = ψ(H b ˜ ˜ ˜ ˜ ˜ ˜ H ˜b s13 = ψHg (Hb , 0)R − gHg (Hb , 0)R, s33 = eg gHg (Hb , 0) − µg and t21 = eb P 0 (R) and f1 (φ1 (0), φ2 (0), φ3 (0)) , f (φ (0), φ (0), φ (0), φ (−1)) f (φ, ν) = (τ10 + µ) 2 1 2 3 1 f3 (φ1 (0), φ2 (0), φ3 (0))
(B.18)
where φ = (φ1 , φ2 , φ3 )T ∈ C. By the Riesz representation theorem, there exists a (3 × 3) matrix, η(θ, µ), (−1 ≤ θ ≤ 0) whose elements are bounded variation functions such that Z 0 dη(θ, µ)φ(θ) f or φ ∈ C([−1, 0], R3 ). (B.19) Lµ (φ) = −1
In fact, we can choose
η(θ, µ) =
(τ10
s11 s12 s13
+ µ) 0 0
0 0
0
0 0
t21 0 0 0 δ(θ) + s33 0 0 0
where δ is the Dirac function defined by 0, δ(θ) = 1,
θ 6= 0;
δ(θ + 1) ,
(B.20)
(B.21)
θ = 0.
For φ ∈ C([−1, 0], R3 ), we define
dφ(θ) , θ ∈ [−1, 0); dθ A(µ)φ(θ) = Z 0 dη(µ, s)φ(s), θ = 0.
(B.22)
−1
and
R(µ)φ(θ) = Then, the system (B.15) is equivalent to
0,
θ ∈ [−1, 0);
f (µ, θ),
θ = 0.
u(t) ˙ = A(µ)ut + R(µ)ut ,
43
(B.23)
(B.24)
where ut (θ) = u(t + θ) for θ ∈ [−1, 0]. For ψ ∈ C([−1, 0], R3 ), we define dψ(s) − , s ∈ (0, 1]; ds A∗ ψ(s) = Z 0 dη T (t, 0)ψ(−t), s = 0.
(B.25)
−1
and a bilinear inner product
¯ hψ(s), φ(θ)i = ψ(0)φ(0) −
Z
0
−1
Z
θ
ξ=0
¯ − θ)dη(θ)φ(ξ)dξ, ψ(ξ
(B.26)
where η(θ) = η(θ, 0). Then, A(0) and A∗ are adjoint operators. By Theorem 4, we know that ±iω0 τ10 are eigenvalues of A(0). Thus, they are also eigenvalues of A∗ . Let us consider q(θ) as the eigenvector of A(0) corresponding to +iω0 τ10 and q ∗ (θ) the eigenvector of A∗ corresponding to −iω0 τ10 . T T 0 0 Then, q(θ) = 1, β3 , γ3 eiθω0 τ1 and q ∗ (s) = D 1, β3∗ , γ3∗ eiω0 τ1 s . From the above discussion, it is easy to know that A(0)q(0) = iω0 τ10 q(0) and A∗ (0)q ∗ (0) = −iω0 τ10 q ∗ (0). That is
τ10
and
τ10
Then,
s11 s12 s13
0 0
s11 0
and
0
0 0
0 0 q(0) + τ1 t21 0 0 s33 0 0 0
0 0
0
s12 0 0 s13 0 s33
0 t21 0
∗ 0 q (0) + τ1 0 0
iω0 − s11
−t21 e−iω0 τ10 0
0 0
−s12
−s13
iω0
0
0
iω0 − s33
44
q(−1) = iω0 τ10 q(0),
∗ 0 ∗ 0 q (−1) = −iω0 τ1 q (0). 0
1
0
β3 = 0 γ3 0
(B.27)
(B.28)
(B.29)
0
−iω0 − s11 −t21 e−iω0 τ1
0
−s12
−iω0
0
−s13
0
−iω0 − s33
Thus, we can easily obtain
1
∗
D1 1, β3∗ , γ3∗
q (s) = where
T
0
∗ β3 = 0 . ∗ γ3 0
T 0 q(θ) = 1, β3 , 0 eiθω0 τ1 ,
and
(B.30)
(B.31)
0
eiω0 τ1 s ,
(B.32)
˜ b , 0)R ˜ ˜ b , 0)R ˜ + gHg (H ˜ b , 0)R ˜ + P (R) ˜ ˜ H ˜ b e−iω0 τ10 s −ψHg (H −ψHb (H eb P 0 (R) , β3∗ = , γ3∗ = β3 = − ¯ g )R ¯ − µg . iω0 iω0 iω0 + eg gHg (0, H In order to assure hq ∗ (s), q(θ)i = 1, we need to determine the value of D1 . From (B.26), we have Z 0Z θ ∗ ∗ q¯∗ (ξ − θ)dη(θ)φ(ξ)dξ, hq (s), q(θ)i = q¯ (0)q(0) − −1 0
Z
ξ=0 θ
Z
T ¯ 1, β¯∗ , γ¯ ∗ e−iω0 τ10 (ξ−θ) dη(θ) 1, β3 , 0 eiω0 τ10 ξ dξ, D 3 3 −1 ξ=0 0 0 0 ˜ H ˜ b 0 0 (−e−iω0 τ10 )q(0), = q¯∗ (0)q(0) − q¯∗ τ10 eb P 0 (R) 0 0 0 = q¯∗ (0)q(0) −
¯ 1 [1 + β3 β¯3∗ + τ10 β¯3∗ eb P 0 (R) ˜ H ˜ b e−iω0 τ10 ]. = D Therefore,
1 , ∗ 0 ¯∗ ¯ ˜ H ˜ b e−iω0 τ10 1 + β3 β3 + τ1 β3 eb P 0 (R) 1 . = ∗ 0 ∗ ¯ ˜ H ˜ b eiω0 τ10 1 + β3 β3 + τ1 β3 eb P 0 (R)
D¯1 =
(B.33)
D1
(B.34)
Using the same notations as in Hassard et al. 1981 [16], we first compute the coordinates to describe the center manifold C2 (0) at µ = 0. Let ut be the solution of Equation (B.24) when µ = 0. Define z0 (t) = hq ∗ , ut i, W (t, θ) = ut (θ) − 2Re(z0 (t)q(θ)) = ut (θ) − z0 (t)q(θ) + z¯0 (t)¯ q (θ) . (B.35) On the center manifold C2 (0), we have
W (t, θ) = W (z0 , z¯0 , θ), 45
(B.36)
where, W (z0 , z¯0 , θ) = W20 (θ)
z02 z¯2 z3 + W11 (θ)z0 z¯0 + W02 (θ) 0 + W30 (θ) 0 + ..., 2 2 6
(B.37)
where z0 and z¯0 are local coordinates for center manifold C2 (0) in the direction of q ∗ and q¯∗ . Note that W is real if ut is real. We only consider real solutions. For solution ut ∈ C2 (0) of (B.15), since µ = 0, we have def z˙0 (t) = iω0 τ10 z0 + q¯∗ (0)f 0, W (z0 , z¯0 , 0) + 2Re(z0 q(θ)) = iω0 τ10 z0 + q¯∗ (0)f0 (z0 , z¯0 ). (B.38)
We rewrite this equation as
z˙0 (t) = iω0 τ10 z0 + p(z0 , z¯0 ),
(B.39)
where, p(z0 , z¯0 ) = q¯∗ (0)f0 (z0 , z¯0 ) = p20 (θ)
z02 z¯2 z 2 z¯0 + p11 (θ)z0 z¯0 + p02 (θ) 0 + p21 (θ) 0 + ..., 2 2 2
(B.40)
Remark 7. In what follows, our objective is to compute coefficients p20 , p11 , p02 and p21 of p(z0 , z¯0 ). These coefficients will be use to find the direction of Hopf bifurcation. T 0 Let us recall, ut (θ) = u1t (θ), u2t (θ), u3t (θ) and q(θ) = 1, β3 , 0 eiω0 τ1 θ . So from (B.35) and (B.37), it follows that ut (θ) = W (t, θ) + 2Re(z0 (t)q(θ)), z¯2 z02 + W11 (θ)z0 z¯0 + W02 (θ) 0 2 2 T T 0 0 + 1, β3 , 0 eiω0 τ1 θ z0 + 1, β¯3 , 0 e−iω0 τ1 θ z¯0 + ... = W20 (θ)
Then,
(1)
z02 z¯2 (1) (1) + W11 (0)z0 z¯0 + W02 (0) 0 + ..., 2 2
u1t (0)
= z0 + z¯0 + W20 (0)
u2t (0)
z z¯ (2) (2) (2) = β3 z0 + β¯3 z¯0 + W20 (0) 0 + W11 (0)z0 z¯0 + W02 (0) 0 + ..., 2 2
u3t (0)
= W20 (0)
2
(3)
(B.41)
2
z02 z¯2 (3) (3) + W11 (0)z0 z¯0 + W02 (0) 0 + ..., 2 2
0
0
(1)
u1t (−1) = e−iω0 τ1 z0 + eiω0 τ1 z¯0 + W20 (−1)
z02 z¯2 (1) (1) + W11 (−1)z0 z¯0 + W02 (−1) 0 + ..., 2 2 2
2
z z¯ 0 0 (2) (2) (2) u2t (−1) = β3 e−iω0 τ1 z0 + β¯3 eiω0 τ1 z¯0 + W20 (−1) 0 + W11 (−1)z0 z¯0 + W02 (−1) 0 + ..., 2 2 (3)
u3t (−1) = W20 (−1)
z02 z¯2 (3) (3) + W11 (−1)z0 z¯0 + W02 (−1) 0 + .... 2 2 46
(B.42)
Also taking into account Eq.(B.18), it follows that
where,
p(z0 , z¯0 ) = q¯∗ f0 (z0 , z¯0 ), = q¯∗ f (0,ut ), ¯ 1 1, β¯3 , 0 f (0, ut ), = τ10 D m1 ¯ 1 1, β¯3 , 0 × m2 , = τ10 D m3
(B.43)
2 r h ˜ b αb u2t (0)u1t (0) + αg u3t (0)u1t (0) − αb u2t (0) + αg u3t (0) γ − β − 2αb H βγ i ˜ − δu2 (0) − P 0 (R)u ˜ 2t (0)u1t (0) − 1 P 00 (R)u ˜ 2t (0)u2 (0) − 1 P 00 (R) ˜ H ˜ b u2 (0) u1t (0) + R 1t 1t 1t 2 2 2 ˜ b , 0)u3t (0)u1t (0) − 1 ∂ g (H ˜ b , 0) u2 (0)u1t (0) + Ru ˜ 2 (0) − gHg (H 3t 3t 2 ∂Hg2 ∂ 2g ˜ ˜ + (Hb , 0) u3t (0)u2t (0)u1t (0) + Ru3t (0)u2t (0) ∂Hb Hg
m1 =
˜ 2t (0)u2 (−1) + 1 eb P 00 (R) ˜ H ˜ b u2 (−1) + ... ˜ 2t (0)u1t (−1) + 1 eb P 00 (R)u m2 = eb P 0 (R)u 1t 1t 2 2 2 ˜ b , 0)u2t (0)u1t (0) + 1 eg ∂ g (H ˜ b , 0) u23t (0)u1t (0) + Ru ˜ 23t (0) m3 = eg gHb (H 2 ∂Hg2 ∂ 2g ˜ ˜ (Hb , 0) u3t (0)u2t (0)u1t (0) + Ru3t (0)u2t (0) + .... + eg ∂Hb Hg
p(z0 , z¯0 ) is given by n h ¯1 r ˜ b αb u2t (0)u1t (0) + αg u3t (0)u1t (0) − αb u2t (0) p(z0 , z¯0 ) = τ10 D γ − β − 2αb H βγ 2 i ˜ − δu21t (0) − P 0 (R)u ˜ 2t (0)u1t (0) − 1 P 00 (R)u ˜ 2t (0)u21t (0) + αg u3t (0) u1t (0) + R 2 1 00 ˜ ˜ 2 1 ∂ 2g ˜ 2 2 ˜ ˜ − P (R)Hb u1t (0) − gHg (Hb , 0)u3t (0)u1t (0) − (Hb , 0) u3t (0)u1t (0) + Ru3t (0) 2 2 ∂Hg2 ∂ 2g ˜ b , 0) u3t (0)u2t (0)u1t (0) + Ru ˜ 3t (0)u2t (0) + (H ∂Hb Hg h io 1 1 0 ˜ 00 ˜ 2 00 ˜ ˜ 2 ¯ + β3 eb P (R)u2t (0)u1t (−1) + eb P (R)u2t (0)u1t (−1) + eb P (R)Hb u1t (−1) + ... 2 2
This expression can be extended to
P 00 (R) ˜ H ˜b z02 n 0 ¯ h rαb β3 2˜ 2 0 ˜ ˜ p(z0 , z¯0 ) = 2τ1 D1 γ − β − 2αb Hb − αb Rβ3 − δ + β3 P (R) − 2 βγ 2 47
io ˜ b P 00 (R)e ˜ −2iω0 τ10 ˜ 3 e−iω0 τ10 + 1 eb H + β¯3 eb P 0 (R)β 2
P 00 (R) ˜ H ˜b z¯02 n 0 ¯ h rαb β¯3 2 ˜ ¯2 0 ˜ ¯ ˜ + 2τ1 D1 γ − β − 2αb Hb − αb Rβ3 − δ + β3 P (R) − 2 βγ 2 io ˜ β¯3 eiω0 τ10 + 1 eb H ˜ b P 00 (R)e ˜ 2iω0 τ10 + β¯3 eb P 0 (R) 2
n h rα (β + β¯ ) b 3 3 0 ¯ 2˜ 0 ˜ ¯ ¯ ˜ + z0 z¯0 τ1 D1 γ − β − 2αb Hb − 2αb Rβ3 β3 − 2δ − β3 + β3 )P (R βγ io 0 0 ˜ b P 00 (R) ˜ + β¯3 eb P 0 (R)Re{( ˜ ˜ b P 00 (R) ˜ − H β¯3 e−iω0 τ1 + β3 eiω0 τ1 )} + eb H z02 z¯0 n 0 ¯ h r (2) (1) (2) ˜ b 2αb W11 τ1 D1 γ − β − 2αb H (0) + 2β3 αb W11 (0) + αb W20 (0) 2 βγ (1) (3) (3) + β3 αb W20 (0) + 2αg W11 (0) + αg W20 (0) − 2 αb2 β32 + 2β3 β¯3 αb2 +
+
(2) 2W11 (0)
+
(2) W20 (0)
(3) (3) ˜ + 2αb αg Rβ3 W11 (0) + αb αg W20 (0)
(1) (1) (2) (1) ˜ 2W11 − 2δ 2W11 (0) + W20 (0) − P 0 (R) (0) + 2β3 W11 (0)
(2) (1) (1) (1) ˜ β¯3 + 2β3 + 2H ˜ b W11 ˜ b W20 + W20 (0) + β¯3 W11 (0) − P 00 (R) (0) + H (0) ∂ 2g (3) (3) (3) (3) ˜ b , 0) 2W11 ˜ b , 0) 2β3 W11 − gHg (H (0) + W20 (0) − (H (0) + W20 (0) ∂Hb Hg 0 0 (2) (1) (2) (1) ˜ W11 + β¯3 eb P 0 (R) (0)e−iω0 τ1 + β3 W11 (−1) + W20 (0)eiω0 τ1 + β¯3 W11 (−1)
io 1 0 0 0 (1) (1) ˜ 2H ˜ b W11 ˜ b W20 eb P 00 (R) (−1)e−iω0 τ1 + H (−1)eiω0 τ1 + 2β¯3 e−2iω0 τ1 + 4β3 2 Comparing the coefficients with (B.40), gives h 00 ˜ ˜ ¯ 1 rαb β3 γ − β − 2αb H ˜ b − α2 Rβ ˜ 2 − δ + β3 P 0 (R) ˜ − P (R)Hb p20 = 2τ10 D 3 b βγ 2 i ˜ b P 00 (R)e ˜ −2iω0 τ10 , ˜ 3 e−iω0 τ10 + 1 eb H + β¯3 eb P 0 (R)β 2 +
p02
h rα β¯ P 00 (R) ˜ H ˜b b 3 2 ˜ ¯2 0 ˜ ¯ ˜ γ − β − 2αb Hb − αb Rβ3 − δ + β3 P (R) − = βγ 2 i ˜ β¯3 eiω0 τ10 + 1 eb H ˜ b P 00 (R)e ˜ 2iω0 τ10 , + β¯3 eb P 0 (R) 2
¯1 p11 = τ10 D
¯1 2τ10 D
h rα (β + β¯ ) b 3 3 ˜ b − 2α2 Rβ ˜ 3 β¯3 − 2δ − β3 + β¯3 )P 0 (R ˜ −H ˜ b P 00 (R) ˜ γ − β − 2αb H b βγ 48
p21
i 0 0 ˜ b P 00 (R) ˜ ˜ , + β¯3 eb P 0 (R)Re{ β¯3 e−iω0 τ1 + β3 eiω0 τ1 } + eb H h r (2) (1) (2) (1) ¯1 ˜ b 2αb W11 = τ10 D γ − β − 2αb H (0) + 2β3 αb W11 (0) + αb W20 (0) + β3 αb W20 (0) βγ (2) (2) (3) (3) + 2αg W11 (0) + αg W20 (0) − 2 αb2 β32 + 2β3 β¯3 αb2 + 2W11 (0) + W20 (0) (3) (3) (1) (1) ˜ 3 W11 + 2αb αg Rβ (0) + αb αg W20 (0) − 2δ 2W11 (0) + W20 (0) (2) (1) (2) (1) ˜ 2W11 − P 0 (R) (0) + 2β3 W11 (0) + W20 (0) + β¯3 W11 (0)
(1) (1) (3) (3) ˜ β¯3 + 2β3 + 2H ˜ b W11 ˜ b W20 ˜ b , 0) 2W11 − P 00 (R) (0) + H (0) − gHg (H (0) + W20 (0) −
∂ 2g (3) (3) ˜ b , 0) 2β3 W11 (H (0) + W20 (0) ∂Hb Hg
0 0 (1) (1) (2) (2) ˜ W11 + β¯3 eb P 0 (R) (0)e−iω0 τ1 + β3 W11 (−1) + W20 (0)eiω0 τ1 + β¯3 W11 (−1) +
1 ˜ e P 00 (R) 2 b
i 0 0 0 (1) (1) ˜ b W20 ˜ b W11 . (−1)eiω0 τ1 + 2β¯3 e−2iω0 τ1 + 4β3 2H (−1)e−iω0 τ1 + H
(B.44)
Since there are W20 and W11 in p21 , we still need to compute them. From (B.24) and (B.35), we have: ˙ W
=
=
def
=
u˙ t − z˙0 q − z¯˙0 q¯,
AW − 2Re(q¯∗ (0)f0 q(θ)),
θ ∈ [−1, 0),
AW − 2Re(q¯∗ (0)f0 q(θ)) + f0 ,
(B.45)
θ = 0,
AW + H(z0 , z¯0 , θ)
where
z02 z¯2 + H11 (θ)z0 z¯0 + H02 (θ) 0 + .... (B.46) 2 2 Substituting the corresponding series into Eq. (B.45), and comparing the coefficients, give H(z0 , z¯0 , θ) = H20 (θ)
(A − 2iω0 τ10 )W20 (θ) = −H20 (θ),
(B.47)
= −H11 (θ).
AW11 (θ) From Eq. (B.45) we know that for θ ∈ [−1, 0),
H(z0 , z¯0 , θ) = −q ∗ f0 q(θ) − q ∗ f¯0 q¯(θ) = −p(z0 , z¯0 )q(θ) − p¯(z0 , z¯0 )¯ q (θ).
(B.48)
Comparing the coefficients with Eq. (B.46), gives H20 (θ) = −p20 q(θ) − p¯02 q¯(θ), 49
(B.49)
H11 (θ) = −p11 q(θ) − p¯11 q¯(θ).
(B.50)
From Eqs. (B.46) and (B.49) and the definition of A, it follows that ˙ 20 (θ) = 2iω0 τ 0 W20 (θ) + p20 q(θ) + p¯02 q¯(θ). W 1
(B.51)
T 0 Notice that q(θ) = 1, β3 , 0 eiω0 τ1 θ . Hence,
i¯ p02 ip20 0 0 iω0 τ10 θ q(0)e + q¯(0)e−iω0 τ1 θ + E3 e2iω0 τ1 θ , (B.52) 0 0 ω 0 τ1 3ω0 τ1 (1) (2) (3) where, E3 = E3 , E3 , E3 ∈ R3 is a constant vector. Similarly, from (B.46) and (B.50)), we obtain W20 (θ) =
W11 (θ) = −
i¯ p11 ip11 0 0 q(0)eiω0 τ1 θ + q¯(0)e−iω0 τ1 θ + E4 , 0 0 ω 0 τ1 ω 0 τ1
(B.53)
(1) (2) (3) where, E4 = E4 , E4 , E4 ∈ R3 is also a constant vector. In what follows, we will seek appropriate E3 and E4 . From the definition of A and (B.46), we obtain Z 0
−1
and
dη(θ)W20 (θ) = 2iω0 τ10 W20 (0) − H20 (0), Z
(B.54)
0
−1
dη(θ)W11 (θ) = −H11 (0),
(B.55)
where, η(θ) = η(0, θ). By (B.45), we have H20 (0) = −p20 q(0) − p¯02 q¯(0) + 2τ10 ×
00 ˜ ˜ rαg β3 ˜ b ) − αb2 Rβ ˜ 32 − (δ + β3 P 0 (R)) ˜ − P (R)Hb (γ − β − 2α H b βγ 2 ˜ 3 e−iω0 τ10 + 1 eb H ˜ b P 00 (R)e ˜ −2iω0 τ10 eb P 0 (R)β 2 0
H11 (0) = −p11 q(0) − p¯11 q¯(0) + 2τ10 ×
,
(B.56)
rαb (β3 + β¯3 ) 2˜ 0 ˜ 00 ˜ ¯ ¯ ˜ ˜ (γ − β − 2αb Hb ) − 2αb Rβ3 β3 − 2δ − (β3 + β3 )P (R) − Hb P (R) βγ . 0 ˜ −iω0 τ10 iω0 τ10 00 ˜ ¯ ˜ eb P (R)Re{(β3 e + β3 e )} + eb Hb P (R) 0 (B.57)
50
Substituting Eqs. (B.52) and (B.56) into Eq. (B.54) with the fact that Z 0 0 0 iω0 τ1 I − eiω0 τ1 dη(θ) q(0) = 0, −1
− iω0 τ10 I −
Z
0
−iω0 τ10
e
−1
(B.58)
dη(θ) q¯(0) = 0.
gives,
h1 0 e2iω0 τ1 dη(θ) E3 = 2τ20 × h2 , 2iω0 τ10 I − −1 0 Z
0
00 ˜ ˜ rαg β3 ˜ b ) − αb2 Rβ ˜ 32 − (δ + β3 P 0 (R)) ˜ − P (R)Hb (γ − β − 2αb H βγ 2 1 00 ˜ −2iω0 τ10 0 ˜ −iω0 τ10 ˜ . + 2 eb Hb P (R)e and h2 = eb P (R)β3 e This leads to 2iω0 − h3 h4 h5 h 1 ˜ H ˜ b e−iω0 τ10 2iω0 eb P 0 (R) E3 = 2 h2 , 0 0 ˜ ˜ 0 0 2iω0 − eg gHg (Hb , 0)R + µg
where h1 =
˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b , h4 = −ψH (H ˜ b , 0)R ˜ + P (R) ˜ where, h3 = ψ(H b ˜ b , 0)R ˜ + gHg (H ˜ b , 0)R. ˜ and h5 = −ψHg (H Solving this system for E3 gives ˜ b , 0)R ˜ + P (R) ˜ −ψHg (H ˜ b , 0)R ˜ + gHg (H ˜ b , 0)R ˜ h1 −ψHb (H 2 (1) h2 2iω0 0 E3 = ∆3 0 ˜ b , 0)R ˜ + µg 0 2iω0 − eg gH (H g
(2)
E3
and
˜ b , 0) + 2δ R ˜ + P 0 (R) ˜ H ˜ b h1 −ψHg (H ˜ b , 0)R ˜ + gHg (H ˜ b , 0)R ˜ 2iω0 − ψ(H 2 ˜ H ˜ b e−iω0 τ10 = eb P 0 (R) h2 0 ∆3 ˜ b , 0)R ˜ + µg 0 0 2iω0 − eg gHg (H
(3)
E3
,
,
˜ b , 0) + 2δ R ˜ + P 0 (R) ˜ H ˜ b −ψH (H ˜ b , 0)R ˜ + P (R) ˜ h1 2iω0 − ψ(H b 2 0 0 −iω τ ˜ H ˜ be 0 1 = eb P (R) 2iω0 h2 = 0, ∆3 0 0 0 51
where, 2iω0 − h3 h4 h5 ˜ H ˜ b e−iω0 τ10 0 2iω0 0 ∆3 = eb P 0 (R) ˜ b , 0)R ˜ + µg 0 0 2iω0 − eg gHg (H
.
Similarly, substituting (B.53) and (B.57) into (B.55) gives h3 −h4 −h5 l 1 ˜ H ˜b 0 −eb P 0 (R) E4 = 2 l2 , 0 0 ˜ ˜ 0 0 eg gHg (Hb , 0)R − µg
rαb (β3 + β¯3 ) ˜ b ) − 2αb2 Rβ ˜ 3 β¯3 − 2δ − (β3 + β¯3 )P 0 (R) ˜ −H ˜ b P 00 (R) ˜ and (γ − β − 2αb H βγ 0 0 ˜ ˜ b P 00 (R). ˜ l2 = eb P 0 (R)Re{( β¯3 e−iω0 τ1 + β3 eiω0 τ1 )} + eb H Hence, ˜ b , 0)R ˜ − P (R) ˜ ψHg (H ˜ b , 0)R ˜ + gHg (H ˜ b , 0)R ˜ l1 ψHb (H 2 (1) , l 0 0 E4 = 2 ∆4 0 ˜ b , 0)R ˜ − µg 0 eg gH (H where l1 =
g
(2)
E4
and
˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b l1 ψHg (0, H ¯ g )R ¯ − gHg (0, H ¯ g )R ¯ ψ(H 2 ˜ H ˜b = −eb P 0 (R) l2 0 ∆4 ˜ b , 0)R ˜ − µg 0 0 eg gHg (H
(3)
E4
where
,
˜ b , 0) − 2δ R ˜ − P 0 (R) ˜ H ˜ b ψH (H ˜ b , 0)R ˜ − P (R) ˜ l1 ψ(H b 2 0 ˜ ˜ = 0, = −e P ( R) H 0 l b b 2 ∆4 0 0 0 h3 −h4 −h5 ˜ H ˜b 0 0 ∆4 = −eb P 0 (R) ˜ b , 0)R ˜ − µg 0 0 eg gHg (H 52
.
Thus, we can determine W20 and W11 from (B.52) and (B.53). Furthermore, g21 in (B.44) can be expressed using the parameters and delay. Thus, we can compute the following values: C2 (0) = µ2
|p02 |2 p21 i 2 p p − 2|p | − + , 20 11 11 2ω0 τ10 3 2
= −
Re{C2 (0)} , Re{ξ 0 (τ10 )}
β2
= 2Re{C2 (0)},
T2
= −
(B.59)
Im(C2 (0)) + µ2 Im(ξ 0 (τ10 )) . ω0 τ10
Appendix B.8. Proof of Lemma 3 It suffices to prove that there exists ω ¯ ∈ R+ such that |P1 |2 − |P0 |2 − |P2 |2 = 0.
(B.60)
Note that equation (B.60) is equivalent to ω ¯6 + ω ¯ 4 (R12 − 2R2 ) + ω ¯ 2 (T12 + R22 − S12 ) − S22 = 0.
(B.61)
By the intermediate value theorem the above equation has at least one positive root for each sign of (T12 + R22 − S12 ). This achieves the proof. Appendix B.9. Proof of Lemma 4 Let 2 2 2 2 F (ω) = |P1 | − |P0 | − |P2 | − 4|P0 |2 |P¯2 |2 .
There exists ω ˆ such that F (ˆ ω ) = 0. In fact, we have F (0) = S22 > 0. By lemma 3, there 2 exists ω ¯ such that |P1 | − |P0 |2 − |P2 |2 = 0. Hence, F (¯ ω ) = −4|P0 |2 |P¯2 |2 < 0 and ω ˆ ∈]0, ω ¯ [. Moreover, F (ω) is a polynomial and then, has a finite number of roots. This concludes the proof.
53