以遲滯型微分方程評估登革熱重症危險性

全文

(1)國立屏東大學應用數學系碩士班碩士論文 Department of Applied Mathematics National Pingtung University Master’s Thesis. 以遲滯型微分方程評估登革熱重症危險性 Evaluating Risks Of Severe Dengue Diseases With Time Delays. 指導教授：鄭昌源博士 Advisor: Dr. Chang-Yuan Cheng 研究生：秦淳 Student: Chuen Chin. 中華民國 107 年 8 月 August 2018.

(2)

(3) 以遲滯型微分方程評估登革熱重症危險性 Evaluating Risks of Severe Dengue Diseases with Time Delays 摘要在本篇研究,我們以遲滯型微分方程評估三種危險性不同族群(M-H-L)之動態,此三個遲滯(dealy)分別為 M 與 H 之接觸時間及族群產生一種和兩種抗體之時間。藉由分析線性化系統的特徵方程式,我們探討平衡點的局部穩定性,週期解及 Hopf bifurcation。我們也在本模式中探討疫苗之效力,失敗率,其乃臨床上一重要課題。吾人的結果顯示,遲滯可導致 Hopf bifurcation。在數值模擬方面,在三個遲滯均不為零及其中一個遲滯不為零的狀況下,可出現週期解。本模式可用於研究登革熱高危險性族群之病例量(disease burden),因其可因此病造成死亡, 此遲滯也可用於探討三劑疫苗注射之最佳時機。. 關鍵字: 登革熱, 遲滯型微分方程, Hopf 分歧, 疫苗效力. I.

(4) Abstract In this study we analyze the dynamical behavior of moderate-high-low(MHL) groups with three delays, which are MH contact time, production time of one and two types of antibodies to elicit. By analyzing the characteristic equations of the linearized systems, we establish the stability of equilibria and existence of Hopf bifurcation when three delays are used as the bifurcation parameters. Vaccine efficacy and failure rate are incorporated into our model and this has a clinical importance. Our results show that individual delay can lead to the Hopf bifurcation under some parameter value conditions. We numerically illustrate periodic solutions develop under the conditions of three nonzero delays and only one of which is nonzero. From the viewpoint of epidemiology, this study has two main results: (1) the used mathematical tools provide us a method to investigate the disease burden and the size of high risk group which can lead to disease-induced mortality; (2) estimating the time delay also provides us a method to determine the timing of individual vaccine dose.. Key words: dengue fever, delay differential equation, Hopf bifurcation, vaccine efficacy. II.

(5) 目錄摘要……………………………………………………………………………………………............ I. Abstract…………………………………………………………………………………………………. II 目錄………………………………………………………………………………………. III. 圖表索引…………………………………………………………………………………... V. Chapter 1. 1. Preface……………………………………………………………………………. 1.1. Introduction……………………………………………………………………. 1. 1.2. Research motivation……………………………………………………….. 2. Literature review…………………………………………………………... 4. Chapter 2 2.1. A theoretical model for the dengue epidemic using delayed differential equation………………………………………... 4. Dynamics in a system with two delays……………………………. 7. Research topics and methods……………………………………….. 9. 3.1. Modelling………………………………………………………………………. 9. 3.2. Research methods…………………………………………………………. 12. Analysis………………………………………………………………………... 13. 4.1. Equilibrium……………………………………………………………………. 13. 4.2. Linearization………………………………………………………………….. 15. 4.3. Local stability and Hopf bifurcation……………………………….. 16. 2.2 Chapter 3. Chapter 4. III.

(6) Chapter 5. Numerical simulation…………………………………………………….. 28. (I). 𝜏1 ≠ 0, 𝜏2 ≠ 0, 𝜏3 ≠ 0………………………………………………... 30. (II). 𝜏1 ≠ 0, 𝜏2 = 0, 𝜏3 = 0……………………………………………….. 33. (III). 𝜏1 = 0, 𝜏2 ≠ 0, 𝜏3 = 0…………………………………………………. 35. (IV). 𝜏1 = 0, 𝜏2 = 0, 𝜏3 ≠ 0…………………………………………………. 38. Discussion……………………………………………………………………... 41. Reference…………………………………………………………………………………………….. 43. Chapter 6. IV.

(7) 圖表索引 Table 1. Numbers of positive roots of 𝐻 ∗ …………………………………………… 15 Table 2. Reference values of parameters……………………………………………... 29. Figure 1. [𝜏1 , 𝜏2 , 𝜏3 ]=[8,14,14]; Dynamics of moderate risk group……. 30 Figure 2. [𝜏1 , 𝜏2 , 𝜏3 ]=[8,14,14]; Dynamics of high risk group…………….. 31. Figure 3. [𝜏1 , 𝜏2 , 𝜏3 ]=[8,14,14]; Dynamics of low risk group……………... 32. Figure 4. [𝜏1 , 𝜏2 , 𝜏3 ]=[8, 0, 0]; Dynamics of moderate risk group……... 33. Figure 5. [𝜏1 , 𝜏2 , 𝜏3 ]=[8, 0, 0]; Dynamics of high risk group………………. 34 Figure 6. [𝜏1 , 𝜏2 , 𝜏3 ]=[8, 0, 0]; Dynamics of low risk group………………... 34. Figure 7. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 14, 0]; Dynamics of moderate risk group……. 35 Figure 8. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 14, 0]; Dynamics of high risk group…………….. 36 Figure 9. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 14, 0]; Dynamics of low risk group………………. 37. Figure 10. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 0, 14]; Dynamics of moderate risk group….. 38. Figure 11. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 0, 14]; Dynamics of high risk group…………... 39. Figure 12. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 0, 14]; Dynamics of low risk group……………. 40. Figure A. Study modelling………………………………………………………………….. 46. V.

(8) Evaluating risks of severe dengue diseases with time delays. VI.

(9) CHAPTER 1 PREFACE 1.1 Introduction Dengue is the most important human viral disease transmitted by arthropod vectors. There are an estimated 50–100 million cases of dengue fever (DF), and 250 000 to 500 000 cases of dengue hemorrhagic fever (DHF) in the world annually. DF and DHF are caused by the dengue viruses, which includes 4 serotypes DEN 1, 2, 3, and 4. Infection with one serotype provides life-long immunity to that virus but not to the others. Dengue viruses are maintained in an urban transmission cycle in tropical and subtropical areas by the mosquito Aedes aegypti. In some areas other Aedes species, such as Ae albopictus and Ae polynesiensis are also vectors. DHF is defined as an acute febrile illness with minor or major bleeding, thrombocytopenia (≤ 105 /μL), and evidence of plasma leakage documented by hemoconcentration (hematocrit increased by at least one-fifth or decreased by the same amount after intravenous fluid therapy), pleural or other effusions, hypoalbuminemia or hypoproteinemia. Dengue shock syndrome (DSS) is defined as DHF with signs of circulatory failure, including narrow pulse pressure (≤20 mm Hg), hypotension, or frank shock. Once shock develops the fatality rate could be high (12% to 44%). DSS and DHF are considered severe dengue viral diseases. The most widely accepted hypothesis for the pathogenesis of DHF in secondary 1.

(10) infections is immune enhancement. Antibody-dependent enhancement (ADE) of dengue viruses, is a process in which the infecting virus is complexed with nonneutralizing antibodies, thus enhancing phagocytosis by mononuclear cells (the primary site of virus replication in man). These monocytes release vasoactive mediators, which result in fulminant immune responses, including the vascular permeability and hemorrhagic manifestations that characterize DHF and DSS [1]. In short, preformed antibody responded by first dengue disease or vaccination, leads to a more severe clinical condition in the consequent infections.. 1.2 Research Motivation A candidate tetravalent dengue vaccine has been assessed in three clinical trials involving more than 35,000 children between the ages of 2 and 16 years in Asian– Pacific and Latin American countries in 2014 and 2015. An unexplained higher incidence of hospitalization for dengue in year 3 among children younger than 9 years of age was noted. In contrast, the risk among children 2 to 16 years of age was lower in the vaccine group than in the control group in another trial [2]. A shorter time interval between first and second dengue infections is associated with protection from clinical illness was found in a school based cohort in Thailand [3]. These phenomena showed intriguing effective antibody concentration and the time lag between two episodes of 2.

(11) dengue infections. To elucidate some aspects of antibodies response is help to successful dengue vaccination. We conduct a model of dengue dynamics among vaccinated people to evaluate their risks.. 3.

(12) CHAPTER 2 LITERATURE REVIEW In addition to aforementioned studies delay differential equations (DDEs) are frequently used in the description of various natural phenomena. The spectrum includes the context of romantic relationships, description of immune reactions, biochemical reactions modelling, tumor growth and treatment, avascular stage of tumor growth and the immunotherapy of cancers.. 2.1 A Theoretical Model For The Dengue Epidemic Using Delayed Differential Equation Application of SIR epidemiological model showed new trends [4]. The SIR model with demography, including births and deaths, was introduced as 𝑑𝑆 𝑑𝑡 𝑑𝐼 𝑑𝑡 𝑑𝑅 𝑑𝑡. = μ–βSI-μS, = βSI-γI-μI, = γI-μR,. μ. natural birth rate, natural death rate;. β. transmission rate;. γ. recovery rate.. S, I, R represent susceptible, infected and recovered population 4.

(13) respectively. If we assume the entire population is susceptible, then the average number of new infections per infectious individual is determined by basic reproduction number 𝑅0 =. 𝛽 𝛾+𝜇. .. The inclusion of demographic dynamics may allow a disease to die out or persist in a population in the long term. For this reason it is important to analyze the dynamics when the system is at equilibrium.. Elhia introduced the time delay to the SIR system as a new model [5]: 𝑑𝑆(𝑡) 𝑑𝑡 𝑑𝐼(𝑡) 𝑑𝑡 𝑑𝑅(𝑡) 𝑑𝑡. 𝐼(𝑡). = Λ− 𝛽S(𝑡)𝑁(𝑡) − 𝑑S(𝑡) − 𝑢(𝑡−𝜏)𝑆(𝑡−𝜏), 𝐼(𝑡). = 𝛽S(𝑡) 𝑁(𝑡) − (𝛾+𝑑+𝜖)I(𝑡), = 𝛾𝐼(𝑡) − 𝑑𝑅(𝑡) + 𝑢(𝑡 − 𝜏)𝑆(𝑡 − 𝜏).. The different states are described by the following parameters: (i) Λ is the recruitment rate of susceptibles; (ii) 𝛽 is the effective contact rate; (iii) 𝑑 is the natural mortality rate; (iv) 𝛾 is the recovery rate; (v) 𝜖 is the disease induced death rate; (vi) 𝑢 is a control, represents the percentage of susceptible individuals 5.

(14) being vaccinated per time unit; (vii) S, I, R, N represents susceptible, infected, recovered and total population. Elhia assumes at time 𝑡 only a percentage of susceptible individuals that have been vaccinated 𝜏 time unit ago, that is, at time 𝑡 − 𝜏, are removed from the susceptible class and added to the recovered class. The introduction of delay 𝜏 provides more information about the disease dynamics. Perez [6] propose that the mosquitoes that bit sick persons at moment t- 𝜏 and become infectious, will have adult and infective descendent at the moment t. This fact is the main one in his concept for the expansion of the epidemic: vertical transmission. 𝜏 is the time that it takes for the eggs to hatch and for larvae to reach the adult state. Taking simplicity he considers the effect of mortality due to Dengue fever equal to zero and Dengue Hemorrhagic Fever is not present. The model is represented following: 𝑑𝑆(𝑡) 𝑑𝑡 𝑑𝐼(𝑡) 𝑑𝑡 𝑑𝑅(𝑡) 𝑑𝑡. = − 𝛽S(𝑡)I(t- 𝜏) – bS(t) + bN, = 𝛽S(𝑡) I(t- 𝜏) – (b+𝛾)I(𝑡), = 𝛾I(𝑡) − b𝑅(𝑡).. b. natural birth rate, natural death rate;. β. transmission rate;. γ. recovery rate; 6.

(15) N. total population.. The basic reproduction number is 𝑅0 =. 𝛽 𝑏+γ. in this model. It has two equilibria: a. disease free equilibrium 𝑃1 = (N, 0, 0) and an endemic equilibrium 𝑃2 = ( 1) ,. γ 𝛽. N 𝑅0. ,. b 𝛽. (𝑅0 -. (𝑅0 - 1) ). One major finding is that the epidemic is under control if the basic. reproduction number of the epidemic, is less than one, and the number of new infections produced during 𝜏 days is less than one too. It means that a delay in production of new infected cases, represented by 𝛽S(𝑡)I(t- 𝜏), could be a factor in the diseases dynamics.. 2.2 Dynamics In A System With Two Delays Dong[7] study the dynamical behavior of a tumor–immune system (T–IS) interaction model with two discrete delays, namely the immune activation delay for effector cells (ECs) and activation delay for helper T cells (HTCs). They establish the stability of two equilibria (tumor-free equilibrium 𝐸 − and immune control equilibrium 𝐸 + ) and the existence of Hopf bifurcations when two delays are used as the bifurcation parameters. Their results exhibit that both delays do not affect the stability of tumor-free equilibrium, but destabilize the immune-control equilibrium and cause periodic solutions. In their analysis three conditions are considered: (i) 𝜏1 = 𝜏2 = 0 (ii) 𝜏1 > 0, 𝜏2 = 0 (iii) 𝜏1 = 0, 𝜏2 > 0. For example while condition (ii) 𝜏1 > 0 and 𝜏2 = 0, the characteristic equation with eigenvalue λ becomes 7.

(16) D (λ) = 𝑃0 (λ) + 𝑃2 (λ) + (𝑃1 (λ) + 𝑃12 (λ)) 𝑒 −𝜆𝜏1 = 0. We rewrite this equation to 𝜆3 + 𝑎2 𝜆2 + 𝑎1 𝜆 + 𝑎0 + (𝑏2 𝜆2 + 𝑏1 𝜆 + 𝑏0 ) 𝑒 −𝜆𝜏1 =0, and take 𝜏1 as a parameter to investigate the Hopf bifurcation. The analytical results show that both delays 𝜏1 and 𝜏2 do not affect the stability of tumor-free equilibrium 𝐸 − , however they can destabilize the immune-control equilibrium 𝐸 + to undergo Hopf bifurcation. Time delay model is a tool to investigate the bifurcation phenomenon.. 8.

(17) CHAPTER 3 RESEARCH TOPICS AND METHODS Dengue virus can cause severe diseases such as dengue hemorrhagic fever and dengue shock syndrome. Antibody-dependent enhancement is a complex immune entity and it leads to morbidity and mortality. A dengue vaccine is developing with a vaccination schedule on 0, 6th and 12th month. Building effective antibody concentration needs an incubation period and the concentration decays as a natural course. In some circumstances high level of such antibody provides protection and the patient manifests a mild disease. An individual gets dengue infection under a low antibody level, he would get a more severe dengue infection. Differential equations with delays could be a useful tool to elucidate this sophisticated postvaccination dynamics.. 3.1 Modelling The SIR model with demography is a benchmark model. The basic model Rodrigues mentioned is a prototype which lacks of control or vaccination. In Elhia [5] model a control(vaccination) 𝑢 influence the susceptible and recovery subgroups. This model does not consider the decay of a vaccine potency and the proportion of effective vaccination response. In Perez [6] model they emphasize the delay term on infected individuals. In the vaccination viewpoints we suppose the delay influences susceptible 9.

(18) subpopulation would be appropriate. To investigate postvaccination dynamics we set up an initial model similar to the aforementioned two models as the following: 𝑑 𝑆(𝑡) 𝑑𝑡 𝑑 𝐼(𝑡) 𝑑𝑡 𝑑 𝑅(𝑡) 𝑑𝑡. = b–βS(t)I(t)-bS(t)-ρS(t-𝜏1 )+σS(t-𝜏2 ), = βS(t)I(t)-(b+γ+∈)I(t), = γI(t)-bR(t)+ρS(t-𝜏1 )-σS(t-𝜏2 ),. However, the vaccinated individuals could be mild or severe in the following infection. We modify this SIR model to “MHL model”; M represents moderate risk group, H represents high risk group and L represents low risk group. Several studies mention that those individuals get mild disease in the 3th and 4th episodes, we categorize these people as the low risk group based on this finding. For simplicity we assume the dengue-related mortality of M and L group is zero. Logistic equation is modelled into the M risk group and 3 delays is employed. The final model is modified as following: M(t) is moderate risk group for severe dengue disease, based on no enough amount of antibody or delayed production of antibody; H(t) is high risk group for severe dengue disease, based on low antibody level and antibody-dependent enhancement(ADE, due to previous heterotypic antibody, which could be recent or years ago); 10.

(19) L(t) is low risk group for severe dengue disease, based on coexisted multiple(≥ 2) serotypes of antibody; 𝑑 𝑀(𝑡) 𝑑𝑡 𝑑 𝐻(𝑡) 𝑑𝑡 𝑑 𝐿(𝑡) 𝑑𝑡. K. = rM(t)(1-. M(t) 𝐾. )–βM(t)H(t)-ρM(t)+σL(t-𝜏3 ). = α + βM(t-𝜏1 )H(t-𝜏1 )-(μ+γ+ε)H(t). (1). = γH(t)+ρM(t-𝜏2 )-σL(t)-μL(t). carrying capacity of population;. r. growth rate of population;. α. previous infected people with only one heterotypic dengue antibody;. β. transmission rate;. γ. recovery rate;. ε. disease-caused death rate;. μ. natural death rate;. ρ. the proportion of vaccinated susceptible individuals obtaining effective antibodies;. σ. the proportion of vaccinated susceptible individuals losing effective antibodies, mainly on account of antibody decay;. 𝜏1. the time span between contact of moderate risk group with high risk group and development of severe diseases of high risk group individuals;. 𝜏2. the time span of moderate risk group developing two heterotypic antibodies, among this group people have low or lowest risk for severe dengue disease; 11.

(20) 𝜏3. the time span of moderate risk group developing one monotypic antibody.. 3.2 Research Topics And Methods Our concerns include basic properties of this model like: (i) equilibria; (ii) local stability of steady states; (iii) existence of periodic solutions; (iv) Hopf bifurcations; (v) numerical simulation. Stability of equilibria are analyzed via linearizing the model and studying the relevant characteristic equation. All negative real parts of eigenvalues indicates the local stability of an equilibrium. About bifurcation phenomenon we focus on the Hopf bifurcation. Similar to the study of Dong [7], we use 𝜏1 , 𝜏2 and 𝜏3 as parameters to investigate the Hopf bifurcation.. 12.

(21) CHAPTER 4 ANALYSIS Preliminaries We denote the phase space of (1) by X=C ([−τ, 0], 𝑅+3 ), the space of continuous functions equipped with the norm ∥φ∥= 𝑠𝑢𝑝−τ≤θ≤0 |φ(θ)|, τ=max{𝜏1 , 𝜏2 , 𝜏3 }. Here, 𝑅+ :={ξ∈R|ξ≥0}. Write φ ≥ 0, for φ ∈X, if φ(θ) ≥ 0 for –τ ≤ θ ≤ 0 and denote φ =(ξ̂1 , ξ̂2 , ξ̂3 ) if 𝜑𝑖 (θ) = ξ𝑖 for θ∈[−τ, 0]. Herein, we consider the initial conditions φ ∈X, and the existence and uniqueness of the solution 𝜑𝑡 = Φt(φ) = (𝑀𝑡 , 𝐻𝑡 , 𝐿𝑡 ) of (1), as long as it exists, follow from the standard theory of functional differential equations [8]. The vector field defined by (1) is Lipschitzian in each compact set, so the initial value problem has a unique solution for all t ≥ 0 as long as the boundedness of solutions holds.. 4.1 Equilibrium We denote 𝐸 0 and 𝐸 ∗ the naïve equilibrium and post-exposure equilibrium respectively. (A) Naïve equilibrium (α = 0): A virgin region to Dengue virus would have its unique characteristics: no immunities 13.

(22) or increased susceptibilities to the various Dengue serotypes, and the existence of a single serotype throughout the epidemic. Under such a circumstance α equals to zero. The system (1) has following equilibrium: 𝐸10 (𝑀10 , 𝐻10 , 𝐿01 ) = ( 0, 0, 0 ), 𝐾(𝑟(𝜎+𝜇)−𝜇𝜌). 𝐸20 (𝑀20 , 𝐻20 , 𝐿02 ) = ( 𝐸30 (𝑀30 , 𝐻30 , 𝐿03 ) = ( 𝐻30 =. 𝑟(𝜎+𝜇). 𝐾𝜌(𝑟(𝜎+𝜇)−𝜇𝜌). , 0,. 𝑟(𝜎+𝜇)2. μ+γ+ε. 1. 𝛽. (𝜎+𝜇). , 𝐻30 ,. (γ𝐻30 +. 𝜌 𝛽. ),. (𝜇 + 𝛾 + 𝜀)) , where. (μ+γ+ε)(r(βK−μ−γ−ε)(𝜎+𝜇)−βμρK). ;. 𝛽 2 𝐾((𝜎+𝜇)(μ+γ+ε)−γσ). 𝐸20 exists if r>. 𝜇𝜌. ,. 𝜇+𝜎. 𝐸30 exists if r>. 𝛽𝜇𝜌𝐾 (𝜇+𝜎)(𝛽𝐾−𝜇−𝛾−𝜀). .. (B) Post-exposure equilibrium (α ≠ 0): The system (1) has post-exposure equilibrium while α ≠ 0, then 𝐸 ∗ = (𝑀∗ , 𝐻 ∗ , 𝐿∗ ) 1. = (𝛽𝐻 ∗ [(μ+γ+ε) 𝐻 ∗ -α] , 𝐻 ∗ ,. 1 σ+μ. [γ𝐻 ∗ +. 𝜌 𝛽𝐻 ∗. ((μ+γ+ε) 𝐻 ∗ -α)] ) .. 𝐻 ∗ is one of the roots of A𝐻 ∗3 + 𝐵𝐻 ∗2 + 𝐶𝐻 ∗ + 𝐷 = 0, where A = 𝛽 2 K[(σ+μ)(μ+γ+ε)- γσ], 14. (2).

(23) B = [-𝛽Kr(μ+γ+ε) +r(μ + γ + ε)2 − 𝛼𝛽 2 𝐾](σ+μ)+ 𝛽μρK(μ+γ+ε), C = 𝛼r(σ+μ)[ 𝛽K-2(μ+γ+ε)]- 𝛼𝛽μρK, D = 𝛼 2 r(σ + μ). There could be one to three positive roots of 𝐻 ∗ , by the Routh-Hurwitz criterion, the numbers of possible positive roots is listed in Table 1. In addition, the post-exposure equilibrium is positive if 𝛼. 𝐻∗ >. μ+γ+ε. >0.. 1. (BC-AD). A. B. +. +. +. 0. +. +. -. 2. +. -. + or -. 2. -. +. +. 1. -. +. -. 3. -. -. + or -. 1. 𝐵. Possible positive roots. Table 1. Numbers of positive roots of 𝐻 ∗. 4.2 Linearization In order to investigate the local stability of the equilibria of system (1), we linearize the 15.

(24) system and obtain the characteristic equation evaluated at an equilibrium 𝐸 ∗ (𝑀∗ , 𝐻 ∗ , 𝐿∗ ), given by the following determinant: − det(. 2𝑟 𝐾. 𝑀∗ − 𝛽𝐻 ∗ − 𝜌 + 𝑟 − 𝜆 𝛽𝐻 ∗ 𝑒 −𝜆𝜏1 𝜌𝑒 −𝜆𝜏2. −𝛽𝑀∗. 𝜎𝑒 −𝜆𝜏3. 𝛽𝑀∗ 𝑒 −𝜆𝜏1 − (μ + γ + ε) − 𝜆 𝛾. ) 0 −𝜎 − 𝜇 − 𝜆. = 0, that is D(λ, 𝜏1 ,. 𝜏2 , 𝜏3 ) = 𝑃0 (λ) + 𝑃1 (λ)𝑒 −𝜆𝜏1 +. 𝑃13 (λ)𝑒 −𝜆(𝜏1 +𝜏3) +. 𝑃23 (λ)𝑒 −𝜆(𝜏2 +𝜏3) + 𝑃123 (λ)𝑒 −𝜆(𝜏1 +𝜏2 +𝜏3) = 0, 𝑃0 (λ) = (λ +. 2𝑟 𝐾. (3). 𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟) (λ + 𝜎 + 𝜇) (λ + μ + γ + ε),. 𝑃1 (λ) = −𝛽𝑀∗ (λ + 𝜎 + 𝜇) (λ +. 2𝑟 𝐾. 𝑀∗ + 𝜌 − 𝑟),. 𝑃13 (λ) = -βγσ𝐻 ∗ , 𝑃23 (λ) = -ρσ(λ+μ + γ + ε), 𝑃123 (λ) = βρσ𝑀∗ .. 4.3 Local Stability And Hopf Bifurcation Now we analyze the local asymptotic stability of 𝐸 ∗ . Case (I) 𝜏1 ≠ 0, 𝜏2 =. 𝜏3 = 0.. D(λ, 𝜏1 , 𝜏2 , 𝜏3) = D(λ, 𝜏1 , 0, 0) 16.

(25) =𝑃0 (λ) + 𝑃1 (λ)𝑒 −𝜆𝜏1 + 𝑃13 (λ)𝑒 −𝜆𝜏1 + 𝑃23 (λ) + 𝑃123 (λ)𝑒 −𝜆𝜏1 = 0. (𝑃0 (λ) + 𝑃23 (λ)) + (𝑃1 (λ) + 𝑃13 (λ) + 𝑃123 (λ)) 𝑒 −𝜆𝜏1 = 0.. At 𝐸 ∗ (𝑀∗ , 𝐻 ∗ , 𝐿∗ ) , the characteristic equation (3) takes the form 𝜆3 + 𝑎2 𝜆2 + 𝑎1 𝜆. + 𝑎0 + (𝑏2 𝜆2 + 𝑏1 𝜆 + 𝑏0 ) 𝑒 −𝜆𝜏1 = 0,. (4). where 𝑎2 =. 2𝑟 𝐾. 𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟 + 𝜎 + 2𝜇 + γ + ε.. 2𝑟. 2𝑟. 𝑎1 = ( 𝐾 𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟)( 𝜎 + 𝜇) + ( 𝐾 𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟) ( 𝜇 + γ + ε) + (𝜎 + 𝜇)(𝜇 + γ + ε)-ρσ. 𝑎0 = (𝜇 + γ + ε)[(. 2𝑟 𝐾. 𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟)( 𝜎 + 𝜇)-ρσ].. 𝑏2 = - 𝛽𝑀∗ . 𝑏1 = - 𝛽𝑀∗ (𝜎 + 𝜇 +. 2𝑟 𝐾. 𝑀∗ + 𝜌 − 𝑟).. 2𝑟. 𝑏0 = - 𝛽𝑀∗ (𝜎 + 𝜇) ( 𝐾 𝑀∗ + 𝜌 − 𝑟) + βσ(ρ𝑀∗ − γ𝐻 ∗ ) .. If the time delay 𝜏1 is able to destabilize 𝐸 ∗ , this can only occur when characteristic roots cross the imaginary axis to the right as increasing 𝜏1 . Let ω > 0 and suppose λ= iω is a pure imaginary root of Eq. (4). Separating the real and imaginary parts, we have -𝜔3 + 𝑎1 ω = (𝑏0 − 𝑏2 𝜔2 ) sinω𝜏1 - 𝑏1 ω cosω𝜏1 , 17.

(26) 𝑎2 𝜔2 - 𝑎0 = (𝑏0 − 𝑏2 𝜔2 ) cosω𝜏1 + 𝑏1 ω sinω𝜏1 .. (5). Squaring and adding both above equations lead to the following equation 𝐹1 (ω) = 𝜔6 + ℎ2 𝜔4 + ℎ1 𝜔2 + ℎ0 = 0 .. (6). ℎ2 = 𝑎2 2 - 2𝑎1 - 𝑏2 2 . ℎ1 = 𝑎1 2 - 2𝑎0 𝑎2 - 𝑏1 2 + 2𝑏0 𝑏2 . ℎ0 = 𝑎0 2 - 𝑏0 2 . Letting 𝑣 = 𝜔2 in Eq. (6) gives the following simplified system G(𝑣) = 𝑣 3 + ℎ2 𝑣 2 + ℎ1 𝑣. + ℎ0 = 0 .. Then 𝑑 𝐺(𝑣) 𝑑𝑣. = 3𝑣 2 + 2ℎ2 𝑣 1. + ℎ1. = 3(𝑣 + 3 ℎ2 )2 -. 𝛥 3. ,. with two roots 𝑣1 and 𝑣2 if Δ > 0, where Δ = ℎ2 2 - 3ℎ1 . 𝑣1 =. −ℎ2 +√𝛥. 𝑣2 =. −ℎ2 −√𝛥. 3. 3. . .. We follow the method in [9] to analyze the existence of positive roots of Eq. (7) 𝐺 ′′ (𝑣) = 6 𝑣+2ℎ2 . 𝐺 ′′ (𝑣1 )= 2√𝛥 > 0, 𝑣1 is a local minimum. 𝐺 ′′ (𝑣2 )= -2√𝛥 < 0, 𝑣2 is a local maximum. 18. (7).

(27) If ℎ0 < 0, since lim 𝐺(𝑣) = ∞, then we obtain that Eq. (7) has at least one positive 𝑣→∞. root. When Δ ≤ 0,. 𝑑 𝐺(𝑣) 𝑑𝑣. ≥ 0. Under this condition 𝐺(𝑣) is monotonically increasing.. Thus, if ℎ0 ≥ 0 and Δ ≤ 0, then Eq. (7) has no positive roots. If ℎ0 ≥ 0, Eq. (7) has positive roots if and only if. Δ > 0, 𝑣1 > 0 and G(𝑣1 ) ≤ 0.. We summarize the conditions in solving Eq(7). (i) If ℎ0 < 0, then Eq. (7) has at least one positive root. (ii) If ℎ0 < 0 and Δ ≤ 0, then Eq. (7) has exactly one positive root. (iii) If ℎ0 ≥ 0 and Δ ≤ 0, then Eq. (7) has no positive roots. (iv) If ℎ0 ≥ 0, then Eq. (7) has positive roots if and only if Δ > 0, 𝑣1 > 0 and G(𝑣1 ) ≤ 0.. Suppose that Eq. (7) has three positive roots, which are defined by 𝑣̃1 , 𝑣̃2 and 𝑣̃3 , respectively. Then Eq. (7) has three positive roots 𝜔𝑘 = √𝑣̃𝑘 , k= 1, 2, 3. Solving Eq. (5) for 𝜏1 yields (𝑗). 𝜏1,𝑘 =. 1 𝜔𝑘. arcos. 2 (𝑏1 − 𝑎2 𝑏2 ) 𝜔𝑘4 + (𝑎0 𝑏2 + 𝑎2 𝑏0 −𝑎1 𝑏1 ) 𝜔𝑘 − 𝑎0 𝑏0 2 2 2 2 (𝑏2 𝜔𝑘 −𝑏0 ) + 𝑏1 𝜔𝑘. +. 2𝑗𝜋 𝜔𝑘. ,. Where k = 1, 2, 3; j = 0, 1, 2, .., and ±i𝜔𝑘 is a pair of pure imaginary roots of Eq. (6).. (𝑗). In the following, we obtain the transversality condition for the Hopf bifurcation at 𝜏1,𝑘 . 19.

(28) Differentiate Eq. (4) with respect to 𝜏1 , we get { 3𝜆2 + 2𝑎2 𝜆. + 𝑎1 + [ ( 2𝑏2 𝜆 + 𝑏1 )𝑒 −𝜆𝜏1 − (𝑏2 𝜆2 + 𝑏1 𝜆 + 𝑑𝜆. 𝑏0 )𝜏1 𝑒 −𝜆𝜏1 ]}𝑑𝜏. 1. = 𝜆 (𝑏2 𝜆2 + 𝑏1 𝜆 + 𝑏0 )𝑒 −𝜆𝜏1 . 𝜆3 + 𝑎2 𝜆2 + 𝑎1 𝜆. From Eq. (4), we have 𝑒 −𝜆𝜏1 = −. 𝑏2 𝜆2 +𝑏1. + 𝑎0. .. 𝜆 + 𝑏0. Thus 3𝜆2 + 2𝑎2 𝜆. 𝑑𝜆. (𝑑𝜏 )−1 = -. 𝜆(𝜆3. 1. + 𝑎2. 𝜆2. + 𝑎1. + 𝑎1 𝜆. + 𝑎0 ). +. 2𝑏2 𝜆 +𝑏1 𝜆(𝑏2 𝜆2 +𝑏1. 𝜆 + 𝑏0 ). -. 𝜏1 𝜆. .. From Eq. (5), we have (𝜔3 − 𝑎1 ω)2 + (𝑎2 𝜔2 − 𝑎0 )2 = (𝑏2 𝜔2 − 𝑏0 )2 + (𝑏1 ω)2 . 𝑑𝜆. (𝑗). Evaluating (𝑑𝜏 )−1 at 𝜏1,𝑘 (i.e. λ = iω, similarly for λ = - iω) and taking the real part, 1. we have 𝑑𝜆. Re [(𝑑𝜏 )−1] |𝜏 1. = Re [= = =. (𝑗) 1 = 𝜏1,𝑘. 3𝜆2 + 2𝑎2 𝜆 𝜆(𝜆3. + 𝑎2. 𝜆2. + 𝑎1. + 𝑎1 𝜆. + 𝑎0. ] + Re [𝜆(𝑏 ). 2𝑏2 𝜆 +𝑏1. 2𝜆. (3𝜔 2 − 𝑎1 − 2𝑎2 𝑖𝜔 )(𝜔 4 − 𝑎1 𝜔2 +𝑎2 𝑖𝜔3 −𝑎0 𝑖𝜔 ) (𝜔 4 −𝑎1 𝜔2 )2 + (𝑎2 𝜔3 −𝑎0 𝜔 )2. 2 +𝑏 1. −. ( 𝑏1 +2𝑏2 𝑖𝜔 )(𝑏1 𝜔2 +𝑏0 𝑖𝜔 −𝑏2 𝑖𝜔3 ). 3𝜔4 + 2(𝑎22 − 𝑏22 −2𝑎1 )𝜔2 + 𝑎12 −2𝑎0 𝑎2 − 𝑏12 + 2𝑏0 𝑏2 𝑏12 𝜔2 + (𝑏2 𝜔2 −𝑏0 )2. 𝐺1 (𝜔 2 ) 𝑏12 𝜔2 + (𝑏2 𝜔2 −𝑏0 )2. .. Therefore, sign {. 𝑑 𝑅𝑒(𝜆) 𝑑𝜏1. } |𝜏. (𝑗) 1 = 𝜏1,𝑘. 𝑑𝜆. = sign { Re[(𝑑𝜏 )−1]} |𝜏 1. (𝑗) 1 = 𝜏1,𝑘. 20. ]. 𝜆 + 𝑏0 ). (𝑏1 𝜔2 )2 +(𝑏2 𝜔3 −𝑏0 𝜔)2.

(29) = sign {𝐺1 (𝜔2 )} . So if 𝐺1 (𝜔2 ) ≠ 0, then the transversality condition holds.. Case (II) 𝜏2 ≠ 0, 𝜏1 = 𝜏3 = 0. D(λ, 𝜏1 , 𝜏2 , 𝜏3) = D(λ, 0, 𝜏2 , 0) = 𝑃0 (λ) + 𝑃1 (λ) + 𝑃13 (λ) + 𝑃23 (λ)𝑒 −𝜆𝜏2 + 𝑃123 (λ)𝑒 −𝜆𝜏2 = 0. (𝑃0 (λ) + 𝑃1 (λ) + 𝑃13 (λ)) + (𝑃23 (λ) + 𝑃123 (λ)) 𝑒 −𝜆𝜏2 = 0.. At 𝐸 ∗ (𝑀∗ , 𝐻 ∗ , 𝐿∗ ) , the characteristic equation (3) takes the form 𝜆3 + 𝑎2 𝜆2 + 𝑎1 𝜆. + 𝑎0 + (𝑏1 𝜆 + 𝑏0 ) 𝑒 −𝜆𝜏2 = 0.. (4)’. Where 2𝑟. 𝑎2 = ( 𝐾 − 𝛽)𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟 + 𝜎 + 2𝜇 + γ + ε. 2𝑟. 2𝑟. 𝑎1 = ( 𝐾 𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟)( 𝜎 + 𝜇) + ( 𝐾 𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟)( 𝜇 + γ + ε) 2𝑟. + (𝜎 + 𝜇)(𝜇 + γ + ε)- 𝛽𝑀∗ ( 𝐾 𝑀∗ + 𝜌 − 𝑟 + 𝜎 + 𝜇). 𝑎0 = ( 𝜇 + γ + ε )(. 2𝑟 𝐾. 2𝑟. 𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟 )( 𝜎 + 𝜇 ) - 𝛽𝑀∗ ( 𝐾 𝑀∗ + 𝜌 −. 𝑟)( 𝜎 + 𝜇)-βγσ𝐻 ∗ . 𝑏1 = - 𝜌𝜎. 21.

(30) 𝑏0 = 𝜌𝜎( 𝛽𝑀∗ − 𝜇 − γ − ε) .. If the time delay 𝜏2 is able to destabilize 𝐸 ∗ , this can only occur when characteristic roots cross the imaginary axis to the right as increasing 𝜏2 . Let ω > 0 and suppose λ= iω is a pure imaginary root of Eq. (4)’. Separating the real and imaginary parts, we have -𝜔3 + 𝑎1 ω = 𝑏0 sinω𝜏2 - 𝑏1 ω cosω𝜏2 , 𝑎2 𝜔2 - 𝑎0 = 𝑏0 cosω𝜏2 + 𝑏1 ω sinω𝜏2 .. (5)’. Squaring and adding both above equations lead to the following equation 𝐹2 (ω) = 𝜔6 + ℎ2 𝜔4 + ℎ1 𝜔2 + ℎ0 = 0 .. (6)’. ℎ2 = 𝑎2 2 - 2𝑎1 . ℎ1 = 𝑎1 2 - 2𝑎0 𝑎2 - 𝑏1 2 . ℎ0 = 𝑎0 2 - 𝑏0 2 . Letting v = 𝜔2 in Eq. (6)’ and repeating the same procedures as aforementioned 𝜏1 ≠ 0(𝜏2 = 𝜏3 = 0) gives the following simplified system G(𝑣) = 𝑣 3 + ℎ2 𝑣 2 + ℎ1 𝑣. + ℎ0 = 0 .. (7). Analyze G(𝑣) and summarize the conditions in solving Eq. (7) producing the same results as 𝜏1 ≠ 0(𝜏2 = 𝜏3 = 0).. 22.

(31) Suppose that Eq. (7) has three positive roots, which are defined by 𝑣̂1 , 𝑣̂2 and 𝑣̂3 , respectively. Then Eq. (7) has three positive roots 𝜔𝑘 = √𝑣̂𝑘 , k= 1, 2, 3. Solving Eq. (5)’ for 𝜏2 yields (𝑗). 𝜏2,𝑘 =. 1 𝜔𝑘. 4 𝑏1 𝜔𝑘 + ( 𝑎2 𝑏0 −𝑎1 𝑏1 ) 𝜔𝑘2 − 𝑎0 𝑏0 2 𝑏0 2 + 𝑏12 𝜔𝑘. arcos. +. 2𝑗𝜋 𝜔𝑘. ,. Where k = 1, 2, 3; j = 0, 1, 2, ….., and ±i𝜔𝑘 is a pair of pure imaginary roots of Eq. (𝑗). (6)’ with 𝜏2,𝑘 .. (𝑗). In the following, we obtain the transversality condition for the Hopf bifurcation at 𝜏2,𝑘 . Differentiate Eq. (4)’ with respect to 𝜏2 , we get 𝑑𝜆. {3𝜆2 + 2𝑎2 𝜆. + 𝑎1 + [ 𝑏1 𝑒 −𝜆𝜏2 − (𝑏1 𝜆 + 𝑏0 )𝜏2 𝑒 −𝜆𝜏2 ]}𝑑𝜏. 2. = 𝜆 (𝑏1 𝜆 + 𝑏0 )𝑒 −𝜆𝜏2 . From Eq. (4)’, we have 𝑒 −𝜆𝜏2 = −. 𝜆3 + 𝑎2 𝜆2 + 𝑎1 𝜆 𝑏1 𝜆 + 𝑏0. + 𝑎0. .. Thus 3𝜆2 + 2𝑎2 𝜆. 𝑑𝜆. (𝑑𝜏 )−1 = 2. 𝜆(𝜆3. + 𝑎2. 𝜆2. + 𝑎1. + 𝑎1 𝜆. + 𝑎0 ). +. 𝑏1 𝜆(𝑏1 𝜆 + 𝑏0 ). -. 𝜏2 𝜆. .. From Eq. (5)’, we have (𝜔3 − 𝑎1 ω)2 + (𝑎2 𝜔2 − 𝑎0 )2 = 𝑏0 2 + (𝑏1 ω)2 . 𝑑𝜆. (𝑗). Evaluating (𝑑𝜏 )−1 at 𝜏2,𝑘 (i.e. λ = iω, similarly for λ = -iω) and taking the real part, 2. we have 𝑑𝜆. Re [(𝑑𝜏 )−1]|𝜏 2. (𝑗) 2 = 𝜏2,𝑘. 23.

(32) = Re [= = =. 3𝜆2 + 2𝑎2 𝜆 𝜆(𝜆3. + 𝑎2. 𝜆2. + 𝑎1. + 𝑎1 𝜆. + 𝑎0. ] + Re [𝜆(𝑏 ). 1. (3𝜔 2 − 𝑎1 − 2𝑎2 𝑖𝜔 )(𝜔 4 − 𝑎1 𝜔2 +𝑎2 𝑖𝜔3 −𝑎0 𝑖𝜔 ) (𝜔 4 −𝑎1 𝜔2 )2 + (𝑎2 𝜔3 −𝑎0 𝜔 )2. 𝑏1. ]. 𝜆 + 𝑏0 ). −. 𝑏1 (𝑏1 𝜔2 +𝑏0 𝑖𝜔) (𝑏1 𝜔2 )2 +(𝑏0 𝜔)2. 3𝜔4 + 2(𝑎22 −2𝑎1 )𝜔2 + 𝑎12 −2𝑎0 𝑎2 − 𝑏12 𝑏12 𝜔2 + 𝑏02. 𝐺2 (𝜔 2 ). .. 𝑏12 𝜔2 + 𝑏0 2. Therefore, sign {. 𝑑 𝑅𝑒(𝜆) 𝑑𝜏2. }|𝜏. (𝑗) 2 = 𝜏2,𝑘. 𝑑𝜆 −1. = sign { Re[(𝑑𝜏 ) ]}|𝜏 2. (𝑗) 2 = 𝜏2,𝑘. = sign {𝐺2 (𝜔2 )} . So if 𝐺2 (𝜔2 ) ≠ 0, then the transversality condition holds.. Case (III) 𝜏3 ≠ 0, 𝜏1 = 𝜏2 = 0. D(λ, 𝜏1 , 𝜏2 , 𝜏3) = D(λ, 0, 0, 𝜏3 ) = 𝑃0 (λ) + 𝑃1 (λ) + 𝑃13 (λ)𝑒 −𝜆𝜏3 + 𝑃23 (λ) 𝑒 −𝜆𝜏3 + 𝑃123 (λ)𝑒 −𝜆𝜏3 = 0. (𝑃0 (λ) + 𝑃1 (λ)) + (𝑃13 (λ) + (𝑃23 (λ) + 𝑃123 (λ)) 𝑒 −𝜆𝜏3 = 0.. At 𝐸 ∗ (𝑀∗ , 𝐻 ∗ , 𝐿∗ ) , the characteristic equation (3) takes the form 𝜆3 + 𝑎2 𝜆2 + 𝑎1 𝜆. + 𝑎0 + (𝑏1 𝜆 + 𝑏0 ) 𝑒 −𝜆𝜏3 = 0. 24. (4)’’.

(33) Where 2𝑟. 𝑎2 = ( 𝐾 − 𝛽)𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟 + 𝜎 + 2𝜇 + γ + ε. 2𝑟. 𝑎1 = ( 𝐾 𝑀∗ + 𝜌 − 𝑟)(-𝛽𝑀∗ + 𝜇 + γ + ε) + 𝛽𝐻 ∗ (𝜇 + γ + ε) + ( 𝜎 + 𝜇) 2𝑟. [( 𝐾 − 𝛽)𝑀∗ + 𝛽𝐻 ∗ + 𝜌 − 𝑟 + 𝜇 + γ + ε]. 𝑎0 = ( 𝜎 + 𝜇) [(. 2𝑟 𝐾. 𝑀∗ + 𝜌 − 𝑟) (−𝛽𝑀∗ + 𝜇 + γ + ε) +β𝐻 ∗ (𝜇 + γ + ε)].. 𝑏1 = - 𝜌𝜎. 𝑏0 = 𝛽𝜎(𝜌𝑀∗ − γ𝐻 ∗ ) − 𝜌𝜎( 𝜇 + γ + ε) .. If the time delay 𝜏3 is able to destabilize 𝐸 ∗ , this can only occur when characteristic roots cross the imaginary axis to the right as increasing 𝜏3 . Let ω > 0 and suppose λ= iω is a pure imaginary root of Eq. (4)’’. Separating the real and imaginary parts, we have -𝜔3 + 𝑎1 ω = 𝑏0 sinω𝜏3 - 𝑏1 ω cosω𝜏3 , 𝑎2 𝜔2 - 𝑎0 = 𝑏0 cosω𝜏3 + 𝑏1 ω sinω𝜏3 .. (5)’’. Squaring and adding both above equations lead to the following equation 𝐹3 (ω) = 𝜔6 + ℎ2 𝜔4 + ℎ1 𝜔2 + ℎ0 = 0 . ℎ2 = 𝑎2 2 - 2𝑎1 . ℎ1 = 𝑎1 2 - 2𝑎0 𝑎2 - 𝑏1 2 . 25. (6)’’.

(34) ℎ0 = 𝑎0 2 - 𝑏0 2 . Letting v = 𝜔2 in Eq. (6)’’ and repeating the same procedures as aforementioned 𝜏1 ≠ 0 (𝜏2 = 𝜏3 = 0) and 𝜏2 ≠ 0 (𝜏1 = 𝜏3 = 0) gives the following simplified system G(𝑣) = 𝑣 3 + ℎ2 𝑣 2 + ℎ1 𝑣. + ℎ0 = 0 .. (7). Analyze G(𝑣) and summarize the conditions in solving Eq. (7) producing the same results as 𝜏2 ≠ 0 (𝜏1 = 𝜏3 = 0). Suppose that Eq. (7) has three positive roots, which are defined by 𝑣̌1 , 𝑣̌2 and 𝑣̌3 , respectively. Then Eq. (7) has three positive roots 𝜔𝑘 = √𝑣̌𝑘 , k= 1, 2, 3. Solving Eq. (5)’’ for 𝜏3 yields (𝑗). 𝜏3,𝑘 =. 1 𝜔𝑘. 4 𝑏1 𝜔𝑘 + ( 𝑎2 𝑏0 −𝑎1 𝑏1 ) 𝜔𝑘2 − 𝑎0 𝑏0 2 𝑏0 2 + 𝑏12 𝜔𝑘. arcos. +. 2𝑗𝜋 𝜔𝑘. ,. Where k = 1, 2, 3; j = 0, 1, 2, ….., and ±i𝜔𝑘 is a pair of pure imaginary roots of Eq. (𝑗). (6)’’ with 𝜏3,𝑘 . (𝑗). In the following, we obtain the transversality condition for the Hopf bifurcation at 𝜏3,𝑘 . Differentiate Eq. (4)’’ with respect to 𝜏3 , we get 𝑒 −𝜆𝜏3 = −. 𝜆3 + 𝑎2 𝜆2 + 𝑎1 𝜆 𝑏1 𝜆 + 𝑏0. + 𝑎0. .. Thus 3𝜆2 + 2𝑎2 𝜆. 𝑑𝜆. (𝑑𝜏 )−1 = 3. 𝑑𝜆. 𝜆(𝜆3. + 𝑎2. 𝜆2. + 𝑎1. + 𝑎1 𝜆. + 𝑎0 ). (𝑗). +. 𝑏1 𝜆(𝑏1 𝜆 + 𝑏0 ). -. 𝜏3 𝜆. .. Evaluating (𝑑𝜏 )−1 at 𝜏3,𝑘 (i.e. λ = iω and similarly for λ = - iω) and taking the real 3. part, we have 𝑑𝜆. Re [(𝑑𝜏 )−1] |𝜏 3. (𝑗) 3 = 𝜏3,𝑘. 26.

(35) = =. 3𝜔4 + 2(𝑎22 −2𝑎1 )𝜔2 + 𝑎12 −2𝑎0 𝑎2 − 𝑏12 𝑏12 𝜔2 + 𝑏02. 𝐺3 (𝜔 2 ). .. 𝑏12 𝜔2 + 𝑏0 2. Therefore, sign {. 𝑑 𝑅𝑒(𝜆) 𝑑𝜏3. }|𝜏. (𝑗) 3 = 𝜏3,𝑘. 𝑑𝜆. = sign { Re[(𝑑𝜏 )−1]} |𝜏 3. (𝑗) 3 = 𝜏3,𝑘. = sign {𝐺3 (𝜔2 )} . So if 𝐺3 (𝜔2 ) ≠ 0, then the transversality condition holds.. 27.

(36) CHAPTER 5 NUMERICAL SIMULATION In clinical practice 𝜏1 , 𝜏2 and 𝜏3 means the time an individual obtains antibody or the incubation period to get disease. All are not equal to zero. The clinical set of parameter values could be as K=500, r=5, α=0.002[10], β= 1[11], γ=0.15[12], ε=0.01[13], μ=0.0073[14], ρ=0.75[15], σ=0.01, 𝜏1 =8[16, 17], 𝜏2 =14[15], 𝜏3 =14. The range of parameters is listed in Table 2. 0.00187 is the prevalence of dengue in 2015[10], we might assume α=0.002 based on the fact of increased prevalence in an endemic area. If 𝜏1 = intrinsic incubation period(in host) + extrinsic intubation period(in virus), it takes a clinical meaning. We apply software Matlab (version R2017b) with program dde23 to simulate numerically.. 28.

(37) Minimum. Maximum. Reference. K carrying capacity r growth rate. 0.000004 (Brasil). α previous one heterotype. 0.00187 (Taiwan 2015 prevalence). β. 0.12-0.24. transmission rate. 0.0045. 0.02 [10]. 100-300. 200-400. /year. /year. [11]. γ recovery rate. 0.00370.0167 (in 2-9 mo). 0.146/day 0.274/day (100/year). 0.329. [12]. ε disease-caused death rate. 2× 10−7 , 2.92× 10−8. 0.005-0.01 (0-1/day). 0.26. [13]. μ natural death rate. 0.0073 (Taiwan). ρ 0.54 vaccine effective % (4 types, age (0.50, or serological status). [14]. 0.60 0.34, (0.66 > 9y,. 0.75, 0.77). 0.45 < 9 y) 0.62-0.79. σ vaccine failure %. Immunity lose in 2-9 months (not vaccine). 𝜏1 MH contact time. 8(2+6). 𝜏2 M two heterotype time. 14 (after dose 164. 𝜏3 M one monotype time. 14. 1). 29. [15] 0.68,. 0.78, 0.90). Immunity lose 0.0055/d (not vaccine) [16, 17]. dose 2). Table 2. Reference values of parameters. 0.75 (0.70,. (after 344 dose 3). (after [15].

(38) We choose various sets of parameter values to perform numerical simulation. (I) 𝜏1 , 𝜏2 , 𝜏3 ≠ 0; [𝜏1 , 𝜏2 , 𝜏3 ]=[8,14,14];. [M(0), H(0), L(0)]=[1,2,3];. K=500, r=2, α=0.01, β= 0.1, γ=0.8, ε=0.11, μ=0.0275, ρ=0.34, σ=0.01. (a) 𝑥1 = M = moderate risk group, t = days,. Figure 1. [𝜏1 , 𝜏2 , 𝜏3 ]=[8,14,14]; Dynamics of moderate risk group From Figure 1-3, it shows that when this set of parameters exists, a stable periodic solution with a simple shape appears.. 30.

(39) (b) 𝑥2 = H = high risk group, t = days,. Figure 2. [𝜏1 , 𝜏2 , 𝜏3 ]=[8,14,14]; Dynamics of high risk group. 31.

(40) (c) 𝑥3 = L = low risk group, t = days,. Figure 3. [𝜏1 , 𝜏2 , 𝜏3 ]=[8,14,14]; Dynamics of low risk group. 32.

(41) (II) 𝜏1 ≠ 0, 𝜏2 =0, 𝜏3 =0; [𝜏1 , 𝜏2 , 𝜏3 ] =[8, 0, 0]; [M(0), H(0), L(0)] =[1,2,3]; K=500, r=2, α=0.01, β= 0.1, γ=0.3, ε=0.11, μ=0.0275, ρ=0.34, σ=0.01. (a) 𝑥1 = M = moderate risk group, t = days,. Figure 4. [𝜏1 , 𝜏2 , 𝜏3 ]=[8, 0, 0]; Dynamics of moderate risk group Under this set of parameters the coefficients of equation (2) and (7) are: A= 0.067031, B=-1.423613, C=0.032169, D=0.0000075; 𝑀∗ = 4.370286, 𝐻 ∗ = 21.2155(we take the largest, the other two are 0.0229, -0.0002), 𝑀∗ = 4.370286, 𝐿∗ = 209.3479; ℎ2 = 0.255143, ℎ1 = -0.451383, ℎ0 = -0.000429, Δ= 1.419246, 𝑣1 = 0.312059, G(𝑣1 )=-0.086052. By the judgement criteria (ii) of equation (7) ℎ0 < 0 is a requirement to obtain a positive ω. Therefore there is periodic solution under this set of parameter values. From Figure 4-6, it shows that when this set of parameters exists, a stable periodic solution appears. 33.

(42) (b) 𝑥2 = H = high risk group, t = days,. Figure 5. [𝜏1 , 𝜏2 , 𝜏3 ]=[8, 0, 0]; Dynamics of high risk group. (c) 𝑥3 = L = low risk group, t = days,. Figure 6. [𝜏1 , 𝜏2 , 𝜏3 ]=[8, 0, 0]; Dynamics of low risk group. 34.

(43) (III) 𝜏2 ≠ 0, 𝜏1 =0, 𝜏3 =0; [𝜏1 , 𝜏2 , 𝜏3 ] =[0, 14, 0]; [M(0), H(0), L(0)] =[1,2,3]; K=500, r=3, α=0.01, β= 0.12, γ=0.146, ε=0.26, μ=0.0275, ρ=0.34, σ=0.01. (a) 𝑥1 = M = moderate risk group, t = days,. Figure 7. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 14, 0]; Dynamics of moderate risk group Under this set of parameters the coefficients of equation (2) and (7) are: A= 0.106533, B=-2.66449, C=0.0609146, D=0.00001125; 𝑀∗ = 3.609165, 𝐻 ∗ = 24.9881(we take the largest, the other two are 0.0231, -0.0002), 𝐿∗ = 130.0101; ℎ2 = -2.450122, ℎ1 = 1.6871857, ℎ0 = 0.001965, Δ= 0.941539, 𝑣1 = 1.1401503, G(𝑣1 )= 0.2227226. 35.

(44) By the judgement criteria (iv) of equation (7) G(𝑣1 ) ≤ 0 is a requirement to obtain a periodic solution. Therefore there is no periodic solution under this set of parameter values. From Figures 7-9 it shows that under this set of parameters the solution of system (1) tends to the equilibrium 𝑀∗ , 𝐻 ∗ and 𝐿∗ respectively.. (b) 𝑥2 = H = high risk group, t = days,. Figure 8. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 14, 0]; Dynamics of high risk group. 36.

(45) (c) 𝑥3 = L = low risk group, t = days,. Figure 9. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 14, 0]; Dynamics of low risk group. 37.

(46) (IV) 𝜏3 ≠ 0, 𝜏1 =0, 𝜏2 =0; [𝜏1 , 𝜏2 , 𝜏3 ] =[0, 0, 14]; [M(0), H(0), L(0)] =[1,2,3]; K=500, r=3, α=0.01, β= 0.12, γ=0.146, ε=0.26, μ=0.0275, ρ=0.34, σ=0.01. (a) 𝑥1 = M = moderate risk group, t = days,. Figure 10. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 0, 14]; Dynamics of moderate risk group Under this set of parameters the coefficients of equation (2) and (7) are: A= 0.106533, B=-2.66449, C=0.0609146, D=0.00001125; 𝑀∗ = 3.609165, 𝐻 ∗ = 24.9881(we take the largest, the other two are 0.0231, -0.0002), 𝑀∗ = 3.609165, 𝐿∗ = 130.0101; ℎ2 = -56.84794, ℎ1 = 812.89781, ℎ0 = 0.0023531, Δ= 792.99522, 𝑣1 = 28.336038, G(𝑣1 )= 141.21162. 38.

(47) By the judgement criteria (iv) of equation (7) G(𝑣1 ) ≤ 0 is a requirement to obtain a periodic solution. Therefore there is no periodic solution under this set of parameter values. From Figures 10-12 it shows that under this set of parameters the solution of system (1) tends to the equilibrium 𝑀∗ , 𝐻 ∗ and 𝐿∗ respectively.. (b) 𝑥2 = H = high risk group, t = days,. Figure 11. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 0, 14]; Dynamics of high risk group. 39.

(48) (c) 𝑥3 = L = low risk group, t = days,. Figure 12. [𝜏1 , 𝜏2 , 𝜏3 ]=[0, 0, 14]; Dynamics of low risk group. 40.

(49) CHAPTER 6 DISCUSSION In this study we adopt differential equation with delays, which is employed in physics, chemistry, biology, medicine and social sciences. The delay provides time lag between drug or vaccine injection and their effects. This is one of the medical examples, therefore there is more flexible and diversified aspects than ordinary differential equations to explain biological phenomena. We use a model similar to the SIR model, but we transform to a MHL model. This model states the risk stages and it is different from the SIR model on the totalized sum. In our model M+H+L ≠ 1 due to there is disease-caused death. In our model vaccine efficacy and failure rate(ρ, σ) are analyzed. The delays also provide a reference to adjust the optimal vaccination timing. If the initial timing of vaccination is month 0, 6th and 12th, it could be adjusted to 0, 1st and 6th or others to get more effective disease control. We analyze its equilibria, local stability and the Hopf bifurcation. We simulated periodic solution in a numerical way. In a clinical scenario 𝜏1 , 𝜏2 and 𝜏3 mean the time an individual obtains antibody or the incubation period to get disease. All are not equal to zero and they would take days to happen. The set of parameter values could be restricted and varied in a large range, e.g. α is the previous infected people with only one heterotypic dengue antibody. It is about 0.002 in some cities in Taiwan, but could 41.

(50) be higher in Thailand or Cuba. While its value equals to zero the equilibrium is a naïve equilibrium and we can analyze the epidemical dynamics in a virgin region to dengue virus invasion. Vaccine efficacy to individual serotype of dengue influence the dynamics of epidemics because it varies from 0.34 to 0.90 with a mean about 0.75. While an epidemic occurs the insult serotype and the vaccine efficacy can lead to a minor or large scale. Due to the diversity of parameters aspect the periodic solution cannot show up always. This is a common biological phenomenon. We take 𝜏1 = intrinsic incubation period(in host) + extrinsic intubation period(in virus), it takes a clinical meaning to realize the disease dynamics. Although it might not be predicted in a periodic way for each delay, the numerical simulation gives us an approximate number of ultimate high risk group. The number is important to clinical physicians to alert.. 42.

(51) REFERENCE [1] Rigau-Perez, JG., Clark GG, Gubler DJ, et al. Dengue and dengue haemorrhagic fever. Lancet 1998; 352: 971–77. [2] Hadinegoro SR, Arredondo‑García JL, Capeding MR, et al. Efficacy and Longterm Safety of a Dengue Vaccine in Regions of Endemic Disease. New England Journal of Medicine 2015(373); 13: 1195-1206. [3] Anderson KB, Gibbons RV, Cummings DAT, et al. A shorter time interval between first and second dengue infections is associated with protection from clinical illness in a school based cohort in Thailand. Journal of Infectious Diseases 2014; 209 (1 February): 360-368. [4] Rodrigues HS. Application of SIR epidemiological model: new trends. International journal of applied mathematics and informatics. 2016; 10: 92-97. [5] Elhia M, Rachik M, and Benlahmar E. Optimal control of an SIR model with delay in state and control variables. ISRN Biomathematics 2013; article ID 403549, 1-7. [6] Perez AS, Rodriguez HA, Sanchez TN, et al. A theoretical model for the dengue epidemic using delayed differential equation: numerical approaches. IWANN 2009; Part I, LNCS 5517: 893-900. [7] Dong Y, Huang G, Miyazaki R, et al. Dynamics in a tumor immune system with time delays. Applied Mathematics and Computation 2015; 252: 99-113. 43.

(52) [8] Hale JK, Verduyn Lunel SM, Introduction to Functional Differential Equations, Springer-Verlag, New York, 1993. [9] Yongli Song, Sanling Yuan. Bifurcation analysis in a predator-prey system with time delay. Nonlinear analysis: Real world application 2006; 7: 265-284. [10] Hsu JC, Hsieh CL, Lu CY. Trend and geographic analysis of the prevalence of dengue in Taiwan, 2010–2015. International Journal of Infectious Diseases 2017; 54:. 43–49.. [11] Ferguson NM, Donnelly CA, Anderson RM. Transmission dynamics and epidemiology of dengue: insights from age-stratified sero-prevalence surveys. Philosophical Transactions of the Royal Society B: Biological Sciences 1999; 354: 757-768. [12] Adams B, Boots M. Modelling the relationship between antibody-dependent enhancement and immunological distance with application to dengue. Journal of Theoretical Biology 2006; 242: 337–346. [13] Durham DP, Ndeffo Mbah ML, Medlock J, et al. Dengue dynamics and vaccine cost-effectiveness in Brazil. Vaccine 2013; 31(37): 3957–3961. [14]. Report. of. Ministry. of. Health. and. Wealth,. ROC,. Taiwan.. https://www.mohw.gov.tw/cp-16-33598-1.html. Accessed on May 27, 2018. [15] Yang Y, Meng Y, Halloran ME, et al. Dependency of Vaccine Efficacy on 44.

(53) Preexposure and Age: A Closer Look at a Tetravalent Dengue Vaccine. Clinical Infectious Diseases 2018; 66(2): 178–184. [16] Halstead SB. Dengue. Lancet 2007; 370: 1644–52. [17] Watts DM, Burke DS, Harrison BA, et al. Effect of temperature on the vector efficiency of Aedes aegypti for dengue 2 virus. The American journal of tropical medicine and hygiene 1987; 36: 143-152.. 45.

(54) 46.

(55)