Open Collections

UBC Theses and Dissertations

Vaccination models in infectious diseases Bai, Fan 2016

Media
24-ubc_2016_september_bai_fan.pdf [ 854.52kB ]
JSON: 24-1.0305652.json
JSON-LD: 24-1.0305652-ld.json
RDF/XML (Pretty): 24-1.0305652-rdf.xml
RDF/JSON: 24-1.0305652-rdf.json
Turtle: 24-1.0305652-turtle.txt
N-Triples: 24-1.0305652-rdf-ntriples.txt
Original Record: 24-1.0305652-source.json
Full Text
24-1.0305652-fulltext.txt
Citation
24-1.0305652.ris

Full Text

`Vaccination models in infectious diseasesbyFan BaiB.Sc., Nankai University, 2006M.Sc., Memorial University of Newfoundland, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Mathematics)The University of British Columbia(Vancouver)June 2016c© Fan Bai, 2016AbstractVaccination is the most effective method of preventing the spread of infectiousdiseases. In this thesis, we develop and apply mathematical models to study vacci-nation. The thesis consists of three main chapters.1, In deciding whether to be vaccinated before the outbreak of the epidemic,people need to consider the risk from vaccination and the probability of be-ing infected. The decision of one individual is indirectly influenced by thedecisions of all other individuals. Because the vaccine coverage levels aredetermined by all decisions of all individuals. We apply game theory to thisscenario, to predict the expected vaccine coverage level. We construct andanalyze four vaccination games. The novelty of this part of work is we intro-ducing the replicator equations to describe the evolutionary processes. Forall games, we are able to predict the vaccine coverage levels accurately.2, In order to prove the uniqueness of Nash equilibrium in games, it is crucialto prove the attack ratios are decreasing functions of the vaccine coveragelevels. We are able to obtain complete results for the cases of homogeneousmixing population. Only partial results are obtained for the case of hetero-geneous mixing population.3, Some insights into the dynamics of malaria infection are obtained. We pro-pose two new malaria models. For delayed malaria transmission models, wecalculate the basic reproduction numberR0, the disease-free equilibrium andthe possible endemic equilibrium. We also analyze the stabilities of the equi-libria. For malaria vaccination model, we calculate the control reproductionnumberRc and the disease-free equilibrium. The threshold of eradication isiidiscussed. We find that, under certain circumstances, the disease of malariacan be eradicated.iiiPrefaceThe introductions in sections 2.1, 2.3.1, 2.3.3 and 4.1 provide background andmotivations. No original results are presented.The research presented in game theory and epidemiology (sections 2.2, 2.3.2),attack ratios and vaccine coverage levels (chapter 3) and two malaria transmissionmodels (sections 4.2 and 4.3) are carried out by the author, under the supervisionof prof. Fred Brauer.The paper ”Uniqueness of Nash Equilibrium in Vaccination Games” has beenaccepted by the Journal of Biological Dynamics.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction and thesis outline . . . . . . . . . . . . . . . . . . . . . 12 Game theory and vaccination . . . . . . . . . . . . . . . . . . . . . . 52.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1.1 Game theory . . . . . . . . . . . . . . . . . . . . . . . . 62.1.1.1 Some important definitions in classic game theory 62.1.1.2 Definitions in evolutionary game theory . . . . 72.1.2 Deterministic epidemic models . . . . . . . . . . . . . . . 92.1.2.1 The simple KermackMcKendrick model and thebasic reproduction numberR0 . . . . . . . . . 92.1.2.2 The final size relation . . . . . . . . . . . . . . 10v2.1.2.3 The next generation matrix approach and the ba-sic reproduction numberR0 . . . . . . . . . . . 122.1.2.4 The deterministic epidemic models with hetero-geneous mixing . . . . . . . . . . . . . . . . . 142.1.2.5 The age-of-infection epidemic model . . . . . . 162.2 Vaccination games . . . . . . . . . . . . . . . . . . . . . . . . . 182.2.1 Vaccination games in homogeneous mixing population . . 192.2.1.1 Vaccination game with perfect vaccine . . . . . 192.2.1.2 Vaccination game with imperfect vaccine . . . . 222.2.2 Vaccination games in heterogeneous mixing population . . 252.2.2.1 Vaccination game with perfect vaccine . . . . . 252.2.2.2 Vaccination game with imperfect vaccine . . . . 282.3 Disease eradication . . . . . . . . . . . . . . . . . . . . . . . . . 322.3.1 Herd immunity: introduction and mathematical formulation 322.3.2 Herd immunity in heterogeneous mixing population . . . 372.3.2.1 Perfect vaccine . . . . . . . . . . . . . . . . . . 382.3.2.2 Imperfect vaccine . . . . . . . . . . . . . . . . 422.3.2.3 Possible extensions . . . . . . . . . . . . . . . 472.3.3 Economic considerations for the eradication . . . . . . . . 472.4 Conclusion of this chapter . . . . . . . . . . . . . . . . . . . . . 493 Mathematical analysis of vaccination models . . . . . . . . . . . . . 513.1 Attack ratios and vaccine coverage levels . . . . . . . . . . . . . 513.1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . 513.1.2 Analysis in homogeneous mixing population . . . . . . . 533.1.2.1 Case of perfect vaccine . . . . . . . . . . . . . 533.1.2.2 Case of imperfect vaccine . . . . . . . . . . . . 553.1.3 Analysis in heterogeneous mixing population . . . . . . . 603.1.3.1 Case of perfect vaccine . . . . . . . . . . . . . 603.1.3.2 Case of imperfect vaccine . . . . . . . . . . . . 643.2 Some extensions to the age-of-infection models . . . . . . . . . . 713.3 Conclusion of this chapter . . . . . . . . . . . . . . . . . . . . . 74vi4 Malaria vaccination model . . . . . . . . . . . . . . . . . . . . . . . 764.1 General frame of modelling the vector-borne disease . . . . . . . 764.1.1 Introduction to vector borne diseases . . . . . . . . . . . 764.1.2 Mathematical modelling of malaria . . . . . . . . . . . . 794.2 Dynamical analysis of a delayed malaria transmission model . . . 834.2.1 Previous work . . . . . . . . . . . . . . . . . . . . . . . 844.2.2 A modified delayed malaria transmission model . . . . . . 854.2.2.1 The mathematical model . . . . . . . . . . . . 864.2.2.2 The basic reproduction number R0 and the en-demic equilibrium . . . . . . . . . . . . . . . . 874.2.2.3 Behavior at the disease-free equilibrium . . . . 904.2.2.4 Behavior at the endemic equilibrium . . . . . . 914.2.2.5 Numerical simulations . . . . . . . . . . . . . . 954.3 Dynamical analysis of a malaria vaccination model . . . . . . . . 984.3.1 Some facts of the malaria vaccine RTS, S . . . . . . . . . 984.3.2 A malaria vaccination model . . . . . . . . . . . . . . . . 994.3.3 The control reproduction numberRc . . . . . . . . . . . 1024.3.3.1 The next generation matrix approach . . . . . . 1044.3.3.2 Another method of approach . . . . . . . . . . 1064.3.4 The stability of the disease-free equilibrium . . . . . . . . 1074.3.5 The threshold of eradication . . . . . . . . . . . . . . . . 1084.4 Conclusion of this chapter . . . . . . . . . . . . . . . . . . . . . 1105 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1125.1 Epidemiological game theory . . . . . . . . . . . . . . . . . . . . 1125.2 Attack ratios and vaccine coverage levels . . . . . . . . . . . . . 1135.3 Contact patterns in heterogeneous mixing population . . . . . . . 1145.4 Some work in malaria transmission models . . . . . . . . . . . . 1145.5 Backward bifurcation analysis for malaria vaccination models . . 116Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122A Similarity of forms of final size relations . . . . . . . . . . . . . . . . 128viiB Asymptotic stability of the disease-free equilibrium . . . . . . . . . 134viiiList of TablesTable 2.1 R0 and herd immunity threshold . . . . . . . . . . . . . . . . 36Table 4.1 Table 1 of vector-borne diseases (cite from [36]) . . . . . . . . 77Table 4.2 Table 2 of vector-borne diseases (cite from [36]) . . . . . . . . 78Table 4.3 Table 3 of vector-borne diseases (cite from [36]) . . . . . . . . 79Table 4.4 Table 4 of vector-borne diseases (cite from [36]) . . . . . . . . 80Table 4.5 Table of variables and parameters for Ross-Macdonald model . 82Table 4.6 Table of variables and parameters for delayed malaria model . 87ixList of FiguresFigure 2.1 SVIR model . . . . . . . . . . . . . . . . . . . . . . . . . . . 24Figure 4.1 The life cycle of malaria parasite (cite from CDC) . . . . . . . 81Figure 4.2 The endemic equilibrium is asymptotically stable . . . . . . . 96Figure 4.3 The endemic equilibrium is asymptotically stable . . . . . . . 97Figure 4.4 The disease-free equilibrium is asymptotically stable . . . . . 98Figure 4.5 The disease-free equilibrium is asymptotically stable . . . . . 99Figure 4.6 Malaria vaccination model . . . . . . . . . . . . . . . . . . . 102xAcknowledgmentsFirst and foremost I would like to thank my advisor Professor Fred Brauer, whointroduced me into this fascinating area of mathematical epidemiology. I am sograteful for his guidance and support throughout the years. It is my honour tohave the opportunity to learn from Fred. This thesis would not have been possiblewithout him.I would like to thank my beloved wife, Xin Zhao, for everything. I am sofortunate to have her love, support and care.I would like to thank my parents for supporting me to chase my dream, foreducating me to be a brave and honest man.I would like to thank Professor Daniel Coombs and Professor Leah Keshet forserving on my advisory committee. I am grateful to Leah for her help and support.I would like to thank Professor David Earn from McMaster University for writ-ing the letter of recommendation for me, for attending my comprehensive exami-nation, and for the valuable advice on my research. I would like to thank Profes-sor Carlos Castillo-Chavez from Arizona State University for reading some of mymanuscripts, for discussing my research.I would like to thank Professor Gordon Slade, Professor Keqin Liu for so manygreat suggestions.I would like to thank Lee Yupitun and Roseann Kinsey for excellent work inhelping us graduate students in Math Department throughout these years. I alsowould like to thank my friends and all classmates at UBC. This is a long list.These years we spent together at UBC will be one of my most precious memories.xiDedicationTo my Family.xiiChapter 1Introduction and thesis outlineMathematical epidemiology has a long history. In 1760, Daniel Bernoulli formu-lated the smallpox model. Much of the fundamental theory was developed between1900 and 1935. One major difference between mathematical epidemiology andmost natural sciences is that people can not perform experiments to verify the the-oretical results. It is impossible and unethical to perform such experiments. Thus,mathematical models become important tools to analyze and compare differentpublic health strategies to plan for epidemics ([20]).There are several types of mathematical epidemic models. Deterministic epi-demic models divide the population under study into compartments and analyzethe time rate of transfer from one compartment to another. This type of models aregoverned by systems of difference equations, ordinary differential equations, delaydifferential equations or integral equations. The epidemic processes described aredeterministic, because the behaviour of a population is determined completely byits history and the rules which are described by models ([20]). If we start with asimple deterministic model, it is convenient to add control measures to this model,such as vaccination, treatment, or a combination of vaccination and treatment. Bycomparing the results of different control strategies, the optimal measure can beobtained. Stochastic epidemic models are studying the transmissions of epidemicsbased on the theory of probability. Several types of stochastic epidemic modelsare: discrete time Markov chain (both time and state variables are discrete), contin-uous time Markov chain (time is continuous, but the state variable is discrete) and1stochastic differential equations (based on a diffusion process). Some unique prop-erties to the stochastic epidemic models include the probability of an outbreak, thequasi-stationary probability distribution, the final size distribution of an epidemicand the expected duration of an epidemic (see Chapter 3 by Linda J.S. Allen in[20]). Network epidemic models (branching processes) are applied to the scenarioof outbreak of epidemics, when homogeneous mixing may not be a good assump-tion. At the beginning of a disease outbreak, there are few number of infectiveindividuals. The initial transmission of infection is a stochastic event, which de-pends on the contact patterns in the population. Once the epidemic has started, itis reasonable to switch to a deterministic epidemic model.We focus on the deterministic epidemic vaccination models in this thesis. Thethesis is divided into four parts.(A). In the first part, we introduce some basic definitions and theory of game the-ory, provide some background knowledge of deterministic epidemic modelsat the beginning. We construct and analyze four vaccination games: vaccina-tion game in homogeneous mixing population with perfect vaccine, vaccina-tion game in homogeneous mixing population with imperfect vaccine, vac-cination game in heterogeneous mixing population with perfect vaccine andvaccination game in heterogeneous mixing population with imperfect vac-cine. For each game, coupling with the replicator equations, we are able todescribe the evolutionary dynamics of the game and obtain the correspond-ing Nash equilibrium. Then we use the deterministic epidemic model togeneralize the final size relation, and further find the attack ratio being animplicit function of the vaccine coverage level. The combination of thesetwo steps of work provides us the better estimation of the expected vaccinecoverage level. Following the analysis of vaccination games, we start the dis-cussions of herd immunity. The herd immunity threshold in homogeneousmixing population has been studied very well; we re-define the herd im-munity in heterogeneous mixing population. The herd immunity thresholdrelations are obtained, for case of perfect vaccine and case of imperfect vac-cine. In the last part of this chapter, we introduce some economic analysis indisease eradication, which is closely related to herd immunity.2(B). In the second part, we consider that, prior to the outbreak of epidemic,whether or not there is a unique Nash equilibrium strategy set for vacci-nation. The necessary condition for uniqueness is that the attack ratio is thedecreasing function of the vaccine coverage level. By using the deterministicepidemic models, we are able to express the attack ratios implicitly in termsof the vaccine coverage levels. For the case of homogeneous mixing popula-tion, we completely prove the monotonicity. For the case of heterogeneousmixing population, we are able to obtain the complete theoretical results forthe case of perfect vaccine. However, for the case of imperfect vaccine, wecan only prove the monotonicity partially. In the final part, we study the finalsize relation of age-of-infection models. By proving the similarity of formsof final size relations. We can easily extend the theoretical results to generalcompartmental epidemic models.(C). In the third part, we formulate and analyze the malaria vaccination model.We first introduce some background knowledge of vector borne diseases,such as: West Nile, dengue fever, zika and malaria. We then introducesome previous modelling work in vector borne diseases, the focus is alsoon malaria. Secondly, we formulate a delayed malaria transmission model,which includes both extrinsic incubation period and intrinsic incubation pe-riod. We calculate the basic reproduction number R0, the disease-free equi-librium and the possible endemic equilibrium (exists if and only if R0 islarger than the unity). Then we analyze the stability of equilibria of thismodel. We prove that ifR0 < 1, the disease-free equilibrium is stable; whileif R0 > 1, the endemic equilibrium is stable for any lengths of the delayperiods. This part of work is the modification of previous work on delayedmalaria transmission models. Several numerical simulations are providedto verify our theoretical results. Thirdly, we propose a malaria vaccinationmodel. We calculate the control reproduction number Rc and the expres-sion of the disease-free equilibrium. We are able to prove that the disease-free equilibrium is asymptotically stable if the control reproduction numberRc < 1. The threshold condition to prevent the outbreak of malaria is pro-vided as well, which is the so-called threshold of eradication.3(D). In the fourth part, based on research in previous chapters, we propose somepossible future work for further considerations, especially for modelling thedisease of malaria.4Chapter 2Game theory and vaccinationThis chapter consists of three main parts. In the first section, we introduce somebasic knowledge and theory in game theory and deterministic epidemic models.In the second section, we describe four types of vaccination games and analyzethese games. In the third section, we extend the discussion of vaccination to thecontext of ”disease eradication”, which is related to the herd immunity. In this sec-tion, we define the herd immunity threshold/relation in the heterogeneous mixingpopulation with two sub-groups.2.1 BackgroundGame theory has been applied to predict human behavior in the context of epi-demiology for over two decades. People choose the best strategy to maximize theirown payoffs, based on the outcomes of different strategies. Before the outbreakof the epidemic, every individual can choose vaccination, or undertake the risk ofinfection. Ultimately, a certain percentage of people goes into the vaccinated class,while others go to the susceptible class. Under these circumstances, the Nash equi-librium is reached; anyone adopting a different strategy will have a lower payoff.The basic method to predict the expected vaccine coverage level by game theoryis coupling game theory and deterministic epidemic models. From the analysisof game theory, an implicit form of the expected vaccine coverage level and theexpected attack ratio is obtained. Through the deterministic epidemic models, the5expected attack ratio can be expressed in terms of the vaccine coverage level. Fi-nally, an accurate estimation of vaccine coverage level can be made. In this section,we introduce some basic concepts and theory in game theory and deterministic epi-demic models, especially the vaccination models.2.1.1 Game theoryGame theory is the formal study of decision-making where several players mustmake choices to maximize their own payoffs. Game theory has been used in manydisciplines such as economics, political science, business, information science, andbiology as well. First paper in modern game theory was written by John von Neu-mann in 1928 in a theory of parlor games. The early work of game theory issummarized in the book ”The Theory of Games and Economic Behavior by Johnvon Neumann and Oskar Morgenstern (see [69]). In 1944, John Nash publishedhis thesis titled Non Cooperative games, to state the concept of equilibrium point(also well-known as Nash Equilibrium). In the 1950s and 1960s, game theory wasbroadly applied to the problems caused by politics and wars. Since 1970s, thegame theory has found its application in economic theory.Game theory has also been applied to the area of evolutionary biology, theideas were first introduced by John Maynard Smith and G. R. Price (see [52]). Theapplication of game theory in biology is called Evolutionary game theory. Since1980s, game theory was also been used to solve the problems in public health.We shall focus on reviewing previous work of game theory and epidemiology inchapter 3.2.1.1.1 Some important definitions in classic game theoryClassic game theory includes two different types of games: cooperative games andnon-cooperative games. For cooperative games, the central concept is the Paretooptimum, which has the definition:Definition 1. An outcome is Pareto optimal if any change of strategies would re-duce the payoff of at least one player.For non-cooperative games, the central concept is the well-known Nash equi-librium, which has the definition:6Definition 2. An outcome is Nash equilibrium if no player can improve its payof-f/utility by unilaterally changing its strategy.In the designed games, players always have multiple strategy to choose to max-imize their payoffs. Under certain circumstances, players may face the dominant ordominated strategies. In these cases, players need to choose the dominant strategyor avoid the dominated strategy.Definition 3. A pure strategy is dominant if it performs equally well or better thanany other strategies. This strategy is strictly dominant if it always performs betterthan any other strategies.Definition 4. A pure strategy is dominated if at least one other pure strategy per-forms equally well or better.In our later analysis in vaccination games, we will repeatedly use Nash equi-librium and show the existence of dominant strategy or dominated strategy.2.1.1.2 Definitions in evolutionary game theoryJohn Maynard Smith and George Price introduced the idea of evolutionary gametheory in 1970s (see [52]). The reproductive fitness is used to identify the utilityfunction (same as the payoff from the game). Successful strategies will spread inthe population. There are two aspects of evolutionary game theory: the first one isthat evolutionary game theory always deals with the populations instead of severalplayers; the second one is that it is more natural for evolutionary game theory tointroduce a dynamics into the system.Remark 1. Part of this section originally comes from the lecture notes of Math564: Evolutionary Dynamics, taught by prof. Christoph Hauert.Evolutionary game dynamics use the replicator equations to describe the evo-lution of the densities of the different strategies.x˙i = xi(pii−〈pi〉).In this replicator equation, xi is the fraction of type i in the population (differenttypes means different species or same specie with different strategies), pii is the7fitness (payoff) of this type and 〈pi〉 is the average fitness in the whole population.From this equation, if the fitness pii is smaller than the average fitness, the densityof type i increases; while if pii is bigger than the average fitness, the density of typei increases. In general, the value of pii depends on the fractions of all types in thepopulation.We consider a population with two types. And we use the following 2×2 gameto represent.ExampleType A Type BType A a bType B c dThis payoff matrix is used to describe interactions of two types: If A interactswith another A, the payoff is a; if A interacts with B, the payoff is b. Similar fortype B, if B interacts with type A, it gains c, and obtains d from interaction withanother B. We havex = x1 = 1− x2,where x1 is the fraction of type A and x2 is the fraction of type B. The payoffs arepiA = ax+b(1− x),piB = cx+d(1− x).The replicator equation isx˙ = x(1− x)[(a−b− c+d)x+b−d].There are three fixed points, two trivial ones are x = 0 and x = 1. The third one isx? =d−ba−b− c+d ,for a > c and d > b or for a < c and d < b.There are four generic cases (see [54] for more discussions):1. Dominance. If a > c and b > d, type A dominates type B. In this case, the8fixed point at x = 1 is stable and the fixed point at x = 0 is unstable. If a < cand b < d, type B dominates type A. In this case, the fixed point at x = 0 isstable and the fixed point at x = 1 is unstable.2. Bistability. If a > c and d > b, the fixed points at x = 0 and x = 1 are bothstable and the fixed point at x? is unstable.3. Coexistence. If a < c and d < b, the fixed points at x = 0 and x = 1 areboth unstable and the fixed point at x? is stable. The population eventuallybecomes a stable mixture of type A and type B. This is the case we will meetin the vaccination games.4. Neutrality. For a = c and b = d.2.1.2 Deterministic epidemic modelsThe content of this subsection is based on the lecture notes by prof. Fred Brauer, thebook mathmatical epidemiology ([20]) edited by Prof. Fred Brauer, prof. PaulineVan den driessche and prof. Jianhong Wu. Some research results are coming froma series of papers by prof. Fred Brauer.2.1.2.1 The simple KermackMcKendrick model and the basic reproductionnumberR0One of the early triumphs of mathematical epidemiology was the formulation ofthe Kermack-McKendrick model. W.O. Kermack and A.G. McKendrick describedthe transmission of communicable diseases with the basic compartmental modelsin a series of papers (see [42], [43], [44]). The general model discussed in thesepapers included dependence on age-of-infection. The special case of the Kerma-ckMcKendrick model is the well-known SIR compartmental model,S′ = −βSII′ = βSI−αIR′ = αI.9In this model, S(t) denotes the number of individuals who are susceptible to thedisease at time t. I(t) denotes the number of individuals who are infected andinfectious. R(t) denotes the number of individuals who have recovered from thedisease and are with full immunity against reinfection. The SIR model is based onthe following assumptions:(1). An average member of the population makes contact sufficient to transmitinfection with βN others per unit time, where N represents total populationsize.(2). There are no disease deaths, and the total population size is a constant N.(3). There is no entry into or departure from the population. In other words, weconsider the case of short time scale.(4). Infectives leave the infective class at rate αI per unit time. This indicatesthat the length of the infective period is distributed exponentially.The central conception of SIR model, or in mathematical epidemiology, is thebasic reproduction number R0. The basic reproduction number R0 is defined asthe expected number of secondary cases produced by the introduction of a singleinfection in the completely susceptible population. If the population is partiallyvaccinated or not completely susceptible, it is appropriate to use the terminologythe control reproduction number Rc. For the SIR model, the basic reproductionnumberR0 =βNα.R0 is used to determine whether the epidemic can invade the population. IfR0 > 1,there is an epidemic. If R0 < 1, the epidemic dies out. For more complicated de-terministic compartmental models, the way to calculateR0 may be more difficult.2.1.2.2 The final size relationFinal size relation is a relation between the basic reproduction number R0 and thesize of epidemic. It can be used to predict how serous will an epidemic can be. Weknow that S(t) is a decreasing non-negative function and has the limit S∞ > 0 as10t→ ∞. The sum of the first two equations in SIR model is(S+ I)′ =−αI.We know that S+ I is a non-negative decreasing function about time t, and tendsto a limit as t→ ∞. And,I∞ = limt→∞ I(t) = 0.Thus, S+ I has the limit S∞. We integrate the sum of the first two equation of SIRmodel from 0 to ∞ givesα∫ ∞0(S(t)+ I(t))dt = S0 + I0−S∞ = N−S∞.Division of the first equation of SIR model by S and integration from 0 to ∞ givesthe final size relationlnS0S∞= β∫ ∞0I(t)dt=βα[N−S∞]= R0[1− S∞N ].To prove the uniqueness of solution of the final size relation, we defineg(x) = lnS0x−R0[1− xN ].We haveg(0+)> 0, g(N)< 0,and g′(x)< 0 if and only if0 < x <NR0.Thus, ifR0 < 1, g(x) is monotone decreasing function with positive value at x= 0+and negative value at x = N. There is a unique solution of g(x) with x < N. IfR0 > 1, g(x) is monotone decreasing from a positive value at x= 0+ to a minimumat x = NR0 , then increase to a negative value at x = N. We have a unique solution11x < NR0 . In either case, the solution to the equation g(x) is unique.Final size relation is crucial in analyzing the deterministic compartmental mod-els. For models with disease deaths, we may obtain the final size inequality insteadof the final size relation.2.1.2.3 The next generation matrix approach and the basic reproductionnumberR0In some models with sub-populations with different susceptibility to infection, itis impossible to calculate the basic reproduction number R0 by following the sec-ondary cases caused by a single infective introduced into the population. The nextgeneration matrix approach is a more general way to calculate the basic reproduc-tion number R0 to the deterministic compartmental models (see [26], [27], [68]).A compartment is called a disease compartment if the individuals therein are in-fected. We assume that there are n disease compartments and m nondisease com-partments. We denote x ∈ Rn and y ∈ Rm be the sub-populations in each of thesecompartments. We also denote by F the rate secondary infections increase theith disease compartment and by Vi the rate disease progression, death and recov-ery decrease the ith compartment. The compartmental model can be written in thefollowing form:x′i = Fi(x,y)−Vi(x,y), i = 1, · · · ,n, (2.1)y′j = g j(x,y), j = 1, · · · ,m. (2.2)The derivation of the basic reproduction number R0 is based on the linearizationof the ODE model about a disease-free equilibrium. We make the following as-sumptions to ensure the existence of the equilibrium and to ensure the model iswell-posed.(1). Fi(0,y) = 0 and Vi(0,y) = 0 for all y≥ 0 and i = 1, · · · ,n.(2). Fi(x,y)≥ 0 for all nonnegative x and y and i = 1, · · · ,n.(3). Vi(x,y)≤ 0 whenever xi = 0 with i = 1, · · · ,n.(4). ∑ni=1Vi(x,y)≥ 0 for all nonnegative x and y.12(5). The disease-free system y′ = g(0,y) has a unique equilibrium that is asymp-totically stable.Assumption (1) ensures that the disease-free set, which consists of all pointsof the form (0,y), is invariant. This also ensures that the disease-free equilibriumis an equilibrium of the full system.Now a single case of infection is introduced into a susceptible population. Theinitial ability of the disease to spread through the population is determined by thelinearization of (2.1). Based on assumption (1), we have∂Fi∂y j(0,y0) =∂Vi∂y j(0,y0) = 0for every pair (i, j). So the linearized equations for the disease compartments, x,are decoupled from the remaining equations and can be written asx′ = (F−V )x, (2.3)where F and V are n×n matrix with entriesF =∂Fi∂x j(0,y0), V =∂Vi∂x j(0,y0).By assumption (5), linear stability of the compartmental model is completely de-termined by the linear stability of the matrix (F−V ) in (2.3).The number of secondary infections produced by a single infected individualcan be expressed as the product of the expected duration of the infectious periodand the rate secondary infections occur. For the general model with n diseasecompartments, these are computed for each compartment for a hypothetical indexcase. The expected time the index case spends in each compartment is given bythe integral∫ ∞0 φ(t,x0)dt, where φ(t,x0) is the solution to the equation (2.3) withF = 0 and nonnegative initial conditions x0:x′ =−V x, x(0) = x0. (2.4)13The solution shows the path of the index case through the disease compartmentsfrom the initial exposure through to death or recovery with the ith component ofφ(t,x0)interpreted as the probability that the index case, which is introduced attime t = 0, is in disease state i at time t. The solution of equation (2.4) isφ(t,x0) = e−Vtx0.Thus, ∫ ∞0φ(t,x0)dt =V−1x0,and the (i, j) entry of the matrix V−1 can be interpreted as the expected time anindividual initially introduced into disease compartment j spends in disease com-partment i. The (i, j) entry of the matrix F is the rate at which secondary infectionsare produced in compartment i by an index case in compartment j. The expectednumber of secondary infections produced by the index case is given by∫ ∞0Fe−Vtx0dt = FV−1x0, (2.5)the matrix K = FV−1 is referred to the next generation matrix for the system at thedisease-free equilibrium. The (i, j) entry of K is the expected number of secondaryinfections in compartment i produced by individuals initially in compartment j.The basic reproduction number R0 is defined as the largest eigenvalue of matrixK, which isR0 = ρ(FV−1).2.1.2.4 The deterministic epidemic models with heterogeneous mixingIn some disease transmission models, not all members of the population makecontacts at the same rate. There are often ”super-spreaders”, who make many con-tacts and are instrumental in spreading disease and in general some members ofthe population make more contacts than others. Based on different activity levels,we divide the whole population into two sub-groups (it is straightforward to ex-tend to more than two sub-groups). The sizes of two sub-populations are N1 andN2, respectively. Each sub-population is divided into susceptibles, infectives, and14removed classes with subscripts to identify the sub-population. Suppose that mem-bers of group i make ai contacts in unit time and that the fraction of contacts madeby a member of group i that is with a member of group j is pi j with (i, j = 1,2).We havep11 + p12 = p21 + p22 = 1.We assume that the size of the whole population is constant. The two-group SIRmodel isS′1 = −[p11a1S1I1N1+ p12a1S1I2N2]I′1 = [p11a1S1I1N1+ p12a1S1I2N2]−α1I1S′2 = −[p21a2S2I1N1+ p22a2S2I2N2]I′2 = [p21a2S2I1N1+ p22a2S2I2N2]−α2I2with the mean infective period in group i is 1αi . It can be proved that (see [13],[16]),S1→ S1(∞)> 0, S2→ S2(∞)> 0,as t→ ∞.The basic reproduction numberR0 can be obtained through the next generationmatrix approach. The final size relations can also be obtained. The simple two-group SIR model can be easily extended to the models with more compartments,such as exposed class, vaccinated class. And this model can be extended to themodels with more than two sub groups. If we consider a population with a generaldistribution of sub-populations, see [19] for reference.In the two-group SIR model, the general mixing pattern is determined by thequantities p11 and p12. The total number of contacts made in unit time by membersof group 1 with members of group 2 is a1 p12N1 and this value must equal the totalnumber of contacts by members of group 2 with members of group 1. The balancerelation isp12a1N2=p21a2N1.15We consider two special cases of mixing pattern:(A). The proportionate mixing. The number of contacts between groups is pro-portional to the relative activity levels. In this case, we havepi j =a jN ja1N1 +a2N2,with i, j = 1,2.(B). The preferred mixing. A fraction pii of each group mixes randomly with itsown group and the remaining members mix proportionately. We havep11 = pi1 +(1−pi1)p1, p12 = (1−pi1)p2,andp21 = (1−pi2)p1, p22 = pi2 +(1−pi2)p2,withpii =(1−pii)aiNi(1−pi1)a1N1 +(1−pi2)a2N2 .Apparently, proportionate mixing is a special case of preferred mixing withpi1 = pi2 = 0.(C). The like-with-like mixing. Members of each group mixes only with mem-bers of the same group. In this case, we will have two separate SIR modelswith homogeneous mixing.2.1.2.5 The age-of-infection epidemic modelThe general epidemic model described by Kermack and McKendrick included adependence of infectivity on the time since becoming infected. This allows infec-tive periods that are not exponentially distributed. The age-of-infection epidemicmodel has a more general form. In age-of-infection model, We denote S(t) thenumber of susceptibles at time t and denote φ(t) the total infectivity at time t. φ(t)is defined as the sum of products of the number of infected members with eachinfection age and the mean infectivity at that infection age. We also make the as-sumption that on the average members of the population make a constant number a16of contacts in unit time. We let B(τ) be the fraction of infected members remaininginfected at infection age τ and let pi(τ) with 0≤ pi(τ)≤ 1 be the mean infectivityat infection age τ . We haveA(τ) = pi(τ)B(τ),the mean infectivity of members of the population with infection age τ . We assumethe total population size is a constant N.The age-of-infection epidemic model isS′(t) = −βSφφ(t) = φ0(t)+∫ t0βS(t− τ)φ(t− τ)A(τ)dτ= φ0(t)+∫ t0[−S′(t− τ)]A(τ)dτ.φ0(t) is the total infectivity of the initial infectives at time t.The basic reproduction number of the age-of-infection model isR0 = βN∫ ∞0A(τ)dτ.To find the final size relation of the age-of-infection model, we have−S′(t)S(t)= βφ0(t)+β∫ t0[−S′(t− τ)]A(τ)dτ.We integrate this equation with respect to t from 0 to ∞ to havelnS0S∞= β∫ ∞0φ0(t)dt +β∫ ∞0∫ t0[−S′(t− τ)]A(τ)dτdt= β∫ ∞0φ0(t)dt +β∫ ∞0A(τ)∫ ∞τ[−S′(t− τ)]dtdτ= β∫ ∞0φ0(t)dt +[S0−S∞]∫ ∞0A(τ)dτ= β [N−S∞]∫ ∞0A(τ)dτ+β∫ ∞0[φ0(t)− (N−S0)A(τ)]dτ= R0[1− S∞N ]−β∫ ∞0[(N−S∞)A(t)−φ0(t)]dt.17If all initial infectives have infection age zero at t = 0, we have φ0(t)= [N−S0]A(t).Further, we have ∫ ∞0[φ0(t)− [N−S0]A(t)]dt = 0.The general final size relation has formula,lnS0S∞=R0(1− S∞N ).If there are initial infectives with infection age greater than zero, we obtain afinal size inequality. If the assumed initial term is small and can be omitted, thefinal size relation still holds.2.2 Vaccination gamesVaccination benefits public health enormously. When a large percentage of a pop-ulation is vaccinated and has immunity to the disease, the part of the populationwho is not vaccinated benefits from the indirect protection. People choose the beststrategy to maximize their own payoffs, based on the outcomes of different strate-gies. Before the outbreak of the epidemic, every individual can choose vaccination,or undertake the risk of infection. Ultimately, a certain percentage of people goesinto the vaccinated class, while others go to the susceptible class. Under thesecircumstances, the Nash equilibrium is reached, anyone adopting a different strat-egy will have a lower payoff. Previous research has indicated that the tendency ofindividuals to optimize self-interest (maximize the payoff) can lead to vaccinationlevel, and this level is suboptimal for the whole community (see [57]). By couplinggame theory and epidemic compartmental models, the expected vaccine coveragelevel can be estimated. This level is lower than the Herd Immunity threshold. Wedemonstrate 4 vaccination games to start the investigate the vaccination programs.We first define the attack ratio in population, this definition will be used throughthe thesis repeatedly.Definition 5. The attack ratio in the population is defined as:pi =Ni−N fNi= 1− N fNi,18where Ni is the size of the whole population, N f is the number of individuals whoare not affected by the epidemic.2.2.1 Vaccination games in homogeneous mixing populationWe first analyze the vaccination games in a homogeneous mixing population. Forevery individuals in the population, there are only two options to choose: to bevaccinated or not to be vaccinated. The scenario we discuss here is the pre-outbreakof epidemic. People should make their decisions prior to the outbreak. All gameswe discuss here are based on self-interest, not the community-interest.2.2.1.1 Vaccination game with perfect vaccineVaccination can offer the protection against the disease. But vaccination also bringsthe cost, which can be the combination of loss of time, monetary cost for beingvaccinated and other side effects of the vaccine. From one individual’s perspec-tive, how to choose the proper strategy to maximize the payoff depends on twofactors: the strategy s/he chooses and other individuals’ choices. We formalize thisvaccination game by setting the elements of the following payoff matrix:Fully-effective Vaccination GameVaccinated individual j Susceptible individual jVaccinated individual i −Cv −CvSusceptible individual i 0 −pipCiThis payoff matrix describes four outcomes of the game: individual i is vacci-nated/not vaccinated and interacts with individual j who is vaccinated/not vacci-nated. The elements in the payoff matrix are the payoffs to individual i. If individ-ual i is vaccinated, the payoff is cost of the vaccination −Cv. This is true whetherthis individual is matched with a vaccinated or not vaccinated individual, becausethe vaccine is fully effective and there is no chance for individual i to get infected.If individual i chooses not to be vaccinated, there are two possibilities. The firstpossibility is to match with a vaccinated individual. In this case, since individualj is protected by the vaccine and will not be infected, individual i will not be in-fected, so the payoff is 0. The second possibility is to match with a not vaccinated19individual. We assume for everyone who is susceptible, the probability of beinginfected is pip, and the cost of infection is −Ci. We also assume that if individual jis infected, individual i will be infected through the contact (to simplify the analy-sis, we assume that each contact generates one more infection case; in most cases,the probability is less than 1, but the analysis is similar and straightforward). Thecosts of vaccination and infection were indicated as negative terms, because theydecrease the payoffs. Similar games can be found in [65]. It is noted that if the costof being vaccinated (Cv) is larger than the cost of being infected (Ci), the strategyof not being vaccinated is the dominant strategy. Since to be vaccinated is morerisky is than being infected.The second step is to calculate the payoffs of individual i choosing differentstrategies. Assume the final fraction of people who take vaccination is p; the frac-tion of people to be vaccinated is xv and the fraction of people not to be vaccinatedis xs, {xv,xs}= {p,1− p} is the equilibrium. We have the payoffs of two strategies:fv =−Cvxv−Cvxs =−Cv,fs =−pipCixs,(2.6)where fv is the payoff of people to be vaccinated and fs is the payoff of people notto be vaccinated. respectively.Then the replicator equations to this vaccination game can be expressed asx˙v = xv( fv− f¯ ),x˙s = xs( fs− f¯ ),(2.7)where f¯ denotes the average payofff¯ = xv fv + xs fs.By solving the replicator equations we find that the equilibrium is xv = 1− CvpipCi ,which implies that the final fraction p of population to be vaccinated is 1− CvpipCi . Ifwe define the relative cost r = CvCi (the relative cost reflects the relative risk of beingvaccinated and being infected), we can rewrite the expected vaccine coverage level20p = 1− rpip . Since pip is determined by the value of p, we rewrite the relation aspip(1− p) = r. (2.8)The relation (2.8) suggests:(1). If the vaccine coverage level p = 0, the left side of (2.8) is pi0, which is theattack ratio to the whole population if the whole population is susceptible.(2). If the vaccine coverage level p ∈ [pc,1], where pc indicates the herd immu-nity threshold value, the left side of (2.8) is 0.(3). Under the assumptions: pi0 > r and pip is a decreasing function about p,there is a unique solution to equation (2.8). We will prove that attack rate ismonotone decreasing with increased vaccination coverage in chapter 3.(4). The Nash equilibrium ”individuals choosing to be vaccinated with probabil-ity p” is convergently stable if it exists, see [30].To estimate the attack ratio pip, we need to use the deterministic epidemic SIRmodel,dSdt=−βSIdIdt= βSI−αIdRdt= αI,(2.9)where S(t) denotes the number of individuals who are susceptible to the disease,I(t) denotes the number of infected individuals and R(t) number of individuals whohave been infected and then removed from the possibility of being infected again,with the total population size a constant N. p is the expected vaccine coverage, sothe initial conditions of this model are S0 = (1− p)N and I0 = 0.We first need to find the final size relation. Following the method we introducedearlier, the final size relation to the simple SIR model can be expressed as:lnS0S∞=βα[S0−S∞].21The ratio of attack ispip =S0−S∞S0,so we have the relation of pip and pln11−pip =βα(1− p)Npip. (2.10)Combining equations (2.8) and (2.10), the more accurate estimation of theexpected vaccine coverage level can be obtained. We observe the equation (2.8)that if individuals’ choice is purely motivated by self interest, the herd immunitycan not be reached. Since the expected vaccine coverage level p < pc and it issuboptimal.2.2.1.2 Vaccination game with imperfect vaccineIn real life, the vaccination may reduce but not completely eliminate susceptibilityto infection, we model this by introducing a factor σ , 0 < σ < 1, in the infectionrate of vaccinated members. σ = 0 means that the vaccine is perfectly effectiveand σ = 1 means that the vaccine has no effect. Other assumptions of the game arethe same as the case of fully-effective vaccination game.We set up the game by setting the entries of the payoff matrix. If individuali is vaccinated and individual j is also vaccinated, then the payoff for individual iis the cost of vaccination plus the risk of being infected, the cost of vaccination isstill −Cv; the probability that individual j getting infected is piv, since individual iis also vaccinated, the probability of being infected by individual j is reduced byσ , and the cost of infection is −Ci. The other 3 entries in the payoff matrix can begenerated in the same way. It is noted that pip is the probability that individual isnot vaccinated but infected. The formal partial-effective vaccination game is:Partial-effective Vaccination GameVaccinated individual Susceptible individualVaccinated individual −Cv−σpivCi −Cv−σpipCiSusceptible individual −pivCi −pipCi22Remark 2. (A). If the relative cost r > (1−σ)max{piv,pip} for any value of p,the strategy of not to be vaccinated is the dominant strategy in this game.(B). It is impossible that strategy of to be vaccinated be the dominant strategy inthis game.The replicator equations to this vaccination game can be expressed asx˙v = xv( fv− f¯ ),x˙s = xs( fs− f¯ ),(2.11)where f¯ denotes the average payoff, payoffs fv and fs are defined as:fv = (−Cv−σpivCi)xv +(−Cv−σpipCi)xs,fs =−pivCixv−pipCixs.(2.12)Solving the replicator equations, we have the expression of fraction pr1−σ = ppiv +(1− p)pip, (2.13)where r represents the relative cost CvCi . By similar analysis, one necessary conditionfor the uniqueness of the Nash equilibrium is that pip and piv are decreasing functionabout the vaccine coverage level p, we will prove this in chapter 3. An equationwhich is similar to (2.13) can be found in [33] by Fine and Clarkson.We now use the deterministic epidemic SVIR model to make more accurateestimations. The SVIR model isdSdt=−βSI,dVdt=−σβV I,dIdt= βSI+σβV I−αI,dRdt= αI,(2.14)with initial values S0 = (1− p)N, V0 = pN, I0 = 0, other parameters are defined asbefore.23Figure 2.1: SVIR modelThe control reproduction numberRc has the expression,Rc =βα(S0 +σV0),and the final size relation can be expressed by the following equations,lnS∞S0=βα[S∞+V∞−S0−V0],lnV∞V0=σβα[S∞+V∞−S0−V0].(2.15)We are able to obtain the formulas of the attack ratios to susceptible class andvaccinated classpip =S0−S∞S0,piv =V0−V∞V0.We go back to equation (2.13), by combing the game theory and the deter-ministic epidemic model, the accurate expected vaccine coverage level can be es-timated.242.2.2 Vaccination games in heterogeneous mixing populationIn this part, we discuss two vaccination games in heterogeneous mixing population,with perfect vaccine and with imperfect vaccine. The assumption in previous partthat the population is perfectly well-mixed is unrealistically simple, due to differentactivity levels, such as: different contact rates, different costs of infection and vac-cination. Each sub-population is divided into susceptibles, vaccinated, infectivesand removed members with subscripts to identify the subpopulation. Mixing ineach subpopulation is assumed to be homogeneous. We first make some assump-tions, which can be found in [16]:(1). The sizes of two sub-populations are N1 and N2, respectively. N1 and N2 areconstants.(2). Group i members make ai contacts in unit time and the fraction of contactsmade by a member of group i is with a member of group j is pi j ,(i, j = 1,2).2.2.2.1 Vaccination game with perfect vaccineWe first consider the simpler case that the vaccine is perfectly effective. To con-struct the game, we assume costs of vaccination for members in sub-group 1 isCv1 and for members in sub-group 2 is Cv2; probabilities of being infected for twogroups are pip1 and piq2 associated with the vaccine coverage level p and q in eachsubgroup, respectively. The formal vaccination game with two-subgroups is:Fully-effective vaccination game with two-subgroupsVaccinatedindividual 1Susceptibleindividual 1Vaccinatedindividual 2Susceptibleindividual 2Vaccinatedindividual 1−Cv1 −Cv1 −Cv1 −Cv1Susceptibleindividual 10 −pip1Ci1 0 −piq2Ci1Vaccinatedindividual 2−Cv2 −Cv2 −Cv2 −Cv2Susceptibleindividual 20 −pip1Ci2 0 −piq2Ci225The explanation of entries of this payoff matrix can be translated easily. Sincethe vaccine is perfect effective for both sub-groups, once individuals are beingvaccinated, the payoff only comes from the cost of vaccination. We defineCr1 =Cv1Ci1,andCr2 =Cv2Ci2,to represent the relative risks for subgroup 1 and subgroup 2, respectively.The final fraction of individuals to be vaccinated in group 1 is p; and the finalfraction of individuals to be vaccinated in group 2 is q; in group i, the fraction of in-dividuals to be vaccinated is xvi and the fraction of individuals not to be vaccinatedis xsi, with i = 1,2. The replicator equations are˙xv1 = xv1( fv1− f¯1),˙xs1 = xs1( fs1− f¯1),˙xv2 = xv2( fv2− f¯2),˙xs2 = xs2( fs2− f¯2),with f¯i denotes the average fitness of group i. Now we need to calculate the payoffs,apparently fv1 and fv2 do not depend on the fractions of individuals to be vaccinatedin two groups, andfv1 =−Cv1,fv2 =−Cv2.fs1 and fs2 depend on the fractions, andfs1 =−pip1Ci1α11−piq2Ci1α12,fs2 =−pip1Ci2α21−piq2Ci2α22,26where αi j is the probability of individual i who is non-vaccinated in group i, tocontact with individual j who is also non-vaccinated in group j,α11 = β11(1− p), α12 = β12(1−q), α21 = β21(1− p), α22 = β22(1−q),where βi j is the fraction of contacts made by a member of group i that is with amember of group j. So thatfs1 =−pip1Ci1β11(1− p)−piq2Ci1β12(1−q),fs2 =−pip1Ci2β21(1− p)−piq2Ci2β22(1−q).Solving the replicator equations, we obtain the equilibriump = 1− Cr1β22−Cr2β21β11β22−β12β211pip1, q = 1− Cr2β11−Cr1β12β11β22−β12β211piq2, (2.16)where Cr1 and Cr2 are relative costs in two groups.Remark 3. (A). Vaccine coverage levels p and q are determined by the mixingpattern, the nature of the disease. But vaccine coverages in the differentgroups are independent.(B). Behavior in group 1 dose not affect behavior in group 2, and vice versa.To find the attack ratios pip1 and piq2, we need to use the two-group SIR model,and the final size relations. The two-group SIR model is:S′1 =−[a1β11S1I1N1+a1β12S1I2N2],I′1 = [a1β11S1I1N1+a1β12S1I2N2]−α1I1,R′1 = α1I1,S′2 =−[a2β21S2I1N1+a2β22S2I2N2],I′2 = [a2β21S2I1N1+a2β22S2I2N2]−α2I2,R′2 = α2I2.27The final size relations can be revealed from the model,lnS1(0)S1(∞)=a1β11N1(1− p)N1−S1(∞)α1+a1β12N2(1−q)N2−S2(∞)α2,andlnS2(0)S2(∞)=a2β21N1(1− p)N1−S1(∞)α1+a2β22N2(1−q)N2−S2(∞)α2.These two equations are implicit equations about pip1 and p, and implicit equationabout piq2 and q. Based on the definitions of attack ratios,pip1 =S1(0)−S1(∞)S1(0),piq2 =S2(0)−S2(∞)S2(0),and combining with the equation in (2.16), the expected vaccine coveragelevels p and q can be obtained.2.2.2.2 Vaccination game with imperfect vaccineWe now consider the vaccination game with imperfect vaccine. The vaccine ispartial-effective to both sub-groups. We introduce the infection rates of vaccinatedmembers σ1 and σ2 in the game, with 0≤ σ1 ≤ 1 and 0≤ σ2 ≤ 1 (for group 1 andgroup 2, respectively). We construct the game by the following payoff matrix,28Partial-effective Vaccination Game with two-subgroupsVaccinatedindividual 1Susceptibleindividual 1Vaccinatedindividual 2Susceptibleindividual 2Vaccinatedindividual 1−Cv1 −σ1pipv1Ci1−Cv1 −σ1pipv1Ci1−Cv1 −σ1piqv2Ci1−Cv1 −σ1piqi2Ci1Susceptibleindividual 1−pipv1Ci1 −pipi1Ci1 −piqv2Ci1 −piqi2Ci1Vaccinatedindividual 2−Cv2 −σ2pipv1Ci2−Cv2 −σ2pipi1Ci2−Cv2 −σ2piqv2Ci2−Cv2 −σ2piqi2Ci2Susceptibleindividual 2−pipv1Ci2 −pipi1Ci2 −piqv2Ci2 −piqi2Ci2Here are some explanations of the payoff matrix:(1). pipv1 is the probability that individual 2 is infected, individual 2 is vaccinatedis the member of group 1. The vaccine coverage level in group 1 is p. Ifindividual 1 is vaccinated, then the probability that s/he being infected isreduced by σi, which depends on which group s/he is in. Similarly, pipi1,piqv2 and piqi2 can be translated in this way.(2). Other notations have same meanings as we defined in the previous sections.To find the final fractions p and q in both sub-groups, we need to find the NashEquilibrium of this game. The first step is to construct the replicator equations,˙xv1 = xv1( fv1− f¯1),˙xs1 = xs1( fs1− f¯1),˙xv2 = xv2( fv2− f¯2),˙xs2 = xs2( fs2− f¯2),29with f¯i denotes the average fitness of group i. The fitnesses/payoffs of each groupcan be expressed asfv1 =−[Cv1 +σ1pipv1Ci1]β11 p− [Cv1 +σ1pipi1Ci1]β11(1− p)− [Cv1 +σ1piqv2Ci1]β12q− [Cv1 +σ1piqi2Ci1]β12(1−q),fs1 =−pipv1Ci1β11 p−pipi1Ci1β11(1− p)−piqv2Ci1β12q−piqi2Ci1β12(1−q),fv2 =−[Cv2 +σ2pipv1Ci2]β21 p− [Cv2 +σ2pipi1Ci2]β21(1− p)− [Cv2 +σ2piqv2Ci2]β22q− [Cv2 +σ2piqi2Ci2]β22(1−q),fs2 =−pipv1Ci2β21 p−pipi1Ci2β21(1− p)−piqv2Ci2β22q−piqi2Ci2β22(1−q),where βi j (i, j = 1,2) reflects the mixing pattern.Solving the replicator equation, we obtain the final fractionsp =Cr1(β11β22+β12β22)1−σ1 −Cr2(β12β21+β12β22)1−σ2 −pipi1(β11β22−β12β21)(pipv1−pipi1)(β11β22−β12β21) , (2.17)q =Cr1(β11β21+β12β21)1−σ1 −Cr2(β11β21+β11β22)1−σ2 −piqi2(β12β21−β11β22)(piqv2−piqi2)(β12β21−β11β22) , (2.18)where Cr1 and Cr2 are relative costs in two groups. The attack ratios pipi1, piqi2,pipv1 and piqv2 are implicit functions about vaccine coverage levels. Thus, to makea better prediction, we need to use the deterministic model to express the attack30ratios. The two-group SVIR model is:S′1 = −[a1β11S1I1N1+a1β12S1I2N2],V ′1 = −σ1[a1β11V1I1N1+a1β12V1I2N2],I′1 = [a1β11S1I1N1+a1β12S1I2N2]+σ1[a1β11V1I1N1+a1β12V1I2N2]−α1I1,R′1 = α1I1,S′2 = −[a2β21S2I1N1+a2β22S2I2N2],V ′2 = −σ2[a2β21V2I1N1+a2β22V2I2N2],I′2 = [a2β21S2I1N1+a2β22S2I2N2]+σ2[a2β21V2I1N1+a2β22V2I2N2]−α2I2,R′2 = α2I2.The final size relations to the SVIR compartmental model arelnS1(0)S1(∞)=a1β11α1N1[N1−S1(∞)−V1(∞)]+ a1β12α2N2 [N2−S2(∞)−V2(∞)],lnV1(0)V1(∞)= σ1[a1β11α1N1[N1−S1(∞)−V1(∞)]+ a1β12α2N2 [N2−S2(∞)−V2(∞)]],lnS2(0)S2(∞)=a2β21α1N1[N1−S1(∞)−V1(∞)]+ a2β22α2N2 [N2−S2(∞)−V2(∞)],lnV2(0)V2(∞)= σ2[a2β21α1N1[N1−S1(∞)−V1(∞)]+ a2β22α2N2 [N2−S2(∞)−V2(∞)]].The attack ratios are defined aspipi1 =S1(0)−S1(∞)S1(0),piqi2 =S2(0)−S2(∞)S2(0),pipv1 =V1(0)−V1(∞)V1(0),piqv2 =V2(0)−V2(∞)V2(0).31The analysis provides us more accurate estimations of the expected vaccinecoverage levels p and q in each sub-population, by coupling the game theory andthe deterministic epidemic model.2.3 Disease eradication2.3.1 Herd immunity: introduction and mathematical formulationHerd immunity is the indirect protection against an infectious disease, which de-rives from the mass vaccination program. When a large percentage of a populationis vaccinated and has immunity to the disease, the part of the population who is notvaccinated benefits from the indirect protection. To further eradicate an infectiousdisease in the population, it is certainly unnecessary to vaccinate every individ-ual. As the percentage of the vaccinated individuals reaching to the herd immunitythreshold, there is little opportunity for disease to spread since there are not manysusceptible individuals, the chains of infection are blocked and the transmissionof disease is prevented. Mathematically speaking, the control reproduction num-ber Rc is decreased below 1 if the vaccine coverage level is higher than the herdimmunity threshold.Mathematical estimation of the herd immunity threshold is obtained from thesimple SIR compartmental model, and the expression is in terms of the basic re-production number R0,pc = 1− 1R0.This equation is constructed based on two assumptions:(1). the population is homogeneous, or well mixed. All individuals have sameactivity level,(2). the vaccine is perfectly effective, or at least highly effective. If the effec-tiveness of the vaccine is relatively low, the expression of the herd immunitythreshold should be also in terms of the vaccine efficacy.The influential paper ”Herd immunity: basic concept and relevance to public healthimmunization practices” [34] was published in 1971 by Fox et al.. This paper intro-duced important mathematical methods to study the herd immunity and obtained32conclusions. Since then, a lot of research work explored this interesting subjectfrom the standpoint of mathematics.The history of herd immunity goes back to nearly a century ago. The term”herd immunity” first appeared in [66] by Topley and Wilson in 1923. AlthoughTopley and Wilson did not distinguish the difference between direct protection andindirect protection which are deriving from vaccination, this term is widely used byepidemiologists since then. The presence of indirect protection from vaccine was inthe 19th century, in [31], Farr noted that ”The smallpox would be distributed, andsometimes arrested, by vaccination, which protected a part of the population...”.Later the ”mass action principle” provides more logical explanation for indirectprotection by herd immunity. During the last century, herd immunity is mainlydiscussed in the context of disease eradication programs, the successes of the globalsmallpox eradication program and polio eradication program are perfect examples.More details on the history of study of herd immunity can be found in [32].The first mathematical model to study the herd immunity was constructed byHamer in 1906, which can be found in [39]. The model he used is discrete, thetime period is defined as the average interval between successive cases in a chainof transmission. Hamer argued that the number of transmissions per measles caseswas a function of the number of susceptibles in the population (see [32]), which isexpressed as:Ct+1 =CtStr, (2.19)where Ct is the number of infectives in time period t, St is number of susceptiblesin time period t, and r is the transmission parameter. To simulate the successivechange over time, Hamer argued that:St+1 +Ct+1 = St +Bt , (2.20)where Bt is the number of susceptibles added in time period t. From time period tto time period t+1, there is total number of St +Bt to become either susceptible orinfective, all cases of infection are recovered from the disease with immunity. Thismodel describes the dynamics of disease spreading, also defines the threshold that33the disease can spread. In (2.19), the disease won’t die out if and only ifCt+1 >Ct , (2.21)which implies thatSt >1r. (2.22)Following the inequalities, we assume the population size is N, then we define theherd immunity threshold as:H = 1− 1rN. (2.23)Before the outbreak of disease, if the vaccination decreases the number of suscep-tibles below the herd immunity threshold, the cases of infection should decline,ultimately die out. Later in the paper [48] by George Macdonald, he formulated”basic case reproduction rate” which is used to determine whether the infection canpersist in the population. By combining (2.23) and the basic case reproduction rateR0, we have the expression for the herd immunity threshold:H = 1− 1R0.Another important theoretical contribution is the Reed Frost model by Lowell Reedand Wade Frost in the 1920s. This model discussed the herd immunity in a closedpopulation.With the developments in differential equations and mathematical epidemiol-ogy, the theory of herd immunity has been better elaborated. More literatures canbe found in [2]. We first assume that the vaccine is perfect. The model we use tofind the herd immunity threshold is the simple SIR compartmental model:S′ =−βSI,I′ = βSI−αI,R′ = αI.(2.24)The population being studied is divided into compartments, namely a susceptibleclass S, an infective class I, and a removed class R, which can be eliminated from34the equations if there is no demographic change to the population. More detailsabout this model can be found in [13] and [14].In this simple SIR model,R0 has the expression,R0 =βNα. (2.25)The condition that the disease can spread among the population is R0 > 1, oth-erwise the disease will die out. If the population is partially vaccinated, and weassume the vaccine coverage level is p, then R0 is decreased to (1− p)R0. Sincethe population is not wholly susceptible, we use the control reproduction numberRc instead of R0. If we could make Rc < 1, such that the disease will die out. ItrequiresRc < 1,(1− p)R0 < 1,p > 1− 1R0.(2.26)When the vaccine coverage level p is higher than 1− 1R0 , the disease will die out.So we havepc = 1− 1R0, (2.27)which is the herd immunity threshold.Remark 4. It is noted that in (2.27), the herd immunity threshold only dependson the basic reproduction number R0. For less contagious infectious disease, therequired herd immunity threshold is relatively low; however, for very contagiousinfectious disease, the required herd immunity threshold is relatively high. Forexample: for Ebola, the estimated basic reproduction numberR0 is between 1.5 to2.5, the herd immunity threshold is estimated between 33% to 60%; for Influenza,R0 is between 1.5 to 1.8, the required herd immunity threshold is between 33%to 44%; but for measles, R0 is between 12 to 15, such that the required herdimmunity threshold is between 92% to 95%. A more complete list regarding basicreproduction numbers R0 and the herd immunity thresholds can be found in 2.1(according to CDC and WHO).In real life, the vaccine is normally imperfect. We model this by including a35Table 2.1: R0 and herd immunity thresholdInfectious diseases R0 Herd immunity thresholdMeasles 12-18 92%-95%Pertussis 12-17 92%-94%Diphtheria 6-7 83%-86%Rubella 6-7 83%-86%Smallpox 5-7 80%-86%Polio 5-7 80%-86%Mumps 4-7 75%-86%SARS 2-5 50%-80%Ebola 1.5-2.5 33%-60%Influenza 1.5-1.8 33%-44%factor σ , with 0 < σ < 1, in the infection rate of vaccinated members. With σ = 0meaning that the vaccine is perfectly effective and σ = 1 meaning that the vaccinehas no effect (see [13] and [14]). We use the SVIR compartmental model to obtainthe herd immunity threshold,dSdt=−βSI,dVdt=−σβV I,dIdt= βSI+σβV I−αI,dRdt= αI,(2.28)where V is the vaccinated class. Other assumptions are the same as the case ofperfect vaccine.If the population is partially vaccinated with the vaccine coverage level p 6= 0,36the control reproduction number isRc =βα(S0 +σV0)=βαN((1− p)+σ p)=R0(1− p+σ p)=R0(1− (1−σ)p).(2.29)IfRc < 1, then we haveR0(1− (1−σ)p)< 1,p >1− 1R01−σ .(2.30)We define the herd immunity threshold with the imperfect vaccine bepc =1− 1R01−σ . (2.31)Remark 5. If the vaccine is imperfect, the herd immunity threshold (2.31) is deter-mined by the basic reproduction numberR0 and the vaccine efficacy σ . If the vac-cine is very low efficient, with σ close to 1, the herd immunity can not be reached.2.3.2 Herd immunity in heterogeneous mixing populationIt is unrealistic to assume that the population is homogeneous. If we consider aheterogeneous mixing population with two subgroups, how to reach the herd im-munity threshold in such non homogeneous mixing population? It is obvious thatin this case, the herd immunity threshold is not a threshold value, but a thresholdrelation that the vaccine coverage levels in every subgroups have to satisfy. If pos-sible, can the herd immunity be achieved by vaccinating only one subgroup? Theresults may be applied to two scenarios.(1). We consider the population being divided into two subgroups. The first sub-group includes children and elderly people. The second subgroup includesother people of this population. The first subgroup may be vulnerable to the37side effects of the vaccine, they are too sick to be vaccinated. In such case,can we only vaccinate the second subgroup to reach the herd immunity? Ifso, the first group will obtain better protection, they do not need to undertakethe risks of the vaccination and infection.(2). We consider in the whole population, there is a certain percentage of thevaccine skeptics. If the population is divided into two subgroups: vaccineskeptics and vaccine believers. The only option to reach the herd immunityis to vaccinate the subgroup of the vaccine believers. We are interested inthe vaccine coverage level in the subgroup of the vaccine believers to reachthe herd immunity.Motivated by such questions, it is valuable to make more accurate estimation of theherd immunity threshold in the heterogeneous mixing population. We still considertwo separate scenarios: perfect vaccine and imperfect vaccine. For each case,we use the compartmental deterministic model to obtain the basic reproductionnumber R0 and the control reproduction number Rc. Then we construct the herdimmunity threshold relation based on these expressions.2.3.2.1 Perfect vaccineThe epidemic models with heterogeneous mixing has been very well studied in[18] and [19]. We first make the following assumptions:(1). Sizes of population in two subgroups are N1 and N2, N1 and N2 are constants,(2). group i members make ai contacts in unit time and that the fraction of con-tacts made by a member of group i that is with a member of group j is pi j,with i, j = 1,2 (see [18]).38The two-group SIR vaccination model isS′1 =−[a1 p11S1I1N1+a1 p12S1I2N2],I′1 = [a1 p11S1I1N1+a1 p12S1I2N2]−α1I1,R′1 = α1I1,S′2 =−[a2 p21S2I1N1+a2 p22S2I2N2],I′2 = [a2 p21S2I1N1+a2 p22S2I2N2]−α2I2,R′2 = α2I2,(2.32)with(1). Mean infective period in group i is 1αi , with i = 1,2,(2). vaccine coverage levels in each subgroup are v1 and v2.We use the next generation matrix approach, which is described in [68], tofind the basic reproduction number R0 and the control reproduction number Rc.F is the vector of rates of appearance of new infections in disease compartmentsand V is the vector of rates of transfer of individuals into/out of the non-diseasecompartments by all other means,F =[a1 p11 S1I1N1 +a1 p12S1I2N2a2 p21 S2I1N1 +a2 p22S2I2N2],andV =[α1I1α2I2].The 2×2 non-negative matrix F and a 2×2 non-singular M-matrix V areF =[a1 p11, a1 p12 N1N2a2 p21 N2N1 , a2 p22],39andV =[α1, 00, α2].We calculate the basic reproduction number R0 as the largest eigenvalue of thematrix FV−1,FV−1 =[a1 p11α1a1 p12N1α2N2a2 p21N2α1N1a2 p22α2].The basic reproduction numberR0 can be expressed as:R0 =a1 p11α1 +a2 p22α2 +√a21 p211α21+a22 p222α22− 2a1a2 p11 p22α1α2 +4a1a2 p12 p21α1α22. (2.33)By using the same method, if v1 > 0 or v2 > 0, we find the control reproductionnumberRc,Rc =a1 p11α1 (1− v1)+a2 p22α2 (1− v2)+(1− v1)(1− v2)√a21 p211α21+a22 p222α22− 2a1a2 p11 p22α1α2 +4a1a2 p12 p21α1α22.(2.34)The expression (2.34) is too complicated and less useful. We want to writeRcin terms of v1, v2 andR0,Rc = (1− v1)(1− v2)R0 + v2(1− v1)a1 p112α1 + v1(1− v2)a2 p222α2. (2.35)To reach the herd immunity, it requires(1− v1)(1− v2)R0 + v2(1− v1)a1 p112α1 + v1(1− v2)a2 p222α2< 1, (2.36)(2.36) is the herd immunity threshold relation which v1 and v2 have to satisfy. Itis noted that in (2.36), the pattern of mixing in the population will affect the herdimmunity.We discuss the herd immunity threshold relation (2.36) in three aspects.(A). Assume v1 = 0, which means the vaccine coverage level in subgroup 1 is 0.40Can we only vaccinate subgroup 2 to reach the herd immunity? We have(1− v2)R0 + v2 a1 p112α1 < 1,(R0− a1 p112α1 )v2 >R0−1.Because the basic reproduction number R0 > 1, and from the expression ofR0, we haveR0 >a1 p112α1 . We writev2 >R0−1R0− a1 p112α1. (2.37)If the subgroup 1 is completely susceptible, the vaccine coverage level insubgroup 2 has to be higher than R0−1R0− a1 p112α1. It is impossible to determinewhether R0−1R0− a1 p112α1is bounded by 0 and 1, so the threshold may not exist. Therelation between a1 p112α1 and 1 plays an important role.(B.) Assume v1 = 1, every individual in subgroup 1 is vaccinated, herd immunitycan definitely be achieved under this circumstance. We now have(1− v2)a2 p222α2 < 1,sov2 > 1− 2α2a2 p22 . (2.38)We observe that the herd immunity must be achieved in this case. Evenmore, if 1 < 2α2a2 p22 , we only need to vaccinate partial members in subgroup1, since the right hand side of this inequality is negative. When the vaccinecoverage level in subgroup 1 is high enough, the herd immunity for the wholepopulation can be reached.(C). Based on the above discussion, the relations between ai pii2αi and 1 are crucial,with i = 1,2. If a1 p112α1 < 1, without vaccinating anyone in subgroup 1, onlythe level of vaccine coverage in subgroup 2 reaches a certain value, the herdimmunity can be achieved; else if a1 p112α1 > 1, even every individuals in sub-41groups 2 are vaccinated, to reach the herd immunity, part of subgroup 1 hasto be vaccinated. We conclude that ai pii2αi is a population mixing thresholdvalue, conditions will be varied due to the relation with 1.We conclude the analysis with the following theorem,Theorem 2.3.1. In the two-subgroups heterogeneous mixing population with per-fect vaccine, if the vaccine coverage levels in two subgroups satisfy the inequality(2.36), the herd immunity can be reached in the whole population. Ifa1 p112α1< 1,when the vaccine coverage level in subgroup 2 reaches a certain percentage, with-out vaccinating anyone in subgroup 1, the herd immunity can be reached in thewhole population; while ifa1 p112α1> 1, to reach the herd immunity, part of sub-group 1 has to be vaccinated.Remark 6. The population mixing threshold value ai pii2αi measures the number ofsecondary infections in subgroup i. If this value is smaller than 1, that means morenew infections occur in another subgroup. In this case, the priority is to vaccinatethe other subgroup.2.3.2.2 Imperfect vaccineWe model the vaccination model in heterogeneous mixing population with imper-fect vaccine by introducing factors σ1 and σ2, with 0 < σ1 < 1 and 0 < σ2 < 1,in the infection rates of vaccinated members. Other assumptions still hold. The42two-subgroups SVIR model isS′1 =−[a1 p11S1I1N1+a1 p12S1I2N2],V ′1 =−σ1[a1 p11V1I1N1+a1 p12V1I2N2],I′1 = [a1 p11S1I1N1+a1 p12S1I2N2]+σ1[a1 p11V1I1N1+a1 p12V1I2N2]−α1I1,R′1 = α1I1,S′2 =−[a2 p21S2I1N1+a2 p22S2I2N2],V ′2 =−σ2[a2 p21V2I1N1+a2 p22V2I2N2],I′2 = [a2 p21S2I1N1+a2 p22S2I2N2]+σ2[a2 p21V2I1N1+a2 p22V2I2N2]−α2I2,R′2 = α2I2.(2.39)If v1 = v2 = 0, the population is wholly susceptible. We find the basic reproductionnumberR0 by using the next generation matrix method,R0 =a1 p11α1 +a2 p22α2 +√a21 p211α21+a22 p222α22− 2a1a2 p11 p22α1α2 +4a1a2 p12 p21α1α22. (2.40)It is noticed that (2.33) and (2.40) have identical formula, this is because R0 isonly determined by the nature of the infectious disease and the population itself.If two subgroups are partially vaccinated, v1 > 0 or v1 > 0, the control repro-43duction numberRc has the expression,Rc =a1 p11α1 (1− v1 +σ1v1)+a2 p22α2 (1− v2 +σ2v2)2+(1− v1 +σ1v1)(1− v2 +σ2v2)√a21 p211α21+a22 p222α22− 2a1a2 p11 p22α1α2 +4a1a2 p12 p21α1α22,=a1 p11α1 (1− (1−σ1)v1)+a2 p22α2 (1− (1−σ2)v2)2+(1− (1−σ1)v1)(1− (1−σ2)v2)√a21 p211α21+a22 p222α22− 2a1a2 p11 p22α1α2 +4a1a2 p12 p21α1α22.(2.41)We express the control reproduction numberRc in (2.41) in terms of v1, v2 andR0,Rc = (1− (1−σ1)v1)(1− (1−σ2)v2)R0+(1−σ2)v2(1− (1−σ1)v1)a1 p112α1 +(1−σ1)v1(1− (1−σ2)v2)a2 p222α2.(2.42)To eradicate the disease, we need to reach the herd immunity, which requiresthe control reproduction numberRc being decreased below 1. We have(1− (1−σ1)v1)(1− (1−σ2)v2)R0+(1−σ2)v2(1− (1−σ1)v1)a1 p112α1 +(1−σ1)v1(1− (1−σ2)v2)a2 p222α2< 1,(2.43)(2.43) is the herd immunity threshold relation that v1 and v2 have to satisfy.We discuss this inequality in four aspects,(A). Assume v1 = 0, which means the vaccine coverage level in subgroup 1 is 0.Then(1− (1−σ2)v2)R0 +(1−σ2)v2 a1 p112α1 < 1,(a1 p112α1−R0)(1−σ2)v2 < 1−R0.SinceR0 > 1,44and from the expression of the basic reproduction numberR0,R0 >a1 p112α1.We havev2 >1−R0(a1 p112α1 −R0)(1−σ2). (2.44)Apparently, comparing with the case of perfect vaccine, the imperfect vac-cine needs higher vaccine coverage level to reach the herd immunity. Thevaccine coverage level value v2 is related to the efficiency factor σ2, but an-other factor σ1 in subgroup 1 has no impact on v2, that is because everyindividual in subgroup 1 is non-vaccinated.(B). Assume that v1 = 1, which means every individual in subgroup 1 is vacci-nated, thenσ1(1−(1−σ2)v2)R0+σ1(1−σ2)v2 a1 p112α1 +(1−σ1)(1−(1−σ2)v2)a2 p222α2< 1.So we have(1−σ2)[σ1(R0− a1 p112α1 −a2 p222α2)+a2 p222α2]v2 > σ1R0+(1−σ1)a2 p222α2 −1.Because 1−σ2 and σ1(R0− a1 p112α1 −a2 p222α2 )+a2 p222α2 are both positive, v2 mustbe larger than a certain value. For the right hand side of the above inequality,if we defineg(σ1) = σ1R0 +(1−σ1)a2 p222α2 −1,we haveg′(σ1) =R0− a2 p222α2 > 0.When σ1 is relatively small, g(σ1) may be negative, in such case, the vaccineis high-efficient for subgroup 1, even the vaccine coverage level in subgroup2 is 0, the herd immunity can still be achieved. However, if σ1 is close to 1,g(σ1) is always positive, in such case, we must vaccinate part of subgroup 2to reach the herd immunity. The infection factors σ1 and σ2 can affect the45herd immunity threshold relations.(C). The population mixing threshold value for subgroup 1 is 1−R0(a1 p112α1−R0)(1−σ2) .If this value is smaller than 1, as the vaccine coverage level in subgroup 2reaching to a certain value, the herd immunity can be reached in the wholepopulation. Else if the value is larger than 1, part of subgroup 1 has to bevaccinated in order to achieve the herd immunity.(D). If the vaccine is more efficient for subgroup 1, or subgroup 1 is more activethan subgroup 2, vaccination in subgroup 1 is more crucial. With priority ofvaccinating subgroup 1, it is easier to reach the herd immunity in the wholepopulation.We conclude the analysis of imperfect vaccine case by the following theorem,Theorem 2.3.2. In the two-subgroups heterogeneous mixing population with im-perfect vaccine, if the vaccine coverage levels in two subgroups satisfy the inequal-ity (2.43), the herd immunity can be reached in the whole population. If1−R0(a1 p112α1 −R0)(1−σ2)< 1,as the vaccine coverage level in subgroup 2 reaching to a certain level, withoutvaccinating any individual in subgroup 1, the herd immunity can be achieved inthe whole population. However, if1−R0(a1 p112α1 −R0)(1−σ2)> 1,to reach the herd immunity, part of subgroup 1 has to be vaccinated. It is easier toachieve the herd immunity if the vaccine efficacy is sufficiently high.Remark 7. The herd immunity threshold relation is determined by the basic re-production numberR0, the vaccine efficacy and the mixing pattern of the two sub-groups. To reach the herd immunity threshold relation in the whole population withvaccinating only one group, it is necessary to focus on the more active subgroup.462.3.2.3 Possible extensionsWe defined the herd immunity threshold in a population with two sub-groups andreach to conclusions that the herd immunity may be reached by vaccinating onlyone subgroup under certain circumstances. Some population mixing threshold val-ues are given. Our results indicate one possible effective public health measure,which is the value of ai pii2αi , with i representing the subgroup. If this value is largerthan 1, this subgroup will generate more new infections, then this subgroup is thepriority of the mass vaccination program. Now we can give positive answers tothe questions in the introduction section. And if the subgroup with children andelderly people can be protected better by the indirect protection from the vaccine,the members of this subgroup do not need to undertake the risk either from theinfection or from the vaccination. Some possible extensions with respect to thissubject are:(A). To extend to multi-subgroup population with the number of subgroups largerthan 2. In such case, the deterministic models are more complicated to ana-lyze, it is more difficult to writeRc in terms ofR0. The paper [19] by Brauerand Watmough studied this type of model in detail.(B). It is noticed that the herd immunity threshold relation is related to the natureof the mixing in the population, such as: population size, number of contactsai and fractions of contacts pi j, which we stated in chapter 2. It is interestingto investigate deeper in this direction. It is possible to first consider somesimpler mixing patterns, such as: proportionate mixing, preferred mixing.2.3.3 Economic considerations for the eradication”Herd immunity” is a central conception in disease eradication. An infectious dis-ease can only be eradicated globally if and only if it is eliminated in every countryin the world (see [8] and [9]). To complete our analysis in herd immunity, we studythe paper [8] by Scott Barrett and introduce some definitions and results in dis-ease eradication, which closely relate with mathematical epidemiology and game47theory. We define the critical level of vaccination pc aspc = 1− 1R0,under the assumptions that the vaccine is highly effective and the population ishomogeneous mixing.Two important concepts from epidemiology are:Definition 6. A disease is eliminated in country i if and only if the vaccine cover-age level pi ≥ pc,Definition 7. A disease is eradicated if and only if pi ≥ pc, for i = 1,2, · · · ,N.Where N denotes the number of countries in the world. Clearly, the eradicationof disease asks for the international cooperation. Within a representative country,we assume that the population is homogeneous mixing and for each individual,there are only two options: to get vaccinated or not to get vaccinated. For individ-uals who decide to get vaccinated, the payoff is −c. c is the cost of vaccine, andincludes the probability of an adverse reaction to the vaccine. For individuals whodecide not to get vaccinated, we define the force of infection, can also be be foundin Anderson and May [3],λ = µR0α(pc− p), for pc ≥ p. (2.45)In (2.45), µ is a parameter which can be eliminated by normalization, α representsthe degree to which the disease can be transmitted internationally. It is obvious thatif the vaccine coverage level p exceeds the critical level of vaccination, the force ofinfection λ = 0. We further let b denote the economic cost of an infection, then forindividuals who choose not to get vaccinated, the payoff is Π = −bR0α(pc− p).The strategy of to be vaccinated will be dominant if bR0α(pc− p) ≥ c. We canfind the competitive equilibrium of vaccination:p0 = 1− bα+ cbαR0. (2.46)It follows that:48Proposition 2.3.3. For c > 0 and pc > p0; the disease will not be eliminated in acompetitive equilibrium autarky.Thus, the elimination of diseases requires the public health control. There isone more definition,Definition 8. A disease is controlled in country i if and only if, for c > 0, pc >pi > p0.We turn to part of public policy. The payoff to society as a whole isΠ=−cpn−bR0α(pc− p)(1− p)n, p < pc,orΠ=−cpn, p≥ pc,where n denotes the population of this country.If the vaccine coverage level is lower than the herd immunity threshold pc, tomaximize the payoff to society, we havepu = 1− bα+ c2bαR0. (2.47)To compare the expression in (2.47) with the repression of pc. If bα ≥ c, thenpu ≥ pc. In this condition, the optimal vaccination rate is pc, the herd immunitycan be reached. So policy should aim for elimination, not control.2.4 Conclusion of this chapterIn the first part of this chapter, we provided some basic knowledge of game theoryand deterministic epidemic models.In the second part of this chapter, we have constructed four vaccination games.The first two games are vaccination games in homogeneous mixing population. Inthese two games, by coupling the classic game theory and the deterministic epi-demic models, we were able to extract more accurate estimations of the expectedvaccine coverage levels. In the vaccination game with perfect vaccine, we obtained49the implicit equations of the vaccine coverage level p,pip(1− p) = r.If the function pip is a decreasing function of p, there exists a unique solution. Themonotonicity is the sufficient condition for uniqueness. We will focus on this proofin the following chapter.For the vaccination game with imperfect vaccine, we obtained the formular1−σ = ppiv +(1− p)pip.Similar to the vaccination game with perfect vaccine, the attack ratios piv and pipof the epidemic being monotone decreasing with increased vaccination coverage isthe sufficient condition for the game having unique Nash Equilibrium strategy set.We have also analyzed the vaccination games in heterogeneous mixing popula-tion. The case we considered in this thesis is that the population can be divided intotwo sub-groups based on different activity levels. It is straightforward to constructvaccination games with more than two groups. The extension will generate payoffmatrices with more dimensions. In these two games, we obtained the expressionsof the expected vaccine coverage levels in implicit forms. In the next chapter, wewill examine whether or not the attack ratios in heterogeneous mixing populationare decreasing functions of the vaccine coverage levels.In the third part of this chapter, we first reviewed some previous work of herdimmunity and generalized the herd immunity threshold value in homogeneous mix-ing population. The main point in this chapter is that we defined the herd immunitythreshold relations in heterogeneous mixing population. The definitions can be eas-ily extended to the case of heterogeneous mixing population with more than twosub-groups.50Chapter 3Mathematical analysis ofvaccination modelsIn this chapter, we demonstrate that the attack ratios can be expressed as implicitfunctions in terms of the vaccine coverage levels. We also prove that the attackratios of the epidemics are decreasing functions with increased vaccine coveragelevels. In the first part of this chapter, we review some work in game theory andepidemiology, and show why the property of monotonically decreasing is crucialin proving the uniqueness of Nash equilibria of games. In the second part, weconsider the case of homogeneous mixing population with perfect vaccine and im-perfect vaccine. In the third part, we consider the case of heterogeneous mixingpopulation with perfect vaccine and imperfect vaccine. Because of the large num-ber of parameters in the mathematical model, we obtain only partial results. In thelast part, we consider the age-of-infection epidemic models. The results for SIRand SVIR models can be extended to the general compartmental epidemic models.3.1 Attack ratios and vaccine coverage levels3.1.1 MotivationThe first paper to discuss game theory and epidemiology is [33] by Fine and Clark-son. They discussed how individuals decide whether or not to accept vaccination51based on their perception of the risks involved. The conclusion in the paper is thatrational informed individuals choose a lower vaccine uptake than would the com-munity if it acted as a whole. The topic was also discussed by Brito et al. in [22] in1991. The application of game theory to epidemiology assumes that the decisionsof individuals are completely rational. It has been proved that if individuals aredriven solely by self-interest, the vaccine coverage level of the whole populationwill be suboptimal, and this level is lower than the herd immunity threshold value(see [62], [57], [56], [35], [10], [11]). It has been mentioned that in the simple vac-cination game for an endemic disease in a homogeneous mixing population, thereis a unique Nash equilibrium behavior ([56]). One crucial condition for unique-ness is that the eventual attack ratio of an individual must decrease strictly withthe vaccine coverage level p ([10]). In chapter (2) of this thesis, we demonstratedfour games. The first two games are vaccination games in homogeneous mixingpopulations. In the vaccination game with perfect vaccine, we have the formulapip(1− p) = r, (3.1)with p the vaccine coverage level, r the relative risk, pip the attack ratio. Whenthe vaccine coverage level p = 0, the left side of the equation is pi0 and this valueshould be larger than r. When the vaccine coverage level p ≥ pc, the left side ofthe equation is 0, which is smaller than r. If we differentiate the function on theleft side with respect to p, the derivative ispi ′p(1− p)−pip,We first assume that if pi ′p = 0, then pi0 = pip for all 0 < p≤ 1. Apparently, we havepi0 > 0, which implies pipc > 0 (pc is the herd immunity threshold). However, ifp≥ pc, the attack ratio has to be 0. This contradiction gives us pi ′p 6= 0. The samemethod can show that it is impossible for pi ′p > 0. If pi ′p is negative, pi ′p(1− p)−pipis negative, the formula (3.1) has a unique solution. Thus, it is important to showthat the attack ratio is a decreasing function of the vaccine coverage level p.52In the vaccination game with imperfect vaccine, we have the formular1−σ = ppiv +(1− p)pip, (3.2)with pip the attack ratio for susceptible group and piv the attack ratio for vaccinatedgroup. If the vaccine coverage level p = 0, the right side is pi0, and this valueshould be larger than r1−σ . If the vaccine coverage level p ≥ pc, the right side ofthe equation is smaller than r1−σ . Because the left side of the equation is a constant,if the derivative of the right side is negative, then we have only one solution. Wetake the derivative of ppiv +(1− p)pip with respect to p to getpiv + ppi ′v−pip +(1− p)pi ′p= (piv−pip)+ ppi ′v +(1− p)pi ′p.We have piv−pip < 0, since the vaccinated individuals are protected by vaccine, theattack ratio will be smaller. If both pi ′v and pi ′p are negative, then the formula (3.2)has a unique solution. It can be proved that pip and piv can not be positive or zero.In this chapter, we will focus on the relation between the vaccine coverage leveland the attack ratio.3.1.2 Analysis in homogeneous mixing population3.1.2.1 Case of perfect vaccineWe use the SIR model to analyze the attack ratio,dSdt=−βSI,dIdt= βSI−αI,dRdt= αI,(3.3)53with initial values S0 = (1− p)N and I0 = 0, other parameters are defined as before.The control reproduction numberRc =βαS0 > 1,if and only if the disease can spread among the population. The final size relationcan be expressed by the followinglnS0S∞=βαS0[1− S∞S0 ]. (3.4)If we take the derivative of S∞ with respect to p, we obtaindS∞d p=1S0− βα1S∞− βα(−N). (3.5)Because the control reproduction numberRc > 1, we have1S0− βα< 0.Since we haveI∞ = limt→∞ I(t) = 0,which means that the number of individuals in infective class increases first andthen decreases. Thus, based on the continuity,βS∞I∞−αI∞ < 0βS∞−α < 0.Thus, we haveβα<1S∞.Now we go back to (3.5) to getdS∞d p> 0. (3.6)54We now investigate how the attack ratio will change if we increase or decreasethe vaccine coverage level p. The attach ratio in SIR model is defined aspip =S0−S∞S0.From the final size relation (3.4), we take the derivative of S∞S0 with respect tothe vaccine coverage level p to haved S∞S0d p=βα (1− S∞S0 )βα S0− S0S∞(−N). (3.7)Since we have S0 > S∞, 1− S∞S0 > 0. And from previous analysis we know thatS∞βα < 1, so (βα − 1S∞ )< 0. Thus, we conclude thatd S∞S0d p > 0. The analysis indicatesthatdpipd p< 0. (3.8)We now reach the conclusions on the change of S∞ and the attack ratio with thevaccine coverage level p. We summarize with two theorems.Theorem 3.1.1. In the SIR compartmental model (3.3), when the vaccine coveragelevel p increases, S∞ will also increase.Theorem 3.1.2. In the SIR compartmental model, if the vaccine coverage levelp increases, the attack ratio decreases. The attack ratio is a decreasing functionabout the vaccine coverage level p.3.1.2.2 Case of imperfect vaccineThe SVIR model we use to describe such condition isdSdt=−βSI,dVdt=−σβV I,dIdt= βSI+σβV I−αI,dRdt= αI,(3.9)55with initial values S0 = (1− p)N, V0 = pN, I0 = 0, other parameters are defined asbefore. Since the population is not completely susceptible, it is accurate to use theterminology control reproduction number instead of basic reproduction number.In this case,Rc =βα(S0 +σV0),and the final size relation can be expressed by the following equations,lnS∞S0=βα[S∞+V∞−S0−V0],lnV∞V0=σβα[S∞+V∞−S0−V0].(3.10)It is easy to find out thatV∞V0=(S∞S0)σ, (3.11)which implies that for every individual, staying in group V is much safer thanstaying in group S.If we change the vaccine coverage level p, how this will affect S∞ and V∞? Wetake the derivatives of S∞ and V∞ with respect to p to have(1S∞− βα)dS∞d p=βαdV∞d p− 11− p ,(1V∞− σβα)dV∞d p=σβαdS∞d p+1p.(3.12)We reorganize equations in (3.12) to obtain two similar equation systems (by dif-ferent substitutions)− 1NdS∞d p=βαdV∞d p (− 1N )+ 1S01S∞− βα,− 1N1− βα S∞− σβα V∞S∞V∞dV∞d p=[σβα V0 +βα S0]− S0S∞S0V0,(3.13)56and− 1NdV∞d p=σβαdS∞d p (− 1N )− 1V01V∞− σβα,− 1N1− σβα V∞− βα S∞S∞dS∞d p=V0−V∞[σβα V0 + βα S0]S0V0.(3.14)It is noted that in (3.13) and (3.14),Rc =σβα V0+βα S0, so the key to find whetherthe derivatives are positive or negative is the relation of V∞V0 ,S∞S0and Rc.We first consider the case that if p= 0. Then the SVIR model can be simplifiedto the simple SIR model, the final size relation of such model islnS0S∞=Rc[1− S∞S0 ],where Rc =βα N. To prove there is a unique solution of the final size relation, wedefineg(x) = lnS0x−Rc[1− xS0 ],we have g(0+) > 0, g(S0) < 0, and g′(x) < 0 if and only if 0 < x < S0Rc . Thedisease can only spread ifRc > 1. In such case, g(x) is monotone decreasing froma positive value at x = 0+ to a minimum at x = S0Rc and then increases to a negativevalue at x = S0. Thus, there is an unique zero S∞ of g(x) with S∞ <S0Rc, which isS0S∞>Rc, p = 0. (3.15)Secondly, we consider the case that p = 1. In such case, we have a simplifiedVIR model, the final size relation of such model islnV0V∞=Rc[1− V∞V0 ],where Rc =σβα N. If we use the same method before to prove the uniqueness ofV∞, we haveV0V∞>Rc, p = 1. (3.16)57Now we consider the case that 0 < p < 1. If we have S0S∞ < Rc, the secondequation in (3.13) yields dV∞d p < 0, then from the first equation in (3.13), this leadsto dS∞d p < 0; the second equation in (3.14) yieldsV0V∞>Rc, finally this assumptionleads to V0V∞ >S0S∞. This is a contradiction to S0S∞ >V0V∞. So S0S∞ >Rc must hold.To investigate that how S0S∞ andV0V∞will change if p is changed, we still use thefinal size relation (3.10). For simplicity, we denote pi1 := S∞S0 and pi2 :=V∞V0,lnpi1 =βα(1− p)N[pi1 + p1− ppi2−11− p ],lnpi2 =σβαpN[pi2 +1− pppi1− 1p ].(3.17)We take the derivatives of pi1 and pi2 with respect to p to obtainσpi1,ppi1=pi2,ppi2,(1− βαS∞− σβα V∞)pi1,ppi1=βαN(pi2−pi1).(3.18)Since pi2 > pi1, we have pi1,p > 0 and pi2,p > 0. From the fundamental knowledgeof calculus, we know thatd( S0S∞ )d p< 0,d(V0V∞ )d p< 0, 0 < p < 1.Now we know that both S0S∞ andV0V∞are monotone decreasing as p increases from0 to 1, and for all p, we have S0S∞ >Rc. If when p = 0+,V0V∞>Rc then for all p,we have V0V∞ >Rc; else if when p = 0+,V0V∞<Rc, as p increases, V0V∞ will be firstsmaller thanRc, then will be bigger thanRc.We have two lemmas,Lemma 3.1.3. In the SVIR model (3.9) with infection factor σ . If the epidemiccan spread among the population initially, we have:(A1) βS0 +σβV0 > α ,(A2) βS∞+σβV∞ < α .58Proof. (1). Since the assumption is that the disease can spread among people,the disease-free equilibrium is unstable, we haveRc =σβα V0 +βα S0 > 1,(2). We know that the endemic equilibrium is asymptotically stable, so βS∞+σβV∞ < α .For the simple SIR model, S∞ can be obtained from the final size relation, wedefine σc as follows: ( S0S∞)σc =Rc = βαN. (3.19)Lemma 3.1.4. In the SVIR model (3.9) with infection factor σ , there exists athreshold value σc, which was defined in (3.19). We have:(A). if the infection factor σ is smaller than σc, V0V∞ < Rc <S0S∞holds when thevaccine coverage level is low, and by increasing the vaccine coverage levelp,Rc < V0V∞ <S0S∞(B). if the infection factor σ is bigger than σc, the inequalityRc < V0V∞ <S0S∞holds.The analysis establishes the result on S∞ and V∞.Theorem 3.1.5. In the SVIR compartmental model (3.9), the infection factor σwill affect S∞ and V∞. There are two conditions:(A). If σ is smaller than σc, there exists a threshold value pc, as vaccine coveragelevel p increases from 0 to pc, S∞ will decrease and V∞ will still increase;while p increases from pc to 1, S∞ and V∞ both increase, where pc is thethreshold value as we defined before.(B). If σ is bigger than σc, when the vaccine coverage level p increases, S∞ willdecrease and V∞ will increase.We define the attack ratios for susceptible class and vaccinated class, pis andpiv, respectively,pis = 1− S∞S0 , piv = 1−V∞V0,since pi1,p > 0 and pi2,p > 0, we know thatdpisd p< 0 anddpivd p< 0. So we reach theresult on the change of attack ratios,59Theorem 3.1.6. In the SVIR compartmental model, if the vaccine coverage levelp increases, for vaccinated group, for non-vaccinated group and for the wholepopulation, the attack ratios decrease. The infection factor σ and the attack ratiosare independent.3.1.3 Analysis in heterogeneous mixing population3.1.3.1 Case of perfect vaccineThe two-group SIR vaccination model isS′1 =−a1 p11S1I1N1+a1 p12S1I2N2,I′1 = a1 p11S1I1N1+a1 p12S1I2N2−α1I1,R′1 = α1I1,S′2 =−a2 p21S2I1N1+a2 p22S2I2N2,I′2 = a2 p21S2I1N1+a2 p22S2I2N2−α2I2,R′2 = α2I2,(3.20)with(1). Mean infective period in group i is 1αi , with i = 1,2,(2). rate of infection in group i is σi, with i = 1,2,(3). vaccine coverage levels in each subgroup are v1 and v2, the initial conditionsare S1(0) = (1− v1)N1, S2(0) = (1− v2)N2, I1(0) = 0 and I2(0) = 0.By using the next generation matrix approach which was introduced in [68],we have the expression of the control reproduction numberRc,Rc =a1 p11α1 (1− v1)+a2 p22α2 (1− v2)+(1− v1)(1− v2)√a21 p211α21+a22 p222α22− 2a1a2 p11 p22α1α2 +4a1a2 p12 p21α1α22.60The final sizer relation of (3.20) is:lnS1(0)S1(∞)=a1 p11N1(1− v1)N1−S1(∞)α1+a1 p12N2(1− v2)N2−S2(∞)α2,lnS2(0)S2(∞)=a2 p21N1(1− v1)N1−S1(∞)α1+a2 p22N2(1− v2)N2−S2(∞)α2.(3.21)We notice that in the situation of proportionate mixing, we have(S1(∞)S1(0))a2=(S2(∞)S2(0))a1,For simplicity, we denote:c1 :=a1 p11α1N1, c2 :=a1 p12α2N2, c3 :=a2 p21α1N1, c4 :=a2 p22α2N2.We take the derivatives of S1(∞) and S2(∞) with respect to v1, to find out ifwe only adjust vaccine coverage level in one subgroup, how this will affect bothsubgroups. We have1S1(0)−dS1(∞)dv1(− 1N1 )S1(∞)= c1[1− dS1(∞)dv1 (−1N1)]+ c2[−dS2(∞)dv1 (−1N1)],−dS2(∞)dv1(− 1N1 )S2(∞)= c3[1− dS1(∞)dv1 (−1N1)]+ c4[−dS2(∞)dv1 (−1N1)].Reorganize the above equations to have,[(c1− 1S1(∞))(c4−1S2(∞))−c2c3]dS1(∞)dv1 (−1N1)= (c1− 1S1(0))(c4−1S2(∞))−c2c3,(3.22)[(c1− 1S1(∞))(c4−1S2(∞))− c2c3]dS2(∞)dv1 (−1N1) = c3(1S1(0)− 1S1(∞)).Because we know that S1(∞)< S1(0) and S2(∞)< S2(0), it is easy to haveS1(0)> S1(∞),1S1(0)<1S1(∞),61andS2(0)> S2(∞),1S2(0)<1S2(∞).To determine that dS1(∞)dv1 anddS2(∞)dv1are positive or negative, we need to analyzethe parameters in (3.22).We state the following lemma, and use the lemma to do further analysis.Lemma 3.1.7. In the two-group SIR model (3.23), we have(A),c1S1(0)> 1,c4S2(0)> 1.(B),c1S1(∞)< 1,c4S2(∞)< 1.Proof. (A). We linearize the system near the disease-free equilibrium to have[I′1I′2]=[α1(c1S1(0)−1) α2c2S1(0)α1c3S2(0) α2(c4S2(0)−1)][I1I2].This disease-free equilibrium is unstable, and the disease can spread amongboth sub-populations. Thus, the 2×2 matrix has two positive eigenvalues,α1(c1S1(0)−1)+α2(c4S2(0)−1)> 0,(c1S1(0)−1)(c4S2(0)−1)> c2c3S1(0)S2(0),which leads toc1S1(0)> 1,andc4S2(0)> 1.62(B). We linearize the system near the endemic equilibrium, we have[I′1I′2]=[α1(c1S1(∞)−1) α2c2S1(∞)α1c3S2(∞) α2(c4S2(∞)−1)][I1I2].The endemic equilibrium is asymptotically stable, so this 2× 2 matrix hastwo negative eigenvalues,α1(c1S1(∞)−1)+α2(c4S2(∞)−1)< 0,(c1S1(∞)−1)(c4S2(∞)−1)> c2c3S1(∞)S2(∞),which lead toc1S1(∞)< 1,andc4S2(∞)< 1.Based on the previous analysis and we go back to (3.22). We reach the con-clusion that dS1(∞)dv1 > 0 anddS2(∞)dv1< 0. The following theorem can be obtained,Theorem 3.1.8. In the heterogeneous SIR model, if we increase the vaccine cover-age level v1 in sub-group 1, S1(∞) and S2(∞) will both increase. Same conclusioncan be drawn if we increase the vaccine coverage level v2.The attack ratios in susceptible class and in vaccinated class are defined as:pip = 1− S1(∞)S1(0) ,andpiv = 1− S2(∞)S2(0) .It is observed that when S1(0) decreases, S1(∞) will increase, thus, pip will de-crease. For piS2 , if we decrease S1(0), then S2(0) will be the same, S2(∞) willincrease, thus, piv will decrease. For the case of increase of S2(0), we will have thesame result.63Now we have the theorem,Theorem 3.1.9. In the heterogeneous SIR model, if we increase the vaccine cover-age level in either sub-group, the attack ratios in every sub-group will decrease.3.1.3.2 Case of imperfect vaccineThe two-group SVIR vaccination model isS′1 =−[a1 p11S1I1N1+a1 p12S1I2N2],V ′1 =−σ1[a1 p11V1I1N1+a1 p12V1I2N2],I′1 = [a1 p11S1I1N1+a1 p12S1I2N2]+σ1[a1 p11V1I1N1+a1 p12V1I2N2]−α1I1,R′1 = α1I1,S′2 =−[a2 p21S2I1N1+a2 p22S2I2N2],V ′2 =−σ2[a2 p21V2I1N1+a2 p22V2I2N2],I′2 = [a2 p21S2I1N1+a2 p22S2I2N2]+σ2[a2 p21V2I1N1+a2 p22V2I2N2]−α2I2,R′2 = α2I2,(3.23)with(1). Mean infective period in group i is 1αi , with i = 1,2,(2). rate of infection in group i is σi, with i = 1,2,(3). vaccine coverage levels in each subgroup are v1 and v2, the initial conditionsare S1(0) = (1− v1)N1, S2(0) = (1− v2)N2, V1(0) = v1N1, V2(0) = v2N2,I1(0) = 0 and I2(0) = 0.64We use the next generation matrix approach, which is described in [68], to find thecontrol reproduction numberRc. The final size relation of (3.23) islnS1(0)S1(∞)=a1 p11α1N1[N1−S1(∞)−V1(∞)]+ a1 p12α2N2 [N2−S2(∞)−V2(∞)],lnV1(0)V1(∞)= σ1[a1 p11α1N1[N1−S1(∞)−V1(∞)]+ a1 p12α2N2 [N2−S2(∞)−V2(∞)]],lnS2(0)S2(∞)=a2 p21α1N1[N1−S1(∞)−V1(∞)]+ a2 p22α2N2 [N2−S2(∞)−V2(∞)],lnV2(0)V2(∞)= σ2[a2 p21α1N1[N1−S1(∞)−V1(∞)]+ a2 p22α2N2 [N2−S2(∞)−V2(∞)]].(3.24)Two equations are revealed from (3.24),(S1(0)S1(∞))σ1 =V1(0)V1(∞),(S2(0)S2(∞))σ2 =V2(0)V2(∞).(3.25)We take the derivatives of S1(∞), V1(∞), S2(∞) and V2(∞) with respect to v1,to find out if we only adjust vaccine coverage level in one subgroup, how this willaffect two subgroups.− 11− v1 −S1,v1(∞)S1(∞)= c1[−S1,v1(∞)−V1,v1(∞)]+ c2[−S2,v1(∞)−V2,v1(∞)],1v1− V1,v1(∞)V1(∞)= σ1c1[−S1,v1(∞)−V1,v1(∞)]+σ1c2[−S2,v1(∞)−V2,v1(∞)],−S2,v1(∞)S2(∞)= c3[−S1,v1(∞)−V1,v1(∞)]+ c2[−S2,v1(∞)−V2,v1(∞)],−V2,v1(∞)V2(∞)= σ2c3[−S1,v1(∞)−V1,v1(∞)]+σ2c4[−S2,v1(∞)−V2,v1(∞)],(3.26)withSi,v1(∞) :=d(Si(∞))dv1, Vi,v1(∞) :=d(Vi(∞))dv1(i = 1,2).65From (3.26) we haveV1,v1(∞)V1(∞)= σ1S1,v1(∞)S1(∞)+σ11− v1 +1v1,V2,v1(∞)V2(∞)= σ2S2,v1(∞)S2(∞).(3.27)The above relation implies that V2,v1(∞) and S2,v1(∞) are both positive or negative,which indicates the change of vaccine coverage level in subgroup 1 having sameimpact to classes V2 and S2 in subgroup 2.We state the following lemma, and use the lemma to do further analysis.Lemma 3.1.10. In the SVIR heterogeneous model (3.23), we have:(A1)c1S1(0)+σ1c1V1(0)> 1,c4S2(0)+σ2c4V2(0)> 1,(3.28)(A2)c1S1(∞)+σ1c1V1(∞)< 1,c4S2(∞)+σ2c4V2(∞)< 1,(3.29)c1, c2, c3 and c4 are defined before.Proof. (1). We linearize the system near the disease-free equilibrium to have[I′1I′2]=[α1(c1S1(0)+σ1c1V1(0)−1) α2(c2S1(0)+σ1c2V1(0))α1(c3S2(0)+σ2c3V2(0)) α2(c4S2(0)+σ2c4V2(0)−1)][I1I2].This equilibrium is unstable, the 2× 2 matrix has two positive eigenvaluesbecause the disease can spread in two subgroups.α1(c1S1(0)+σ1c1V1(0)−1)+α2(c4S2(0)+σ2c4V2(0)−1)> 0,(c1S1(0)+σ1c1V1(0)−1)(c4S2(0)+σ2c4V2(0)−1)> (c2S1(0)+σ1c2V1(0))(c3S2(0)+σ2c3V2(0)),which lead toc1S1(0)+σ1c1V1(0)> 1,66andc4S2(0)+σ2c4V2(0)> 1;(2). We linearize the system near the final disease-free equilibrium, we have[I′1I′2]=[α1(c1S1(∞)+σ1c1V1(∞)−1) α2(c2S1(∞)+σ1c2V1(∞))α1(c3S2(∞)+σ2c3V2(∞)) α2(c4S2(∞)+σ2c4V2(∞)−1)][I1I2].The endemic equilibrium is asymptotically stable, so this 2× 2 matrix hastwo negative eigenvalues,α1(c1S1(∞)+σ1c1V1(∞)−1)+α2(c4S2(∞)+σ2c4V2(∞)−1)< 0,(c1S1(∞)+σ1c1V1(∞)−1)(c4S2(∞)+σ2c4V2(∞)−1)> c2c3(S1(∞)+σ1V1(∞))(S2(∞)+σ2V2(∞)),which lead toc1S1(∞)+σ1c1V1(∞)< 1,andc4S2(∞)+σ2c4V2(∞)< 1.We now substitute expressions in (3.27) into the first and the third equationsin (3.26),( 1S1(∞)− c1(1+ σ1V1(∞)S1(∞) ))S1,v1(∞)− c2(1+σ2V2(∞)S2(∞))S2,v1(∞)= c1(σ1V1(∞)1− v1 +V1(∞)v1)− 11− v1 ,( 1S2(∞)− c4(1+ σ2V2(∞)S2(∞) ))S1,v1(∞)− c3(1+σ1V1(∞)S1(∞))S1,v1(∞)= c3(σ1V1(∞)1− v1 +V1(∞)v1).(3.30)67We extract the expression of∂S2∂v1,(( 1S2(∞)− c4(1+σ2V2(∞)S2(∞)))(( 1S1(∞)− c1(1+σ2V1(∞)S1(∞)))− c2c3(1+σ1V1(∞)S1(∞))(1+σ1V2(∞)S2(∞)))∂S2∂v1= c3( V1(∞)v1S1(∞)− 11− v1).(3.31)Based on the equations in (3.25), we know that the right hand side of (3.31) ispositive. And from the second part of proof of previous lemma, we know that theleft hand side of (3.31) is also positive, which gives that∂S2∂v1> 0. (3.32)The above analysis implies that if the vaccine coverage level in subgroup 1 is in-creased, for subgroup 2, both S2(∞) and V2(∞) will increase. We now investigatehow S1(∞) and V1(∞) will be affected, and how the change will relate to the vac-cine efficacy. Through some simple calculations, we have the left hand side of theexpression of∂S1∂v1be1S1(∞)S2(∞)[1− c4(S2(∞)+σ2V2(∞))− c1(S1(∞)+σ1V1(∞))+(c1c4− c2c3)(S2(∞)+σ2V2(∞))(S1(∞)+σ1V1(∞))].(3.33)The right hand side of the expression of∂S1(∞)∂v1isN1S2(∞)[(c2c3− c1c4)(S2(∞)+σ2V2(∞))(σ1V1(∞)S1(0)+V1(∞)V1(0))+ c1(σ1V1(∞)S1(0)+V1(∞)V1(0))+c4S1(0)(S2(∞)+σ2V2(∞))− 1S1(0) ].(3.34)68By observing the left hand sides of (3.31) and (3.33), they have the same form, sowe only need to analyze the positive or negative of (3.34), to rewrite (3.34) asN1S2(∞)[(c2c3− c1c4)(S2(∞)+σ2V2(∞))(σ1V1(∞)S1(0)+V1(∞)V1(0))+ c1(σ1V1(∞)S1(0)+V1(∞)V1(0))+c4S1(0)(S2(∞)+σ2V2(∞))− 1S1(0) ]=N1S2(∞)S1(0)[(c2c3− c1c4)(S2(∞)+σ2V2(∞))(S1(∞)+σ1V1(∞))+ c1(σ1V1(∞)+S1(∞))+ c4(σ2V2(∞)+S2(∞))−1]+N1S2(∞)S1(0)[(c2c3− c1c4)(S1(∞)+σ1V1(∞))+ c1](S1(0)V1(∞)V1(0)−S1(∞)).(3.35)If σ1 = 1, thenS1(0)V1(∞)V1(0)= S1(∞),since the vaccine has no effect. In such case, we have∂S1(∞)∂v1=−N1S1(∞)S1(0)< 0. (3.36)Due to continuity, we have that when the vaccine is low-effective, σ1 is close to 1,then if the vaccine coverage level in subgroup 1 is increased, S1(∞) will decrease.If we substitute (3.36) into the first equation of (3.27), we have1V1(∞)∂V1(∞)∂v1= σ11S1(∞)∂S1(∞)∂v1+σ11− v1 +1v11V1(∞)∂V1(∞)∂v1=−σ1 1S1(∞)N1S1(∞)S1(0)+σ11− v1 +1v1∂V1(∞)∂v1=V1(∞)v1> 0.This is the case that σ1 = 1, which vaccine has the lowest efficacy. If the vaccineefficacy is upgraded, people get better protection, so∂V1(∞)∂v1> 0,69holds for all 0 < σ < 1.If σ1 = 0, which is the case that the vaccine is perfect effective, thenN1S2(∞)[(c2c3− c1c4)(S2(∞)+σ2V2(∞))(σ1V1(∞)S1(0)+V1(∞)V1(0))+ c1(σ1V1(∞)S1(0)+V1(∞)V1(0))+c4S1(0)(S2(∞)+σ2V2(∞))− 1S1(0) ]=N1S2(∞)S1(0)[(c2c3− c1c4)(S2(∞)+σ2V2(∞))(σ1V1(∞)+V1(∞)S1(0)V1(0))+ c1(σ1V1(∞)+V1(∞)S1(0)V1(0))+ c4(S2(∞)+σ2V2(∞))−1]=N1S2(∞)S1(0)[(c2c3− c1c4)(S2(∞)+σ2V2(∞))(σ1V1(0)+S1(0))+ c1(σ1V1(0)+S1(0))+ c4(S2(∞)+σ2V2(∞))−1].(3.37)If σ2 is also close to 0, then (3.37)≈ S1S2(∞)N1(0)[(c2c3− c1c4)(S2(0)+σ2V2(0))(σ1V1(0)+S1(0))+ c1(σ1V1(0)+S1(0))+ c4(S2(0)+σ2V2(0))−1].From previous lemma, we know that∂S1(∞)∂v1> 0.We establish the conclusion,Theorem 3.1.11. In the heterogeneous SVIR model with imperfect vaccine, if wekeep the vaccine coverage of subgroup 2 in a steady level, and only change thevaccine coverage level of subgroup 1, then:(1). When the vaccine has low effectivity to subgroup 1, if vaccine coverage levelof subgroup 1 increases, then S1(∞) will decrease, and V1(∞) will increase;both S2(∞) and V2(∞) increase.(2). When the vaccine has high effectivity to subgroup 1, if vaccine coverage levelof subgroup 1 increases, S1(∞), V1(∞), V1(∞) and V2(∞) will increase.70Remark 8. We notice that in (3.34), ifσ1V1(∞)+V1(∞)S1(0)V1(0)=1− c4(S2(∞)+σ2V2(∞))(c2c3− c1c4)(S2(∞)+σ2V2(∞))+ c1 , (3.38)there is the change of dynamic behavior of S1(∞), we call (3.38) the critical-relation equation. we know that both σ1V1(∞)+V1(∞)S1(0)V1(0)and 1− c4(S2(∞)+σ2V2(∞)) are positive. If c2c3−c1c4 < 0, (c2c3−c1c4)(S2(∞)+σ2V2(∞))+c1 maybe negative. In such case, the critical-relation equation does not exit, which meansthe relation between σ1 and σ2 does not affect the change of S1(∞), only the valueof σ1 matters. This conclusion is reasonable, because if c2c3− c1c4 < 0 holds, wehave c1 and c4 much bigger than c2 and c3, that means two subgroups have fewcontacts with each other. To the extent, we have two separate homogeneous mixingpopulation. If c2c3− c1c4 is close to 0, or c2c3− c1c4 ≥ 0, the critical-relationequation may exist, that means two different subgroups have frequent contact, thenthe relation of σ1 and σ2 will affect the results.For attack ratios, deriving from the final size relation, we have the conclusion,Theorem 3.1.12. In the heterogeneous SVIR model, if the vaccine coverage levelin subgroup 1 increases, all attack ratios decrease; for the increase of vaccinecoverage level in subgroup 2, the results hold. The infection rates σ1 and σ2 areunrelated to the change of attack ratios.3.2 Some extensions to the age-of-infection modelsFor general compartmental models which may include Exposed class, Quarantineclass, Treatment class, etc., we can always use the general age-of-infection models.By proving the similarity of final size relations obtained from the age-of infectionmodels and from the SIR or SVIR models in homogeneous mixing population andheterogeneous mixing population, we can extend our results in previous sections tothe general compartmental models. Same as in [19], we consider a population witha distribution of subpopulation identified by the variable ζ ranging over a ’state’space Ω, having sizes N(t,ζ ), respectively, at time t. Suppose that each groupζ member makes a(ζ ) contacts sufficient to transmit infection per unit time. We71assume that the fraction of contacts made by a member of group ζ with a memberof group η is p(ζ ,η), then: ∫Ωp(ζ ,η)dη = 1,for each ζ . The total number of contacts made in unit time by member of group ζwith members of group η isa(ζ )p(ζ ,η)N(t,ζ ),and by balance relation we havea(ζ )p(ζ ,η)N(t,ζ ) = a(η)p(η ,ζ )N(t,η),with ζ 6= η .We assume that there is no permanent movement between groups and that thereare no disease deaths, so that N(t,ζ ) is a constant function N(ζ ) of t for all ζ .(A). A multi-group age of infection model with general mixing would beS′(t,ζ ) =−a(ζ )S(t,ζ )∫Ωp(ζ ,η)φ(t,η)N(η)dηφ(t,ζ ) = φ0(t,ζ )+∫ ∞0[−S′(t− τ,ζ )]A(τ,ζ )dτ.(3.39)The initial conditions are S0(ζ ) = (1− p(ζ ))N(ζ ) and the vaccine is per-fectly effective.It is possible to obtain a system of final size relations, which islnS(0,ζ )S(∞,ζ )= a(ζ )∫Ωp(ζ ,η)N(η)(∫ ∞0φ0(t,η)dt+(S(0,η)−S(∞,η))∫ ∞0A(η)dt)dη .(3.40)(B). We now investigate the case which the vaccine are imperfect, and we modelthis by including factors σ(ζ ), where 0 < σ(ζ )< 1. The multi-group age of72infection model with general mixingS′(t,ζ ) =−a(ζ )S(t,ζ )∫Ωp(ζ ,η)φ(t,η)N(η)dηV ′(t,ζ ) =−σ(ζ )a(ζ )V (t,ζ )∫Ωp(ζ ,η)φ(t,η)N(η)dηφ(t,ζ ) = φ0(t,ζ )+∫ ∞0[−S′(t− τ,ζ −V ′(t− τ,ζ )]A(τ,ζ )dτ.(3.41)The final size relations of the model (3.41) arelnS(0,ζ )S(∞,ζ )= a(ζ )∫Ωp(ζ ,η)N(η)(∫ ∞0φ0(t,η)dt +(N(η)−S(∞,η)−V (∞,η))∫ ∞0A(η)dt)dη ,lnV (0,ζ )V (∞,ζ )= σ(ζ )a(ζ )∫Ωp(ζ ,η)N(η)(∫ ∞0φ0(t,η)dt +(N(η)−S(∞,η)−V (∞,η))∫ ∞0A(η)dt)dη .(3.42)Remark 9. The detailed work of how to generate the final size relations can befound in Appendix A.If we compare the final size relations (3.40) and (3.42) with previous formulasof final size relations, it is clear that the forms of final size relations are similar. Thefinal size relation to the homogeneous mixing age-of-infection model is straight-forward. The generality of the final size relation formulas guarantee that we caneasily extend our results from previous sections to the general compartmental mod-els. More details of similar work can be found in [47].To set an example, we consider the standard SEIR model and show this modelcan be expressed in the form of age-of-infection model. The SEIR model is:S′ =−βS(I+ εE)E ′ = βS(I+ εE)−κEI′ = κE−αIR′ = αI.(3.43)Compartment E is the exposed class of individuals having infectivity reduced by afactor ε . We assume there is no removed members at the start of the epidemic andthere are no disease deaths, the population size is a constant N.73We first define:φ = εE + I.As we have shown in the introduction section, we need to determine the kernelA(τ) of the age-of-infection model. We denote u(τ) be the fraction of infectedmembers with infection age τ who are not yet infective and v(τ) the fraction ofinfected members who are infective. The rate at which members become infectiveat infection age τ is κu(τ). We haveu′(τ) =−κu(τ), u(0) = 1v′(τ) = κu(τ)−αv(τ), v(0) = 0.(3.44)The solution of this ODE system isu(τ) = e−κτ ,v(τ) =κκ−α [e−ατ − e−κτ ].(3.45)Thus, we have the expression of A(τ):A(τ) = εe−κτ +κκ−α [e−ατ − e−κτ ].From this example, it is clearly that we are able to express the standard SEIRmodel in the form of age-of-infection model. Based on the analysis that the formsof final size relations are similar. It is straightforward to have the conclusion that,for an epidemic such as influenza (the infected individuals undergo a period with noinfectivity or partial infectivity, then enter into the infectious class), the vaccinationgame still has a unique Nash equilibrium strategy set. Since the attack ratios forsusceptible individuals and for vaccinated individuals are both decreasing functionsof the vaccine coverage levels.3.3 Conclusion of this chapterMotivated by the analyses of vaccination games in previous chapter, it is neces-sary to prove the uniqueness of Nash equilibria of vaccination games. The mostimportant part of the proof is to show that the attack ratio can be expressed as the74function of the vaccine coverage level, and decreases as the vaccine coverage levelincreases. We investigated all four cases. For the cases of homogeneous mixingpopulation, through the SIR model and the SVIR model, we obtained the completeresults. For the cases of heterogeneous mixing population, through the two-groupSIR model and SVIR model, we were able to obtain partial results. In the lastpart of this chapter, we analyzed the age-of-infection models. The purpose was toprove, for age-of-infection models, the attack ratios are decreasing functions of thevaccine coverage levels. The generalization of the final size relations of the age-of-infection models indicates the similarity comparing with the final size relations ofthe classic SIR model and the SVIR model, for both cases of homogeneous mixingpopulation and of heterogeneous mixing population. Thus, we can further extendprevious results on attack ratios to the case of epidemic with more complicatedstages, such as influenza.Remark 10. The case of heterogeneous mixing population with imperfect vaccineactually includes the other three cases. The reasons we organized this chapter inthe above structure are: firstly, the results in this case is partial, we can not provecompletely; secondly, this layout is consistent with the structure in chapter 2. Theresults in the other three cases are complete.75Chapter 4Malaria vaccination modelIn this chapter, we apply the deterministic epidemic models to study the dynamicsof transmission of malaria. In the first part, we demonstrate some facts about thevector-borne diseases and the disease of malaria. We also review some previouswork on mathematical modelling, especially some work on modelling the trans-mission of malaria. We then formulate and analyze a delayed malaria transmissionmodel, which includes both extrinsic incubation period and intrinsic incubationperiod. We also conduct several numerical experiments to verify our theoreticalresults. In the last part of this chapter, we formulate a malaria vaccination model,calculate the control reproduction number Rc, analyze the stability of the disease-free equilibrium. We then propose the threshold condition of eradication.4.1 General frame of modelling the vector-borne disease4.1.1 Introduction to vector borne diseasesAccording to the statistical data from CDC, vector-borne diseases account for morethan 17% of all infectious diseases. Vector-borne diseases cause more than 1 mil-lion deaths per year. Some major types of vector-borne diseases include: malaria(most prevalent vector-borne disease whose vectors are the mosquitoes), Denguefever (second most prevalent mosquito-borne disease globally), Lyme disease, WestNile disease (which is caused by West Nile virus), yellow fever, Zika (Zika virus76can be spread through sexual transmission, this type of vector-borne disease differsfrom the other types because vectors and hosts have similar life spans). A relativelycomplete list of vector-borne diseases is shown in the following tables and can alsobe found in [36].Table 4.1: Table 1 of vector-borne diseases (cite from [36])Disease (Pathogen) Animal reservoirs GeographicaldistributionVectorPathogen: Viruses TogaviridaeChikungunya Primates, humans Africa, Asia MosquitoesRoss River fever Marsupials, humans Australia, SouthPacificMosquitoesMayaro Birds South America MosquitoesOnyong-nyong fever Not known Africa MosquitoesSindbis fever Birds Asia, Africa,Australia, Eu-rope, AmericasMosquitoesEastern equine encephalomyelitis Birds Americas MosquitoesWestern equine encephalomyelitis Birds, rabbits Americas MosquitoesBarmah Forest Not known Americas MosquitoesPathogen: Viruses FlaviviridaeDengue fever (serotypes 1-4) Primates, humans tropical coun-triesMosquitoesYellow fever Primates, humans Africa, SouthAmericaMosquitoesKyasanur Forest disease Primates, rodents, camels India, Saudi Ara-biaTicksOmsk haemorrhagic fever Rodents Asia TicksJapanese encephalitis Birds Asia MosquitoesMurray Valley encephalitis Birds Australia MosquitoesRocio Birds South America MosquitoesSt. Louis encephalitis Birds Americas MosquitoesWest Nile encephalitis Birds Asia, Africa,North America,EuropeMosquitoesTick-borne encephalitis Rodents Europe, Asia Ticks77Table 4.2: Table 2 of vector-borne diseases (cite from [36])Disease (Pathogen) Animal reservoirs GeographicaldistributionVectorPathogen: Viruses BunyaviridaeSandfly fever Not known Europe, Africa,AsiaSandfliesRift Valley fever Not known Africa MosquitoesLa Crosse encephalitis Rodents North America MosquitoesCalifornia encephalitis Rodents North America,Europe, AsiaMosquitoesCrimean-Congo haemorrhagic fever Rodents, sheep Europe, Asia,AfricaTicksOropouche fever Not known Central andSouth AmericaMidges,mosquitoesPathogen: Viruses RhabdoviridaeVesicular stomatitis Cattle, horses, pigs Global Phlebotomusflies, mosquitoesBluetongue Cattle, sheep, goats Global Culicoides fliesIn the late 19th century, it was discovered that the mosquitoes are able to trans-mitting diseases as a vector. Since then, arthropods have been connected with manyother diseases ([36], [51]). The most common type of transmission by vectors isbiological transmission, which includes four types: propagative transmission, cy-clopropagative transmission, cyclodevelopmental transmission, vertical and directtransmission. There are plenty of complex factors that can influence the process oftransmission, such as climate, landscape, immigration of population.78Table 4.3: Table 3 of vector-borne diseases (cite from [36])Disease (Pathogen) Animal reservoirs GeographicaldistributionVectorPathogen: BacteriaPlague (Yersinia pestis) Rodents Global FleasTularaemia (Francisella tularen-sis)Rabbits, rodents North America,Europe, AsiaTicks, tabanid fliesQ fever (Rickettsia) Ungulates Global TicksRocky Mountain spotted fever(Rickettsia rickettsii)Rabbits, rodents, dogs Western hemi-sphereTicksMurine typhus (Rickettsia typhi) Rats Global TicksBoutonneuse fever (Rickettsiaconor)Dogs, rodents Europe, Africa TicksQueensland tick typhus (Rick-ettsia australis)Rodents Australia TicksSiberian tick typhus (Rickettsiasiberica)Rodents Asia TicksScrub typhus (Orientia tsutsuga-mushi)Rodents Asia, Australia MitesRelapsing fever (Borrelia species) Rodents Global Ticks and liceLyme disease (Borrelia burgdor-feri)Rodents North America,EuropeTicks4.1.2 Mathematical modelling of malariaMalaria is a vector-borne infectious disease, it is transmitted by infected femaleAnopheles mosquitoes. Malaria can cause symptoms like fever, fatigue, vomitingand headache. It may cause death in severe cases. According to the statistics ofCDC, there are 3.4 billion people live in areas at risk of malaria transmission in106 countries and territories (Centers for Disease Control and Prevention, [23]).It is estimated that in 2013 malaria caused 198 million clinical cases, and morethan 500,000 deaths by WHO (Centers for Disease Control and Prevention, [23]).The malaria-endemic countries mainly include countries in Sub-Saharan Africa,and other countries in south America, south Asia, southeast Asia, Middle-east andparts of Oceania.79Table 4.4: Table 4 of vector-borne diseases (cite from [36])Disease (Pathogen) Animal reservoirs GeographicaldistributionVectorPathogen: ProtozoaMalaria (Plasmodium spp.) Primates, humans Global Anopheline mosquitoesAfrican trypanosomiasis (Try-panosoma rhodesiense)Ungulates Africa Glossina flies(Trypanosoma gambiense) Pigs, ungulates Africa Glossina fliesNagana (Trypanosoma brucei) Ungulates Africa Glossina fliesChagas disease (Trypanosomacruzi)Dogs, cats, opossum Western hemi-sphereTriatomid bugsLeishmaniasis (Leishmania spp.) Dogs, rodents Asia, Africa, Eu-rope, Central andSouth AmericaPhlebotomus fliesPiroplasmosis (Babesia spp.) Ungulates Global TicksPathogen: FilariaBancroftian filariasis (Wuchereriabancrofti)Humans Global MosquitoesBrugian filariasis (Brugia malayi) Humans, primates Asia MosquitoesThe disease is caused by parasitic protozoans of genus Plasmodium. Thereare four species of Plasmodium that can infect humans, which are: P. falciparum,P. vivax, P. ovale and P. malariae. Among them, P. falciparum is responsible formost of the clinical cases and deaths. Another species P. Knowlesi rarely causesinfection in humans. The life cycle of malaria parasite is shown in the picture 4.1(from CDC).Mathematical modeling is one of the most powerful tools to study the transmis-sion of malaria, and is essential for designing the effective public health measuresto prevent the transmission of the disease. Many of the mathematical work onmalaria arises in the study of malaria begun by R.A.Ross ([59]) in 1911. Dr.Rosswas awarded the second Nobel Prize in Medicine for his demonstration of the dy-namics of the transmission of malaria between humans and mosquitoes. In 1957,George Macdonald integrated the study of malaria transmission ([48]). Since then,malaria has been studied from many angles of mathematical epidemiology. A sci-80Figure 4.1: The life cycle of malaria parasite (cite from CDC)entific review on mathematical models of malaria [50] by Mandal et al. elaboratedcurrent mathematical research on malaria transmission.The so-called Ross-Macdonald model is defined as follows (see [1]):x′ =abMNy(1− x)− rxy′ = ax(1− y)−µy(4.1)where the variables and parameters are stated in table 4.5.The threshold value (or basic reproduction number) R0 =ma2bµr. Ross arguedthat in order to eradicate malaria, the only thing we need to do is reducing thenumber of female Anopheles mosquitoes below a critical threshold level.In 1957, George Macdonald modified the model (4.1) by including the incu-81Table 4.5: Table of variables and parameters for Ross-Macdonald modelVariable/Parameterx proportion of the human population being in-fectedy proportion of the mosquito population being in-fectedN sizes of population of humansM sizes of population of mosquitoesr recovery rate in human. 1r is the average periodfor human to stay in infective classµ death rate in mosquito. 1µ is the expected lifespanof mosquitoesa biting rate on human by a single mosquitob proportion of infected bites on human that pro-duce an infectionbation period τ in mosquitoes, the model is described as (see [1]):x′ =abMNy(1− x)− rxz′ = ax(1− y− z)−axˆ(1− yˆ− zˆ)e−µτ −µzy′ = axˆ(1− yˆ− zˆ)e−µτ −µy(4.2)withxˆ = x(t− τ), yˆ = y(t− τ), zˆ = z(t− τ).And z(t) is the fraction of exposed mosquitoes at time t, y(t) is the fraction ofinfectious mosquitoes at time t. The threshold value of this model (4.2) isR0 =ma2be−µτµr.Important qualitative conclusion that imagicides are more effective than larvicidescan be drawn from this expression of basic reproduction number.Other aspects of mathematical modeling of malaria transmission include:• Heterogeneity in human populations. Anderson and May introduced an age82structure model in [3]. They considered the human density as a function ofage and time.• Immunity. People may obtain the immunity from continuous exposure to thedisease. Malaria mainly causes deaths in children who are under 5 years old.Koella described the importance of introduction of immunity in mathemati-cal models in [45]. Immunity is included in the mathematical model in twoways ([50]), one is to introduce a separate immune class (see [28], [6]), theother is to incorporate an immunity function (details can be found in [50]).• Impact of environmental change, such as global warming and so on. The en-vironmental factors can change the living pattern of the mosquitoes, furtherwill change the patterns of vector-borne diseases. Some papers discussed theimpact from this angle, (see [40] and [58]).• Impact of social and economic factors, such as: migration. Several papershave discussed the impact from this viewpoint (see [61], [7] and [67]). It isnoticed that the malaria candidate vaccine, Mosquirix (trademark name ofRTS,S), is likely to be put on market soon. It is possible to incorporate thevaccination class in the model. Furthermore, the introduction of game theorycan be very useful to evaluate/predict the expected attack ratios (population-wise and individual-wise).• Some stochastic models (see [25] and [55]).4.2 Dynamical analysis of a delayed malariatransmission modelIn this section, we concentrate on the impact of incubation periods in humans and inmosquitoes. This subject has been discussed in several papers. George Macdonaldincorporated the incubation period in mosquitoes. Other models which are relatedcan be found in [1] and [3].834.2.1 Previous workAnother paper [60] by Ruan et al. formulated the malaria transmission model withincubation periods both in humans and in mosquitoes, they analyzed the stabilityof equilibria and conducted the sensitivity analysis of R0. The so-called delayedRoss-Macdonald model for malaria transmission is (see [60])dXdt=−rX(t)+ aHb[H−X(t− τ1)]Y (t− τ1)e−rτ1 ,dYdt=−µY (t)+ aHcX(t− τ2)[M−Y (t− τ2)]e−µτ2 ,(4.3)with X(t) and Y (t) the numbers of infected humans and infected mosquitoes attime t. H and M are the population sizes of humans and mosquitoes, respectively.And τ1 and τ2 are the incubation periods in humans and in mosquitoes, respec-tively. This model was not able to describe the spread of malaria accurately. Webriefly point out two major problems. The first problem of this model (4.3) is theexpressions of cases of new infection. It is assumed that the human population isdivided into susceptible class and infectious class, and the mosquito population isdivided into susceptible class, exposed class and infectious class. The cases of newinfection aH b[H −X(t − τ1)]Y (t − τ1) is not correct, since Y (t) includes both ex-posed mosquitoes and infectious mosquitoes. However, only the contact betweeninfectious mosquitoes and susceptible humans can cause new infection. The sec-ond problem is the result of the stability at endemic equilibrium.The corresponding re-scaled system isdxdt=−rx(t)+abm[1− x(t− τ1)]y(t− τ1)e−rτ1 ,dydt=−µy(t)+acx(t− τ2)[1− y(t− τ2)]e−µτ2 ,(4.4)with x(t) and y(t) the proportions of infected humans and infected mosquitoesat time t, respectively; and m the number of mosquitoes per human. The basicreproduction number to this model is:R0 =a2bcme−rτ1e−µτ2rµ.84The trivial disease-free equilibrium is (0,0) and the possible endemic equilibriumis ( R0−1R0 +ace−µτ2µ,R0−1R0 +abme−rτ1r).The author then proved that ifR0 and a2bcm≤ acr+2ru, this endemic equilibriumis stable in [60]. Apparently the condition a2bcm≤ acr+2ru makes no sense. Andnumerically it was showed that the lengths of incubation periods do not affect thestability, which is also a contradiction to the theoretical results. Other discussionson the results of this paper [60] can be found in [51] by Martcheva et al., they alsomentioned that the special but important case when only the extrinsic incubationperiod is taken into account was not discussed.We formulate a new delayed malaria transmission model in the next part.4.2.2 A modified delayed malaria transmission modelWe consider the human population and mosquito population are both divided intothree compartments: susceptible class, exposed class and infective class. The in-fective class of mosquito population is defined as the mosquitoes with sporozoitesin the salivary glands, and the exposed class of mosquito population is the cate-gory of mosquitoes that are infected, but have not yet developed sporozoites (see[45]). The epidemiological translations of exposed class and infective class of hu-man population can be made in a similar way. Only contacts between susceptiblehumans (mosquitoes) and infective mosquitoes (humans) cause cases of new in-fections in humans (mosquitoes). All new-borns are susceptible to the disease.Humans recover from the infection without immunity. For simplicity, we considerthe sizes of total populations of humans and mosquitoes are constants and denoteNh and Nv to represent. With the time scale short enough, we ignore the birth rateand death rate in human population, and the disease death is not considered. Wehave the relationΛv = µvNvto keep the size of population of mosquito constant.854.2.2.1 The mathematical modelWe have the malaria transmission model including two incubation periods:S′h =−bThvNhShIv + rIh + rEh,E ′h =bThvNhShIv− bThvNh Sh(t− τh)Iv(t− τh)e−rτh− rEh,I′h =bThvNhSh(t− τh)Iv(t− τh)e−rτh− rIh,S′v = Λv−bTvhNhSvIh−µvSv,E ′v =bTvhNhSvIh− bTvhNh Sv(t− τv)Ih(t− τv)e−µvτv−µvEv,I′v =bTvhNhSv(t− τv)Ih(t− τv)e−µvτv−µvIv.(4.5)Where the variables and parameters are stated in table 4.6.Rate of infection bThvNh ShIv in (4.5) is actually bNvNhThvSh IvNv , where bNvNhis thetotal bites on one human in unit time, b NvNh Thv is the total infection caused by thebites in unit time, this value is equivalent with the contact rate. The terms e−rτhand e−µvτv are survival functions of humans and mosquitoes, respectively. Thedifference is that a fraction e−rτh of exposed humans will recover from infectionand go to the susceptible class, while a fraction of e−µvτv of exposed mosquitoeswill die. It is noted that the infected individuals may recover from the diseaseduring the incubation period (see [63]; also discussed in [60]).We now define the proportions of exposed humans, infectious humans, exposedmosquitoes and infectious mosquitoes:x(t) =Eh(t)Nh, y(t) =Ih(t)Nh, z(t) =Ev(t)Nv, w(t) =Iv(t)Nv. (4.6)86Table 4.6: Table of variables and parameters for delayed malaria modelVariable/ParameterSh size of susceptible class in human populationEh size of exposed class in human populationIh size of infected class in human populationSv size of susceptible class in mosquito populationEv size of exposed class in mosquito populationIv size of infected class in mosquito populationNh size of human populationNv size of mosquito populationΛv birth rate in mosquitosµv death rate in mosquitoesr recovery rate in humansb biting rate, the number of bites made by amosquito in unit timeThv probability that a bite by a mosquito will infect asusceptible humanTvh probability that a bite by a susceptible mosquitoof an infected human will infect the mosquitoτh incubation period in humansτv incubation period in mosquitoesThen we obtain the following modelx′ = bThvm(1− x− y)w−bThvm(1− x(t− τh)− y(t− τh))w(t− τh)e−rτh− rx,y′ = bThvm(1− x(t− τh)− y(t− τh))w(t− τh)e−rτh− ry,z′ = bTvh(1− z−w)y−bTvh(1− z(t− τv)−w(t− τv))y(t− τv)e−µvτv−µvz,w′ = bTvh(1− z(t− τv)−w(t− τv))y(t− τv)e−µvτv−µvw,(4.7)with m = NvNh the number of mosquitoes per human.4.2.2.2 The basic reproduction numberR0 and the endemic equilibriumFor the disease transmission by vector, we need to consider one complete loop”host-vector-host”. If we consider that one single infective individual is intro-87duced into the population, this individual will be bitten and infectbTvh NvNhrsuscep-tible mosquitoes; during the extrinsic incubation period τv, a fraction e−µvτv of themosquitoes survive. Thus, the expected number of secondary cases of infectedmosquitoes by a single infectious human isbTvh NvNhre−µvτv . For the second part ofthe loop, the infectious mosquitoes will bite susceptible humans. For each infec-tious mosquito, the expected number of secondary cases would bebThvµve−rτh . Sothe expected number of secondary cases of humans by introducing a single in-fectious human isb2TvhThv NvNhrµve−µvτv−rτh . Thus, we define the basic reproductionnumberR0:R0 =b2mTvhThvrµve−µvτv−rτh . (4.8)We consider the next generation matrix approach, which is described in [68], tofind the basic reproduction number R0. F is the vector of rates of appearance ofnew infections in disease compartments and V is the vector of rates of transfer ofindividuals into/out of the non-disease compartments by all other means,F =[bThvIve−rτhbTvhNhNvIhe−µvτv],andV =[rIhµvIv].The 2×2 non-negative matrix F and a 2×2 non-singular M-matrix V areF =[0, bThve−rτhbTvhNhNve−µvτv , 0],andV =[r, 00, µv].Then we haveFV−1 =[0 bThve−rτh 1µvbTvhNhNve−µvτv 1r 0].88We calculate the basic reproduction number R0 as the largest eigenvalue of thematrix FV−1,R0 =√b2TvhThv NvNhrµve−µvτv−rτh =√b2TvhThvmrµve−µvτv−rτh . (4.9)However, this approach views the vector transmission as two separate generation, itis reasonable to square it to obtain the expression in (4.8). To find the expressionsof the endemic equilibrium, we need to make the equations in (4.7) equal to 0. Itis easy to find the trivial disease-free equilibrium.We have the following lemma,Lemma 4.2.1. The system (4.7) has at most two equilibria, and we have:(A). IfR0 ≤ 1, the system (4.7) has a unique disease-free equilibrium (0,0,0,0);(B). IfR0 > 1, the system (4.7) has two equilibria, the unique disease-free equi-librium (0,0,0,0) and the endemic equilibrium (x?,y?,z?,w?), withx? =R0−1R0 +bTvhµv e−rτh(1− e−rτh),y? =R0−1R0 +bTvhµv e−rτhe−rτh ,z? =R0−1R0 +bmThvr e−µvτv(1− e−µvτv),w? =R0−1R0 +bmThvr e−µvτve−µvτv .(4.10)From the expressions of the endemic equilibrium, the endemic equilibriumexists if and only ifR0 > 1.Remark 11. The expressions in (4.10) show that at the endemic equilibrium, theratios y?x? =e−rτh1− e−rτh andw?z? =e−µvτv1− e−µvτv are fixed and determined by the incu-bation periods τh and τv, respectively. We will use these relations to analyze thedynamical behaviors at the endemic equilibrium.Remark 12. As mentioned in [60], the basic reproduction number R0 is a de-creasing function of both time delays. And it is easy to prove the prevalence of89disease x?, y?, z? and w? are also decreasing functions of both time delays. It iseasy to verify this both analytically and numerically.The model (4.7) is a four-dimensional dynamical system of delay differentialequations. For the behaviour at the disease-free equilibrium, since the expressionof the equilibrium is simple, it is convenient to use the traditional method to an-alyze it. We will linearize the system of DDEs at the equilibrium and obtain thecharacteristic equation, then analyze the distribution of roots of this characteristicequation. For the behaviour at the endemic equilibrium, through the observationof the asymptotic behavior of solutions near the endemic equilibrium, we are ableto transform the four-dimensional delay differential equation system into a two-dimensional ordinary differential equation system. The obtained system of ODEsis easier to study.4.2.2.3 Behavior at the disease-free equilibriumWe first discuss the stability of the disease-free equilibrium (0,0,0,0). The lin-earized system of (4.7) at (0,0,0,0) isx′ = bThvmw−bThvmw(t− τh)e−rτh− rx,y′ = bThvmw(t− τh)e−rτh− ry,z′ = bTvhy−bTvhy(t− τv)e−µvτv−µvz,w′ = bTvhy(t− τv)e−µvτv−µvw.(4.11)The characteristic equation isdetλ + r, 0, 0, −bThvm(1− e−λτh−rτh)0, λ + r, 0, −bThvme−λτh−rτh0, −bTvh(1− e−λτv−µvτv), λ +µv, 00, −bTvhe−λτv−µvτv , 0, λ +µv= 0,(4.12)or(λ + r)(λ +µv)[(λ + r)(λ +µv)−b2mThvTvhe−rτh−µvτve−(τh+τv)λ ] = 0. (4.13)90Apparently, two negative roots of the characteristic equation (4.13) are −r < 0,−µv < 0, the other two roots satisfy the following equation(λ + r)(λ +µv)−b2mThvTvhe−rτh−µvτve−(τh+τv)λ = 0. (4.14)Equation (4.14) has identical form with equation (5) in [60], and this equation hasbeen studied completely in [60]. The analytical results of (4.14) is summarized as(see [60] for the complete analysis):• If R0 < 1, all roots of (4.14) have negative real parts for all positive τh andτv;• if R0 = 1, equation (4.14) has a zero root and a root with negative real partfor all positive τh and τv;• if R0 > 1, equation (4.14) has a unique positive real root for all positive τhand τv, and has a root with negative real part for small τh and τv.By combining two negative real roots −r and −µv, the distribution of roots ofcharacteristic equation (4.13) is:• If R0 < 1, all roots of (4.13) have negative real parts for all positive τh andτv. We have a stable equilibrium;• ifR0 = 1, equation (4.13) has a zero root and other roots with negative realpart for all positive τh and τv. The disease-free equilibrium is a degenerateequilibrium of codimension one and is stable except in one direction;• if R0 > 1, equation (4.13) has a unique positive real root for all positive τhand τv, and has other roots with negative real part for small τh and τv. Thus,the disease-free equilibrium is unstable, and it is a saddle.4.2.2.4 Behavior at the endemic equilibriumWe now analyze the dynamical behavior at the endemic equilibrium. It is observedthat for the model (4.7), the sum of the first two equations is(x+ y)′ = bThvm(1− x− y)w− r(x+ y),= bThvm(1− (x+ y))w− r(x+ y).91The sum of the last two equations is(z+w)′ = bTvh(1− z−w)y−µv(z+w),= bTvh(1− (z+w))y−µv(z+w).We denote thateh(t) := x(t)+ y(t), ev(t) := z(t)+w(t). (4.15)The biological meaning of eh(t) is the fraction of humans who are infected (infec-tious or not infectious) at time t, ev(t) represents the fraction of mosquitoes whichare infected (infectious or not infectious) at time t. By observing the expressionsof endemic equilibrium in (4.10), the asymptotic behavior of w(t) near w? is likee−µvτvev(t), and the asymptotic behavior of y(t) near y? is like e−rτheh(t), whichare:y(t)≈ e−rτheh(t), w(t)≈ e−µvτvev(t)with t→∞. Thus, we have an approximated ordinary differential equation system:e′h = bThvm(1− eh)e−µvτvev− reh,e′v = bTvh(1− ev)e−rτheh−µvev,(4.16)with the endemic equilibrium:e?h =R0−1R0 +bTvhµv e−rτh,e?v =R0−1R0 +bmThvr e−µvτv.(4.17)Without analyzing the complicated system of DDEs (4.10), we can study the sys-tem of ODEs (4.16) instead.We linearize the system (4.16) at the endemic equilibrium (e?h,e?v) to have[e′he′v]=[−bmThveve−µvτv− r, bmThv(1− eh)e−µvτvbTvh(1− ev)e−rτh , −bTvhehe−rτh−µv][ehev]. (4.18)92The trace of the matrix in (4.18) is −bmThve?ve−µvτv − r− bTvhe?he−rτh − µv < 0 if0 < e?v < 1 and 0 < e?h < 1. The determinant of the matrix in (4.18) is(bmThve?ve−µvτv + r)(bTvhe?he−rτh +µv)−b2TvhThvm(1− e?h)(1− e?v)e−µvτv−rτh=: α+β + γ,(4.19)withα = bmThve−µvτvµve?v +b2mThvTvhe−µvτv−rτhe?v=(bmThve−µvτvµv +b2mThvTvhe−µvτv−rτh) R0−1R0 +bmThvr e−µvτv=(bmThve−µvτvµv +b2mThvTvhe−µvτv−rτh) (R0−1)rµv(R0 +bmThvr e−µvτv)rµv= b2mThvTvhe−µvτv−rτh− rµv= rµv(R0−1),(4.20)β = bTvhe−rτhre?h +b2mThvTvhe−µvτv−rτhe?h=(bTvhe−rτhr+b2mThvTvhe−µvτv−rτh) R0−1R0 +bTvhµv e−rτh=(bTvhe−rτhr+b2mThvTvhe−µvτv−rτh) (R0−1)rµv(R0 +bTvhµv e−rτh)rµv= b2mThvTvhe−µvτv−rτh− rµv= rµv(R0−1),(4.21)and we have γ =−b2mThvTvhe−µvτv−rτh + rµv =−rµv(R0−1).Thus, the determinant of the matrix is rµv(R0−1)> 0 if and only if R0 > 1.Based on the calculations on the trace and the determinant of the matrix in (4.18),the endemic equilibrium (e?h,e?v) is stable. Any positive starting point (eh,0,ev,0),eventually, converges to the endemic equilibrium (e?h,e?v). At the endemic equilib-rium, e?h consists of fraction of e−rτh infectious humans and fraction of 1− e−rτhexposed humans; e?v has the similar structure. Thus, we argue that the endemicequilibrium (x?,y?,z?,w?) of the original model (4.7) is also stable.93We conclude the analysis by the following theorem,Theorem 4.2.2. If R0 < 1, the model (4.7) has a unique equilibrium (0,0,0,0),and this disease-free equilibrium is stable. IfR0 > 1, the model (4.7) has two equi-libria, the disease-free equilibrium (0,0,0,0) is unstable and the endemic equilib-rium (x?,y?,z?,w?) is asymptotically stable.Remark 13. If we consider three special cases:(1). The extrinsic incubation period is taken into account while the intrinsic in-cubation period is not, which is τh = 0 and τv > 0. This model has beendiscussed in several publications, including [45]. The re-scaled system is athree dimensional system of DDEs.(2). The intrinsic incubation period is taken into account while the extrinsic in-cubation period is not, which is τv = 0 and τh > 0. This model is less inter-esting.(3). Both τh = τv = 0, which gives the well-known Ross-Macdonald Model,y′ = bThvm(1− y)w− ry,w′ = bTvh(1−w)y−µvw, (4.22)with y(t) the fraction of infective humans and w(t) the fraction of infectivemosquitoes. This model (4.22) is identical with the model (4.2), only withdifferent notations. There are no time lags of incubation periods. The basicreproduction numberR0 of model (4.22) isR0 =b2mTvhThvrµv,it is also straightforward to find the endemic equilibrium of (4.22), which is(y?,w?) =( R0−1R0 +bTvhµv,R0−1R0 +bmThvr).By linearizing the system (4.22) at the disease-free equilibrium and the en-demic equilibrium, it is easy to show the stabilities.94Our result holds for all these three special cases, it can be verified by pickingspecial value sets of τh and τv. Some comments on these three special cases can byfound in [51] By Martcheva et al..4.2.2.5 Numerical simulationsTo verify our theoretical results in previous sections, we conduct several numericalexperiments. We try to compare our results with the results from [60] by Ruan etal., so we choose the identical sets of parameter values:(A1). r = 0.05, recovery rate in human varies from 0.01 to 0.05 per day;(A2). b = 0.2, biting rate on humans by a single mosquito is from 0.2 to 0.5 perday;(A3). Thv = Tvh = 0.5;(A4). m = 2, this value can be varied due to different environments;(A5). µv = 0.05 and µv = 0.105, the range of death rate of mosquitoes is 0.05 to0.5 per day;(A6). τh = 15 days, the incubation period in human can be as long as severalmonths;(A7). τv = 9 days.More details of sources of the parameters can be found in [60] by Ruan et al..By picking two sets of values of µv, we have two different values of the basicreproduction numberR0. The calculations give us:(B1). if µv = 0.05,R0 = 2.4096. The trivial disease-free equilibrium is (0,0,0,0),and the endemic equilibrium is (0.2217,0.1985,0.1030,0.1812),(B2). if µv = 0.105,R0 = 0.5376. The trivial disease-free equilibrium is (0,0,0,0).We simulate the process of malaria transmission by MATLAB, and have twosets of initial values:(x0,y0,z0,w0) = (0.001,0.001,0.001,0.001), t ≤ 0,95Figure 4.2: The endemic equilibrium is asymptotically stable(x0,y0,z0,w0) = (0.2,0.2,0.2,0.2), t ≤ 0.The numerical results are summarized as:(i) For µv = 0.05 and the initial values (x0,y0,z0,w0)= (0.001,0.001,0.001,0.001).We haveR0 > 1, the endemic equilibrium is asymptotically stable, which isshown in Figure 4.2.(ii) For µv = 0.05 and the initial values (x0,y0,z0,w0) = (0.2,0.2,0.2,0.2). Wehave R0 > 1, the endemic equilibrium is asymptotically stable, which isshown in Figure 4.3.(iii) For µv = 0.105 and the initial values (x0,y0,z0,w0)= (0.001,0.001,0.001,0.001).We haveR0 < 1, the disease-free equilibrium is asymptotically stable, whichis shown in Figure 4.4.(iv) For µv = 0.105 and the initial values (x0,y0,z0,w0) = (0.2,0.2,0.2,0.2). We96Figure 4.3: The endemic equilibrium is asymptotically stablehave R0 < 1, the disease-free equilibrium is asymptotically stable, which isshown in Figure 4.5.(v) It is easy to numerically show that, for R0 > 1, the final proportions of ex-posed humans, infected humans, exposed mosquitoes and infected mosquitoesare decreasing functions of both delays.In paper [60], Ruan et al. mentioned that the value of τv, incubation periodin mosquitoes, is more sensitive to the prevalence of malaria. It can be verifiedby picking different parameters. The stability at the endemic equilibrium is notaffected by the lengths of the delay lags, as long as the basic reproduction numberR0 > 1.97Figure 4.4: The disease-free equilibrium is asymptotically stable4.3 Dynamical analysis of a malaria vaccination model4.3.1 Some facts of the malaria vaccine RTS, SThe RTS, S malaria vaccine candidate (Mosquirix) is currently the most advancedin development globally. This vaccine candidate is developed through a partnerbetween GlaxoSmithKline Biologicals (GSK) and the PATH Malaria Vaccine Ini-tiative (MVI) (from WHO). RTS,S is a vaccine against Plasmodium falciparum.A recent report from Malaria Vaccine Initiative showed that the vaccine candidatecan help to protect children and infants from infection for at least 3 years after firstvaccination (see [49] for more details of the test results of the vaccine).98Figure 4.5: The disease-free equilibrium is asymptotically stable4.3.2 A malaria vaccination modelThe newly developed malaria vaccine has low effectivity. This kind of vaccine iscalled ”disease-modifying vaccine”, which permits infection, lowers the durationof infection and reduces viral load. The vaccine may also have effects, such as:increase the recovery rate, increase the acquired immunity rate and reduce the rateof infection. In this model we divide the initial population into four subgroups,which are: people who never receive the vaccine; people who receive the vaccinebut do not take the vaccine; people who take the vaccine but the vaccine wanesover time, and the vaccine becomes ineffective before the outbreak of the disease;people who receive the vaccine and take it, and the vaccine does not wane overtime. The first three subgroups are unvaccinated groups and the last subgroup is thevaccinated group. For individuals in the vaccinated group, they may have a reducedrate of infection, increased life expectancy and faster recovery. So the durationof infection for vaccinated individuals may increase because of fewer deaths, or99decrease because of higher recovery rates (see [64] for reference). By studyingthe paper [64] by Robert Smith?, we now use the compartmental deterministicmodel to simulate the process of spread. In this vaccination model, the wholehuman population is divided into vaccinated group and unvaccinated group. Foreach group, there are susceptible class, infective class and recovery class. Themosquito population is divided into susceptible class and infected class.We list some parameters of this vaccination model,(A1), Subscripts m, u, v represent population of mosquito, population of unvacci-nated individuals and population of vaccinated individuals, respectively.(A2), For the individuals who enter the population, there is a fraction of p of themwill take the vaccine. Among the vaccine accepters, there is a fraction ofε of people who actually take the vaccine. Thus, there is a fraction of ε pof individuals coming to the vaccinated class, and a fraction of 1− ε p ofindividuals coming to the unvaccinated class.(A3), The vaccine wanes over time at rate ω .(A4), Natural death rates are µm and µh for mosquitoes and humans, respectively.(A5), γu and γv are death rates due to malaria for unvaccinated individual and vac-cinated individual, respectively.(A6), The unvaccinated individual recovers from infection at rate hu without tem-porary immunity; while vaccinated individual recovers from infection with-out temporary immunity at rate hv.(A7), The unvaccinated individual recovers from infection with temporary immu-nity at rate αu; while vaccinated individual recovers from infection with tem-porary immunity at rate αv.(A8), The individual from recovery class will lose immunity at rate δu and δv, forunvaccinated individuals and for vaccinated individuals, respectively.(A9), The infection factor is σ2, which means the probability for a vaccinated in-dividual to be infected is decreased by σ2, and 0≤ σ2 ≤ 1.100(A10), It is possible that the probability that one mosquito biting an unvaccinated in-dividual then getting infected is greater than the probability that this mosquitobiting a vaccinated individual then getting infected. We introduce σ1 with0≤ σ1 ≤ 1. When σ1 = 0, it means that the vaccine will protect susceptiblemosquito from being infected by a vaccinated human. This is the indirectprotection for humans; while if σ1 = 1, the vaccine has no effect of this type.(A11), The vaccine related parameters are: σ1 and σ2, death rates γu and γv whichdue to the disease, and other recovery rates. We will later discuss howthese vaccine related parameters will affect the dynamical behaviours of themodel.The malaria vaccine model isS′m = Λm−bTmhNhSm(Iu +σ1Iv)−µmSm,I′m =bTmhNhSm(Iu +σ1Iv)−µmIm,S′u = (1− ε p)Λh−bThmNhSuIm−µhSu +ωSv +huIu +δuRu,S′v = ε pΛh−σ2bThmNhSvIm−µhSv−ωSv +hvIv +δvRv,I′u =bThmNhSuIm−µhIu− γuIu−huIu−αuIu +ωIv,I′v = σ2bThmNhSvIm−µhIv− γvIv−hvIv−αvIv−ωIv,R′u = αuIu−δuRu−µhRu +ωRv,R′v = αvIv−δvRv−µhRv−ωRv.(4.23)The total population sizes areNh = Su +Sv + Iu + Iv +Ru +Rv,andNv = Sm + Im,101Figure 4.6: Malaria vaccination modelwithdNhdt= Λh−µhNh− γuIu− γvIv,anddNvdt= Λm−µmNv.It is noted that in the model (4.23), the contact rate depends on the sizes ofmosquito population and human population. During the transmission, the contactrate changes in terms of time. This may brings the backward bifurcation behaviouroccurring at the control reproduction numberRc.4.3.3 The control reproduction numberRcWe first consider the case when a single infectious human is introduced into thepopulation, how many secondary cases of infections will this introduction pro-duce. Because the population is not wholly susceptible, part of the populationis vaccinated, it is more appropriate to use the terminology ”control reproductionnumber”. We will interpret the control reproduction number Rc in two different102ways, mathematically and biologically. In the absence of infection, we haveIm = Iu = Iv = Ru = Rv = 0.Then the sizes of mosquito population and human population satisfy the relationsΛm = µmSm,(1− ε p)Λh = µhSu−ωSv,ε pΛh = µhSv +ωSv.(4.24)Solving these equations we are able to obtain the values of S?m, S?u and S?v , whichareS?m =Λmµm,S?u =[µh(1− ε p)+ω]Λh(µh +ω)µh,S?v =ε pΛhµh +ω.We notice thatS?uS?v=µh(1− ε p)+ωε pµh.If the vaccine wane is not taken into account, thenS?uS?v=1− ε pε p.It is shown from these expressions that the rate of wane will affect the initial frac-tion of humans who take the vaccine. If the rate of wane increases, there will belower fraction of vaccinated individuals in human population initially.Further, we find the initial sizes of human population and mosquito population,respectively.N?h = S?u +S?v =[µh(1− ε p)+ω]Λh(µh +ω)µh+ε pΛhµh +ω=Λhµh, (4.25)andN?v = S?m =Λmµm. (4.26)103Before the outbreak of the disease, the sizes of populations of human and mosquitoare fixed. There are certain fraction of humans being vaccinated. We now introducea single infection into this steady state (at disease-free equilibrium).4.3.3.1 The next generation matrix approachWe first consider the next generation matrix approach, which is described in [68],to find the basic reproduction numberR0. We define the disease compartments{Im, Iu, Iv}. (4.27)F is the vector of rates of appearance of new infections in disease compart-ments and V is the vector of rates of transfer of individuals into/out of the non-disease compartments by all other means,F =bTmhNhSm(Iu +σ1Iv)bThmNhSuImσ2 bThmNh SvIm ,andV = µmImµhIu + γuIu +huIu +αuIu−ωIvµhIv + γvIv +hvIv +αvIv +ωIv ,we define piu := µh + γu + hu +αu and piv := µh + γv + hv +αv to simplify the ex-pression, so we haveV = µmImpiuIu−ωIvpivIv +ωIv .The 3×3 non-negative matrix F and a 3×3 non-singular M-matrix V areF =0, bTmhN?h S?m, σ1bTmhN?hS?mbThmN?hS?u, 0, 0σ2 bThmN?h S?v , 0, 0 ,104andV =µm, 0, 00, piu, −ω0, 0, piv +ω ,which gives usV−1 =1µm , 0, 00, 1piu ,ωpiu(piv+ω)0, 0, 1piv+ω ,Then we haveFV−1 =0, bTmhN?hpiu S?m,bTmhN?hS?mωpiu(piv+ω) +σ1bTmhN?hS?m1piv+ωbThmN?hS?u1µm , 0, 0σ2 bThmN?h S?v1µm , 0, 0 .The control reproduction number is determined by the largest eigenvalue of FV−1,which isRc =√bTmhN?hpiuS?mbThmN?hS?u1µm+σ2bThmN?hS?v1µm(bTmhN?hS?mωpiu(piv +ω)+σ1bTmhN?hS?m1piv +ω).Or we can write the control reproduction number asRc =bTmhN?hpiuS?mbThmN?hS?u1µm+σ2bThmN?hS?v1µm(bTmhN?hS?mωpiu(piv +ω)+σ1bTmhN?hS?m1piv +ω)=S?uN?hb2TmhThmN?hS?m1piuµm+S?vN?hσ2bThm1µm(bTmhN?hS?mωpiu(piv +ω)+σ1bTmhN?hS?m1piv +ω).(4.28)It is noticed that if there is no wane of the vaccine, the control reproduction numbercan be expressed asRc =S?uN?hbTmhN?hpiuS?mbThmµm+S?vN?hσ1σ2bThmµmbTmhN?hS?m1piv,by setting ω = 0.1054.3.3.2 Another method of approachWe explain the expression (4.28) from the perspective of biology. We consider aninfected mosquito is introduced into the population, there is a probability of S?uN?hthat individuals bitten by this infected mosquito are unvaccinated, and there is aprobability of S?vN?hthat individuals bitten by this infected mosquito are vaccinated.We now discuss it in two different conditions• Unvaccinated individuals are bitten and become infected, and the infectedindividuals will infect more susceptible mosquitoes. The number of expectedcases of secondary infection is:S?uN?hbTmhN?hpiuS?mbThmµm.• Vaccinated individuals are bitten by the infected mosquito, the vaccinated in-dividuals can become infected and then infect more susceptible mosquitoessince the vaccine is not perfectly effective. During the path ”human tomosquito”, the vaccinated individuals may lose the immunity due to the wan-ing vaccine. The transform between the unvaccinated class and the vacci-nated class can be expressed by the following system of ordinary differentialequations:v′ =−ωv−pivv,u′ =+ωv−piuu,(4.29)where ω is the rate of wane of the vaccine, piu and piv are previously defined.We assume the initial values are v(0) = 1 and u(0) = 0. The analytic solutionof the equation (4.29) isv(t) = e−(piv+ω)t ,u(t) =ωpiu−piv−ω [e−(ω+piv)t − e−piut ]. (4.30)The average times of staying in infective class, for individuals who alwaysstay vaccinated and for individuals who transform from vaccinated to unvac-106cinated, are ∫ ∞0e−(piv+ω)tdt=1piv +ω,(4.31)and ∫ ∞0ωpiu−piv−ω [e−(ω+piv)t − e−piut ]dt=wpiu(piv +ω).(4.32)We summarize the discussions in these two different cases, we can obtain the iden-tical result as the expression in (4.28).4.3.4 The stability of the disease-free equilibriumWe find and analyze the disease-free equilibrium. The disease-related classes inhumans or mosquitoes are either infectious or recovered. We setIm = Iu = Iv = Ru = Rv = 0.We haveS?m =Λmµm,S?u =[µh(1− ε p)+ω]Λh(µh +ω)µh,S?v =ε pΛhµh +ω.The size of human population at the disease-free equilibrium isN?h = S?u +S?v =[µh(1− ε p)+ω]Λh(µh +ω)µh+ε pΛhµh +ω=Λhµh. (4.33)Thus, the disease-free equilibrium of the model (4.23) is{S?m, I?m,S?u,S?v , I?u , I?v ,R?u,R?v}= {Λmµm,0,[µh(1− ε p)+ω]Λh(µh +ω)µh,ε pΛhµh +ω,0,0,0,0}(4.34)107We now linearize the system (4.23) at the disease-free equilibrium, the Jacobianmatrix isJ =−µm, 0, 0, 0, −bTmhN?h S?m, −σ1 bTmhN?h S?m, 0, 00, −µm, 0, 0, bTmhN?h S?m, σ1bTmhN?hS?m, 0, 00, −bThmN?h S?u, −µh, ω, hu 0, σu, 00, −σ2 bThmN?h S?v , 0, −µh−ω, 0, hv, 0, σv0, bThmN?h S?u, 0, 0, −piu, ω, 0, 00, σ2 bThmN?h S?v , 0, 0, 0, −piv−ω, 0, 00, 0, 0, 0, αu, 0, −σu−µh, ω0, 0, 0, 0, 0, αu, 0, −σv−µh−ω,with piu := µh + γu +hu +αu and piv := µh + γv +hv +αv.The eigenvalues of matrix J are−µm,−µh,−µh−ω ,−σu−µh,−σv−µh−ω ,and other three ones are also the eigenvalues of the matrix J1,J1 =−µm, bTmhN?h S?m, σ1bTmhN?hS?mbThmN?hS?u, −piu, ωσ2 bThmN?h S?v , 0, −piv−ω .In appendix B, we prove that the eigenvalues of matrix J1 are all negative byusing the Routh-Hurwitz conditions. Furthermore, all eigenvalues of matrix J arenegative.We summarize the results by the following theorem,Theorem 4.3.1. The malaria vaccination model (4.23) has the disease-free equi-librium (4.34). And this disease-free equilibrium (4.34) is asymptotically stable ifthe control reproduction numberRc < 1.4.3.5 The threshold of eradicationWe have proved previously that if the control reproduction number Rc < 1, thedisease-free equilibrium is asymptotically stable. In this case, the disease of malariacan not spread among the population. Thus, to eradicate the infectious disease, thenecessary condition is to decrease the control reproduction numberRc below 1.108Remark 14. The condition ofRc < 1 is sufficient to control the disease and erad-icate it, only in the case that there is no backward bifurcation existing at Rc = 1.If atRc = 1, there is a backward bifurcation, even whenRc < 1, a stable endemicequilibrium may exist, depending on the initial conditions of the population.If the vaccine coverage level is p = 0, we have the basic reproduction number,since the population is wholly susceptible,R0 =b2TmhThmN?h1piuµm.If the vaccine coverage level is p = 1, which means everyone in the population isvaccinated, we use the control reproduction numberRv, andRv has the expression,Rv = σ2bThm1µm(bTmhN?hS?mωpiu(piv +ω)+σ1bTmhN?hS?m1piv +ω). (4.35)It has been argued that, if the population is partially vaccinated, the control repro-duction number is a linear combination of R0 and Rv. We are able to verify thatthis result still holds for the case of vector-borne-disease, which is in (4.28),Rc =S?uN?hR0 +S?vN?hRv=S?uN?hb2TmhThmN?h1piuµm+S?vN?hσ2bThm1µm(bTmhN?hS?mωpiu(piv +ω)+σ1bTmhN?hS?m1piv +ω).(A). We notice that 1µh+γu+hu+αu and1µh+γv+hv+αv are the durations of the infec-tious periods for individuals from unvaccinated group and from vaccinatedgroup, respectively.(B). For individuals in vaccinated group, the values of hv and αv are decreasedbecause of the protection from vaccine, while the value of γv is increased. Itis hoped that the sum of µh + γv +hv +αv is smaller than µh + γu +hu +αu.But if the primary effect of the vaccine is to keep infectious individual to livelonger time but not die, it is possible the opposite situation happens ([64]).(C). Besides the difference in durations of infectious periods between R0 andRv, another change is that Rv is decreased by the term σ1σ2. These are two109functions of the vaccine, which are protections, for susceptible individualsand susceptible mosquitoes, against infection of malaria.We denote the initial vaccine coverage levelIP =S?vN?hfor simplicity.To reach the herd immunity threshold, it is necessary to make Rc < 1. Toachieve this goal, the weakest requirement is that Rv < 1, which means if thepopulation is fully vaccinated, the disease can not spread and will die out. IfRv >1, even if every individual in the population is vaccinated, the epidemic still canspread.IfRv < 1, we further have(1− IP) ·R0 + IP ·Rv < 1IP >R0−1R0−Rv .(4.36)Since IP = ε p µhµh+ω . It is better to improve the quality of vaccine by decreasing thewaning rate of the vaccine.We conclude the analysis by the following theorem,Theorem 4.3.2. For the malaria vaccination model (4.23), if the initial vaccinecoverage levelp >R0−1R0−Rv ,andRv < 1, the epidemic of malaria will die out.4.4 Conclusion of this chapterIn this chapter, we have studied the vector-borne diseases, especially malaria. Inthe first part, we stated some facts about the vector-borne diseases. In the secondpart, we discussed the disease of malaria and some previous mathematical mod-elling work on the transmission of malaria. In the third part, we formulated a110delayed malaria transmission model, which includes both extrinsic incubation pe-riod and intrinsic incubation period. We calculated the basic reproduction numberR0, the disease-free equilibrium and the possible endemic equilibrium (exists ifand only if R0 is larger than the unity). Then we analyzed the equilibria of thismodel. We also proved that ifR0 < 1, the disease-free equilibrium is stable; whileif R0 > 1, the endemic equilibrium is stable for any lengths of the delay periods.In the last part of this chapter, we formulated and analyzed a malaria vaccinationmodel. The control reproduction numberRc and the expression of the disease-freeequilibrium were calculated. We provided two different methods to calculate thecontrol reproduction number Rc and proved that control reproduction number forpartial vaccinated population can be written as the linear combination of the ba-sic reproduction number R0 and the control reproduction number Rv for whollyvaccinated population. This theory was initially proposed for the analysis of HIVvaccines. We later proved that the disease-free equilibrium is asymptotically sta-ble if the control reproduction number Rc < 1. Finally, we provided the sufficientcondition for the eradication of malaria in a population. It is difficult to calculatethe endemic equilibrium for this vaccination model, because there might exist twoendemic equilibria (one is asymptotically stable, and the other is unstable) whenRc < 1. And there is one stable endemic equilibrium when Rc > 1. Under thiscondition, we will have a backward bifurcation behaviour atRc = 1.111Chapter 5Future work5.1 Epidemiological game theoryThe analysis of vaccination games in heterogeneous mixing population providedus the expressions of the expected vaccine coverage levels in each sub-groups. Itcan be shown that the expected vaccine coverage levels are determined by the con-tact patterns of the two sub-groups. A paper ”An SIS Epidemiology Game withTwo Sub-populations” written by Timothy Reluga ([56]) analyzed the rational-expectation equilibria in an epidemiological game with two interacting sub-populationsof equal size. The transmission dynamics are described by the SIS epidemic model.Individuals are only allowed to invest in daily prevention measures. The results ofthis paper showed that disassortative mixing may lead to multiple Nash equilibria,and the interior Nash equilibria are always unstable. If the mixing is positivelyassortative, there is a unique globally stable Nash equilibrium.It is worthwhile to investigate deeper for the case of heterogeneous mixingpopulation (two sub-groups for example) which we discussed in chapter 2. Inchapter 2, we have the expressions of the expected vaccine coverage levels for thecase of perfect vaccine,p = 1− Cr1β22−Cr2β21β11β22−β12β211pip1, q = 1− Cr2β11−Cr1β12β11β22−β12β211piq2,where Cr1 and Cr2 are relative costs in two groups. Because the values of fractions112are Cr1β22−Cr2β21β11β22−β12β21 andCr2β11−Cr1β12β11β22−β12β21 , they can be any numbers, even be negative val-ues. However, through the analysis in chapter 3, since the function (1− p)pip1 ismonotonically decreasing, we are able to conclude that there is at most one Nashequilibrium. Otherwise, we will have a ”dominant” game. For the case of imperfectvaccine, there may exist more than one Nash equilibrium. We need to theoreticallyprove the uniqueness for these two games. The patterns of mixing in the populationare the keys.5.2 Attack ratios and vaccine coverage levelsIn chapter 3, we formulated the deterministic compartmental vaccination models,in homogeneous mixing population and in heterogeneous population, respectively.The final size relations can be obtained from the models. By analyzing the final sizerelations, we proved that the attack ratios can be expressed in terms of the vaccinecoverage levels, and also proved that the attack ratios are decreasing functions ofthe vaccine coverage levels. We provided the mathematical proof that the vaccinecan protect both susceptible population and vaccinated population better, this is theso called direct and indirect protection stemming from vaccine. The probabilityof every individual to survive through the epidemic is greater by increasing thevaccine coverage level. Some possible future work related to this part of work are:(1). Compare the measures of derivatives to conduct the sensitivity analysis. Ifone derivative is bigger than the other, from the standpoint of public healthmanagement, to vaccinate one certain subgroup is the priority, since sameamount of vaccination can lower the attack ratio more.(2). The analysis in case of heterogeneous mixing population only derives ana-lytical results for the vaccine efficacy factor which is close to 1 and closeto 0. Due to many parameters of the model, there are no explicit formulasto be given. By other methods, we may obtain more accurate results, onepossible way is to consider proportional mixing or preferred mixing, insteadof general mixing.(3). We prove that in heterogeneous mixing population, if two subgroups aremainly separate, the vaccine efficacy for subgroup 1 does not affect subgroup1132, vice versa.5.3 Contact patterns in heterogeneous mixing populationContact patterns are very important in analyzing transmissions of epidemics. Thereare a lot of related papers, such as [53] and [41]. To advance the research in epi-demiological games in heterogeneous mixing population, the base is to understandthe theory of contact/mixing patterns.5.4 Some work in malaria transmission modelsI propose two possible extensions for malaria transmission models.(1). It is necessary to include the birth rate and death rate in humans. From ournumerical experiments, if the initial fractions of infected humans or infectedmosquitoes are small, the basic reproduction number R0 > 1, it will take500−600 days to reach a stable state. The birth/death rates in humans shouldbe considered for such a long period. By including the disease death, we willobtain a general model of malaria transmission with incubation periods. Oneof the difficulties in analyzing such a model is that the re-scaled system is ahigh-dimensional system of delay differential equations.(2). We consider the heterogeneity in mosquitoes. The mosquito population isdivided into two age groups, the first group is with age from 0 to T − τv,the second group is with age from T − τv to T , where τv is the extrinsicincubation period in mosquitoes, and T is the possible maximum age formosquitoes. For the first group, a susceptible mosquito will be infected bybiting infective humans, then stays in the exposed class for an incubationperiod of τv before entering the infective class. For the second group, asusceptible mosquitoes will also be infected by contacting with infective hu-mans, but the difference is it will stay in the exposed class and then die. Weassume that the age density function of mosquito is D(x), with 0 ≤ x ≤ Trepresenting the age. This density function D(x) is static, does ont depend114on the natural time. We have ∫ T0D(x)dx = 1. (5.1)The fraction of mosquitoes with age between 0 and T − τv isφ =∫ T−τv0D(x)dx. (5.2)If an infective human is bitten by susceptible mosquitoes, there is a proba-bility of φ that this individual is bitten by mosquitoes from the first group.The fraction of 1−φ of mosquitoes will die before becoming infectious. Itis reasonable to have the basic reproduction number be φR0, where R0 hasthe formula (4.8).The assumption of static age density function is oversimplified. We considerthe age density function of mosquitoes depends on age and time and thefunctionρ(a, t) (5.3)denotes the density of individuals of age a at time t. The total population isdefined asP(t) =∫ T0ρ(a, t)da. (5.4)The fractions can be calculated. It has been proved that the density functionρ(a, t) satisfy the McKendrick equation:ρa(a, t)+ρt(a, t)+µaρ(a, t) = 0. (5.5)This equation is also called the von Foerster equation, and function µa ≥ 0is called the mortality function or death modulus. More details can be foundin chapter 7 of [13] by Brauer and Castillo-Chavez. Some empirical data onthe age structure of mosquitoes is needed to perform numerical experiments.1155.5 Backward bifurcation analysis for malariavaccination modelsThe research of malaria vaccination model in this thesis only examined the condi-tion of eradication, which relates to the control reproduction number Rc and thedisease-free equilibrium. It is interesting and important to perform backward bi-furcation analysis for this model. If there is backward bifurcation for some setsof parameters, even the control reproduction number Rc < 1, there is an endemicequilibrium for certain initial conditions. This certainly brings difficulties in publichealth management, sinceRc < 1 is not enough to eradicate the epidemic.We continue our analysis based on the malaria vaccination model we proposedin previous chapter. At the endemic equilibrium, the following equation must besatisfied,Λh−µhNeh = γuIeu− γvIev , (5.6)with Neh the size of the population at the endemic equilibrium, which consists ofboth vaccinated classes and unvaccinated classes. Ieu is the number of unvaccinatedand infective individuals, Iev is the number of vaccinated but infective individuals.It is apparently that the value of Neh depends on the values of Ieu and Iev . If we assumethat the disease-induced death rates γu and γv are both 0, we haveΛh = µhNeh .As reported by WHO, there were 214 million cases of malaria and 438000 deathsin 2015. In a closed population, the disease-induced death rate can be neglected.However, by considering the disease-induced death rates, the stability of endemicequilibrium may change dramatically. For example, a paper [24] by Chitnis, Cush-ing and Hyman proved that, for a mathematical malaria model with exposed period,in the absence of disease-induced death, the transcritical bifurcation R0 = 1 is su-percritical (forward); and numerically proved that for larger values of the disease-induced death rate, a subcritical (backward) bifurcation is possible atR0 = 1.To start our analysis on the possible impact of vaccination, we analyze the mostsimple case as the starting point. We set several parameters to be 0 in the previousmalaria vaccination model and propose the Ross-Macdonald malaria vaccination116model,S′m = Λm−bTmhNhSm(Iu +σ1Iv)−µmSm,I′m =bTmhNhSm(Iu +σ1Iv)−µmIm,S′u = (1− ε p)Λh−bThmNhSuIm−µhSu +huIu,S′v = ε pΛh−σ2bThmNhSvIm−µhSv +hvIv,I′u =bThmNhSuIm−µhIu−huIu,I′v = σ2bThmNhSvIm−µhIv−hvIv.(5.7)To reduce the dimension of this system, we define the proportions of popula-tions in mosquitoes and in humans,(1). sm =SemNm, im =IemNm.(2). su =SeuNh, iu =IeuNh.(3). sv =SevNh, iv =IevNh.Thus, we obtain the following system of equations,i′m = bTmhsm(iu +σ1iv)−µmim,s′u = (1− ε p)µh−bmThmsuim−µhsu +huiu,s′v = ε pµh−σ2bmThmsvim−µhsv +hviv,i′u = bmThmsuim−µhiu−huiu,i′v = σ2bmThmsvim−µhiv−huiv.(5.8)In this simple system, there is no transmission between unvaccinated individu-als and vaccinated individuals. At the endemic equilibrium, we havesu + iu = 1− ε p,sv + iv = ε p.The second and third equations in the system (5.8) can be removed. By setting all117equations in (5.8) be zero, we can find the expressions of the endemic equilibrium.For the purpose of simplicity, we still use im, iu and iv to represent the fractions atthe endemic equilibrium.Equations (4) and (5) in (5.8) give ussuσ2sv=(hu +µh)iu(hv +µh)iv,1− ε− iuε− iv =σ2(hu +µh)iu(hv +µh)iv.Reorganize this equation, we haveiu =(1− ε p)ivP(ε p− iv)+ iv , (5.9)withK =σ2(hu +µh)hv +µh.From equation (1) in (5.8), we haveim =bTmh(iu +σ1iv)bTmh(iu +σ1iv)+µm. (5.10)We now analyze equation (5) from the system (5.8) to obtainσ2b2mThmTmhε p(1− ε pK(ε p− iv)+ iv +σ1)= [σ2b2mThmTmh +bTmh(hv +µh)](1− ε pK(ε p− iv)+ iv +σ1)iv +µm(hv +µh),σ2b2mThmTmhε p(1− ε p+σ1Kε p+(1−K)σ1iv)= [σ2b2mThmTmh +bTmh(hv +µh)](1− ε p+σ1Kε p+(1−K)σ1iv)iv+ µm(hv +µh)[Kε p+(1−K)iv].Thus, we obtain a quadratic equation for iv,Ai2v +Biv +C = 0, (5.11)118withA := [σ2b2mThmTmh +bTmh(hv +µh)](1−K)σ1,B := [σ2b2mThmTmh +bTmh(hv +µh)](1− ε p+σ1Kε p)−σ2b2mThmTmhε p(1−K)σ1 +µm(hv +µh)(1−K),C := µm(hv +µh)Kε p−σ2b2mThmTmhε p(1− ε p+σ1Kε p).It is noticed that in this model, the control reproduction numberRc has the expres-sion,Rc = (1− ε p) mb2ThmTmh(µh +hu)µm+σ1σ2ε pmb2ThmTmh(µh +hv)µm.The open problem is whether there are two solutions to equation (5.11) whenRc < 1. If there are two solutions under this condition, whether these two solutionsare bounded by 0 and 1. If the answers are yes, we have a backward bifurcation atRc = 1.After solving this problem, it is straightforward to consider a more generalmalaria vaccination with immunity. We consider that the unvaccinated individualsrecovers from infection with temporary immunity at rate αu > 0 and the vaccinatedindividuals recovers from infection with temporary immunity at rate αv > 0. Wealso assume that the individuals from recovery class will lose immunity at rate δuand δv, for unvaccinated class and for vaccinated class, respectively. Thus, a more119general malaria vaccination model isS′m = Λm−bTmhNhSm(Iu +σ1Iv)−µmSm,I′m =bTmhNhSm(Iu +σ1Iv)−µmIm,S′u = (1− ε p)Λh−bThmNhSuIm−µhSu +δuRu,S′v = ε pΛh−σ2bThmNhSvIm−µhSv +δvRv,I′u =bThmNhSuIm−µhIu−αuIu,I′v = σ2bThmNhSvIm−µhIv−αvIv,R′u = αuIu−δuRu−µhRu,R′v = αvIv−δvRv−µhRv.(5.12)We still set1. sm =SemNm, im =IemNm,2. su =SeuNh, iu =IeuNh, ru =ReuNh,3. sv =SevNh, iv =IevNh, rv =RevNh,to obtain the re-scaled model,i′m = bTmhsm(iu +σ1iv)−µmim,s′u = (1− ε p)µh−bmThmsuim−µhsu +δuru,s′v = ε pµh−σ2bmThmsvim−µhsv +δvrv,i′u = bmThmsuim−µhiu−αuiu,i′v = σ2bmThmsvim−µhiv−αviv,r′u = αuiu−δuru−µhru,r′v = αviv−δvrv−µhrv.(5.13)It has been argued that in epidemic models with multiple groups and asym-metry between groups or multiple interaction mechanisms it is possible to have120a very different bifurcation behavior at R = 1 ([21], [29], [37], [38] and [46]).The malaria vaccination models are very likely to exhibit backward bifurcationbehaviours.121Bibliography[1] Anderson R.M. (ed.), 1982, Population Dynamics of Infetious Disease:Thoery and Applications, Chapman and Hall, London. → pages 81, 82, 83[2] Anderson R.M., May R.M., 1985, Vaccination and herd immunity to infec-tious disease, Nature 9: 318-323. → pages 34[3] Anderson R.M., May R.M., 1991, Infectious diseases of humans: dynamicsand control, Oxford University Press, London. → pages 48, 83[4] Arino J., Brauer F., van den Driessche P., Watmough J., Wu J., 2007, Afinal size relatin for epidemic models, Mathematical Biosciences and En-gineering 4: 159-175. → pages[5] Arino J., Brauer F., van den Driessche P., Watmough J., Wu J., 2008, Amodel for influenza with vaccination and antiviral treatment, Journal ofTheoretical Biology 253: 118-130. → pages[6] Aron J.L., 1988 Mathematical modeling of immunity to malaria, MathBiosci 90:385-396. → pages 83[7] Bailey N., 1982, The Biomathematics of malaria London, Charles Griffinand Co Ltd, London. → pages 83[8] Barrett S., 2003, Global disease eradication, Journal of the European Eco-nomics Association 1: 591-600. → pages 47[9] Barrett S., 2013, Economic considerations for the eradication endgame,Phil Trans R Soc B 368: 20120149. → pages 47[10] Bauch C., Earn D., 2004, Vaccination and the theory of games, Proc. NatlAcad. Sci. USA 101: 13391-13394. → pages 52122[11] Bauch C., Galvani A.P., Earn D.J.D., 2003, Group interest versus self-interest in smallpox vaccination policy, PNAS 100: 10564-10567.→ pages52[12] Brauer F., Castillo-Chavez C., 2010, Mathematical Models in PopulationBiology and Epidemiology, Springer. → pages 15, 35, 36, 115[13] Brauer F., 2005, The KermackMcKendrick epidemic model revisited, Math-ematical Biosciences and Engineering 198: 119-131. → pages 15, 35, 36,115[14] Brauer F., 2006, Some simple epidemic models, Mathematical Biosciencesand Engineering 3: 1-15. → pages 35, 36[15] Brauer F., 2008, Age-of-infection and the final size relation, MathematicalBiosciences and Engineering 5: 681-690. → pages[16] Brauer F., 2008, Epidemic Models with Heterogeneous Mixing and Treat-ment, Bulletin of Mathematical Biology 70: 1869-1885. → pages 15, 25[17] Brauer F., 2008, Age of infection and the Final Size Relation, MathematicalBiosciences and Engineering 5: 681-690. → pages[18] Brauer F., 2008. Epidemic Models with Heterogeneous Mixing and Treat-ment, Bulletin of Mathematical Biology 70: 1869-1885. → pages 38[19] Brauer F., Watmough J., 2009, Age of infection epidemic models with het-erogeneous mixing, Journal of Biological Dynamics 3: 324-330. → pages15, 38, 47, 71[20] Brauer F., van den Driseeche P., Wu J. (Eds.), 2008, mathmatical epidemi-ology, Springer. → pages 1, 2, 9[21] Brauer F., 2004, Backward bifurcation in simple vaccination models, J.Math. Anal. Appl. 298: 418-431. → pages 121[22] Brito D.L., Sheshinski E., Intriligator M.D., 1991, Externalities and com-pulsory vaccinations, Journal of Public Economics 45: 69-91. → pages52[23] Centers for Disease Control and Prevention, 2015, Fast Facts aboutMalaria, http://www.cdc.gov/malaria/about/facts.html. → pages 79123[24] Chitnis N., Cushing J.M., Hyman J.M., 2006, Bifurcation analysis of amathematical model for malaria transmission, SIAM J. APPL. MATH. 67:24-45. → pages 116[25] Dangerfield C.E., Ross J.V., Keeling M.J., 2009, Integrating stochasticityand network structure into an epidemic model, J R Soc Interface 6:761-774. → pages 83[26] Diekmann O., Heesterbeek J.A.P., 2000, Mathematical Epidemiology ofInfectious Diseases, Wiley, Chichester. → pages 12[27] Diekmann O., Heesterbeek J.A.P., Metz J.A.J., 1990, On the definition andthe computation of the basic reproductive ratioR0 in models for infectiousdiseases in heterogeneous populations, J. Math. Biol. 28: 365-382. →pages 12[28] Dietz K., Molineaux L., Thomas A., 1974, A malaria model tested in theAfrican savannah, Bull World Health Organ 50: 347-357. → pages 83[29] Dushoff J., Huang W., Castillo-Chavez C., 1998, Backwards bifurcationsand catastrophe in simple models of total diseases, J. Math. Biol. 36:227-248. → pages 121[30] Eshel I., 1996, On the changing concept of population stability as a reflec-tion of a changing problematics in the quantitative theory of evolution, J.Math. Biol. 34:485-510. → pages 21[31] Farr W., 1984, Second annual report of the Registrar-General of Births,Deaths and Marriages of England and Wales. → pages 33[32] Fine P., 1993, Herd Immunity: History, Theory, Practice, EpidemiologicReview 15: No.2. → pages 33[33] Fine P. E. M., Clarkson J. A., 1986, Individual versus public prioritiesin the determination of optimal vaccination policies, American Journal ofEpidemiology 124: 1012-1020. → pages 23, 51[34] Fox J., Elveback L., Scott W., Gatewood L., Ackerman E., 1971, Herd im-munity: Basic concept and relevance to public health immunization prac-tices, Am J Epidemiology 94: 179-189. → pages 32[35] Galvani A. P., Regula T. C., Chapman G. B., 2007, Long-standing influenzavaccination policy is in accord with individual self-interest but not with theutilitarian optimum, PNAS 104: 5692-5697. → pages 52124[36] Gubler D. J., 2009, Vector-borne diseases, Rev. sci. tech. Off. int. Epiz.28(2): 583-588. → pages ix, 77, 78, 79, 80[37] Hadeler K.P., Castillo-Chavez C., 1995, A core group model for diseasetransmission, Math. Biosci. 128: 41-55. → pages 121[38] Hadeler K.P., van den Driessche P., 1997, Backward bifurcation in epi-demic control, Math. Biosci. 146: 15-35. → pages 121[39] Hamer WH., 1906, Epidemic disease in England - the evidence of variabil-ity and of persistency of type, Lancet. → pages 33[40] Hay S.I., Myers M.F., Burke D.S., Vaughn D.W., Endy T., Ananda N.,Shanks G.D., Snow R.W., Rogers D.J., 2000, Etiology of interepidemicperiods of mosquito-borne disease, PNAS 97: 93359339. → pages 83[41] Hethcote H. W., Van Ark J. W., 1987, Epidemiological models for het-erogeneous populations: proportionate mixing, parameter estimation, andimmunization programs, Mathematical Biosciences 84: 85-118. → pages114[42] Kermack W.O., McKendrick A.G., 1927, A contribution to the mathemati-cal theory of epidemics, Proc. Royal Soc. London 115: 700-721. → pages9[43] Kermack W.O., McKendrick A.G., 1932, A contribution to the mathemati-cal theory of epidemics, part. II, Proc. Royal Soc. London 138: 55-83. →pages 9[44] Kermack W.O., McKendrick A.G., 1932, A contribution to the mathemat-ical theory of epidemics, part. III, Proc. Royal Soc. London 141: 94-112.→ pages 9[45] Koella J C., 1991, On the use of mathematical models of malaria transmis-sion, Acta Tropica 49: 1-25. → pages 83, 85, 94[46] Kribs-Zaleta C.M., Velasco-Hernndez J.X., 2000, A simple vaccinationmodel with multiple endemic states, Math. Biosci. 164: 183-201. → pages121[47] Ma J., Earn D., 2006, Generality of the Final Size Formula for an Epidemicof a Newly Invading Infectious Disease, Bulletin of Mathematical Biology68: 679-702. → pages 73125[48] Macdonald, G., 1957, The Epidemilogy and Control of Malaria, OxfordUniversity Press, London. → pages 34, 80[49] Malaria Vaccine Initiative, 2015, RTS, S,http://www.malariavaccine.org/malaria-and-vaccines/first-generation-vaccine/rtss. → pages 98[50] Mandal S., Sarkar R.R., Sinha S., 2011, Mathematical models of malaria -a review, Malaria Journal, 10: 202. → pages 81, 83[51] Martcheva M., Prosper O., 2012, Unstable Dynamics of Vector-Borne Dis-eases: Modeling Through Delay Differential Equations, Dynamic Modelsof Infectious Diseases: Vol 1. Vector Borne Diseases Springer, 43-75. →pages 78, 85, 95[52] Maynard S. J. and Price G.R., 1973, The logic of animal conflicts, Nature246:15-18. → pages 6, 7[53] Nold Annett, 1980, Heterogeneity in Disease-Transmission Modeling,Mathematical Biosciences 52: 227-240. → pages 114[54] Nowak M. A. and Sigmund K., 2004, Evolutionary dynamics of biologicalgames, Science 303:793799. → pages 8[55] Parham P.E., Michael E., 2011, Outbreak properties of epidemic models:The roles of temporal forcing and stochasticity on pathogen invasion dy-namics, Journal of Theoretical Biology 271:1-9. → pages 83[56] Reluga T.C., 2008, An SIS epidemiology game with two subpopulations, JBiol Dyn. 3: 515-531. → pages 52, 112[57] Reluga T.C., Bauch C.T., Galvani A.P., 2006, Evolving public perceptionsand stability in vaccine uptake, Math Biosci. 204(2):185-98. → pages 18,52[58] Rogers D.J., Randolph S.E., Snow R.W., Hay S.I., 2002. Satellite imageryin the study and forecast of malaria, Nature 415: 710715. → pages 83[59] Ross R., 1911, The prevention of Malaria, 2nd ed., (with Addendum), JohnMurray, London. → pages 80[60] Ruan S., Xiao D., Beier J C., 2008, On the delayed Ross-Macdonald Modelfor Malaria Transmission, Bulletin of Mathematical Biology 70: 1098-1114. → pages 84, 85, 86, 89, 91, 95, 97126[61] Sachs J., Malaney P., 2002, The economic and social burden of malaria,Nature 415:680-685. → pages 83[62] Shim E., Chapman G.B., Townsend J.P., Galvani A.P., 2012, The influenceof altruism on influenza on influenza vaccination decisions, J. R. Soc. In-terface 0: 1-10. → pages 52[63] Smith D.L., McKenzie F.E., 2004, Statics and dynamics of malaria infec-tion in Anopheles mosquitoes, Malaria Journal 3:13. → pages 86[64] Smith? R., 2007, Could Low-Efficacy Malaria Vaccines Increase Sec-ondary Infections in Endemic Areas?, (In: Mathematical Modeling of Bio-logical Systems, Volume II. → pages 100, 109[65] Tassier T., 2013, The Economics of Epidemiology, Springer. → pages 20[66] Topley W., Wilson GS., 1923, The spread of bacterial infection. The prob-lem of herd immunity, Journal of Hygiene 21: 21-243. → pages 33[67] Torres-Sorando L., Rodriguez D.J., 1997, Models of spatio-temporal dy-namics in malaria, Ecol Model 104:231-240. → pages 83[68] van den Driessche P., Watmough J., 2002, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease trans-mission, Math. Biosci. 180: 29-48. → pages 12, 39, 60, 65, 88, 104[69] von Neumann J. and Morgenstern O., 1944, Theory of Games and Eco-nomic Behavior, Princeton University Press, Princeton. → pages 6127Appendix ASimilarity of forms of final sizerelationsWe let S(t) denote the number of susceptibles at time t and let φ(t) be the totalinfectivity at time t, defined as the sum of products of the number of infectedmembers with each infection age and the mean infectivity for that infection age.We assume that on average members of the population make a constant number aof contacts in unit time. We let B(τ) be the fraction of infected members remaininginfected at infection age τ and let pi(τ) with 0≤ pi(τ)≤ 1 be the mean infectivityat infection age τ . Then we haveA(τ) = B(τ)pi(τ),the mean infectivity of members of the population with infection age τ . We assumeno disease deaths, so the total population is a constant, which is N0, the vaccinecoverage level is p, and assume the vaccine is perfectly effective, thus we have theinitial number of susceptibles is (1− p)N0.The age-of-infection model isS′ =− aN0Sφ (A.1)128φ(t) = φ0(t)+∫ t0aN0S(t− τ)φ(t− τ)A(τ)dτ= φ0(t)+∫ t0[−S′(t− τ)]A(τ)dτ.(A.2)The control reproduction number isRc = (1− p)a∫ ∞0A(τ)dτ, (A.3)we can use the expression of Rc to estimate how large the vaccine coverage level pis to reach the herd immunity.Now we are looking for the final size relation, we write the model asS′S=− aN0φ=− aN0(φ0(t)+∫ t0[−S′(t− τ)]A(τ)dτ),then we integrate with respect to t from 0 to ∞ to obtainlnS0S∞=aN0∫ ∞0φ0(t)dt +aN0∫ ∞0∫ t0[−S′(t− τ)]A(τ)dτ)dt=aN0∫ ∞0φ0(t)dt +aN0∫ ∞0A(τ)∫ ∞τ[−S′(t− τ)]dtdτ=aN0∫ ∞0φ0(t)dt +a[S0−S∞N0]∫ ∞0A(τ)dτ,since in the assumptions, we have that S0 = (1− p)N0 and other part of thepopulation has fully immunity, and φ0(t) is the total infectivity of members of thepopulation who were infected at t = 0 at time t, so φ0(t) = 0 for any t ≥ 0. Thus,the final size relation can be expressed aslnS0S∞= a[S0−S∞N0]∫ ∞0A(τ)dτ. (A.4)It is clearly that this final size relation has similar formula with the final sizerelation of the classic SIR model.We now consider the case that the vaccine is imperfect, which is modeled byincluding a factor σ , with 0≤ σ ≤ 1, in the infection rate of vaccinated members,129with σ = 0 meaning that the vaccine is perfectly effective and σ = 1 meaning thatthe vaccine has no effect. It is described by the following age-of-infection modelwith imperfect vaccineS′ =− aN0Sφ (A.5)V ′ =−σ aN0Vφ (A.6)φ(t) = φ0(t)+∫ t0aN0S(t− τ)φ(t− τ)A(τ)dτ+∫ t0σaN0V (t− τ)φ(t− τ)A(τ)dτ= φ0(t)+∫ t0[−S′(t− τ)]A(τ)dτ+∫ t0[−V ′(t− τ)]A(τ)dτ.(A.7)To obtain the final size relation we haveS′S=− aN0φ=− aN0(φ0(t)+∫ t0[−S′(t− τ)]A(τ)dτ+∫ t0[−V ′(t− τ)]A(τ)dτ),andV ′V=−σ aN0φ=−σ aN0(φ0(t)+∫ t0[−S′(t− τ)]A(τ)dτ+∫ t0[−V ′(t− τ)]A(τ)dτ),then we integrate with respect to t from 0 to ∞ to obtainlnS0S∞=aN0∫ ∞0φ0(t)dt +aN0∫ ∞0∫ t0[−S′(t− τ)]A(τ)dτ)dt + aN0∫ ∞0∫ t0[−V ′(t− τ)]A(τ)dτ)dt=aN0∫ ∞0φ0(t)dt +aN0∫ ∞0A(τ)∫ ∞τ[−S′(t− τ)]dtdτ+ aN0∫ ∞0A(τ)∫ ∞τ[−V ′(t− τ)]dtdτ=aN0∫ ∞0φ0(t)dt +a[S0−S∞N0]∫ ∞0A(τ)dτ+a[V0−V∞N0]∫ ∞0A(τ)dτ,130andlnV0V∞= σaN0∫ ∞0φ0(t)dt +σaN0∫ ∞0∫ t0[−S′(t− τ)]A(τ)dτ)dt +σ aN0∫ ∞0∫ t0[−V ′(t− τ)]A(τ)dτ)dt= σaN0∫ ∞0φ0(t)dt +σaN0∫ ∞0A(τ)∫ ∞τ[−S′(t− τ)]dtdτ+σ aN0∫ ∞0A(τ)∫ ∞τ[−V ′(t− τ)]dtdτ= σaN0∫ ∞0φ0(t)dt +σa[S0−S∞N0]∫ ∞0A(τ)dτ+σa[V0−V∞N0]∫ ∞0A(τ)dτ,it is easy to find out thatV0V∞= (S0S∞)σ . (A.8)The final size relation is consistent with SVIR model, the analysis on dependenceof initial values will still hold.We now consider a population with a distribution of subpopulation identifiedby the variable ζ ranging over a ’state’ space Ω, having sizes N(t,ζ ), respectively,at time t. Suppose that each group ζ member makes a(ζ ) contacts sufficient totransmit infection in unit time. We assume that the fraction of contacts made by amember of group ζ that is with a member of group η is p(ζ ,η), then:∫Ωp(ζ ,η)dη = 1,for each ζ . The total number of contacts made in unit time by member of group ζwith members of group η isa(ζ )p(ζ ,η)N(t,ζ ),and by balance relation we havea(ζ )p(ζ ,η)N(t,ζ ) = a(η)p(η ,ζ )N(t,η),with ζ 6= η .And we assume that there is no permanent movement between groups and thatthere are no disease deaths, so that N(t,ζ ) is a constant function N(ζ ) of t for allζ .131A multi-group age of infection model with general mixing would beS′(t,ζ ) =−a(ζ )S(t,ζ )∫Ωp(ζ ,η)φ(t,η)N(η)dη (A.9)φ(t,ζ ) = φ0(t,ζ )+∫ ∞0[−S′(t− τ,ζ )]A(τ,ζ )dτ. (A.10)The initial conditions are S0(ζ ) = (1− p(ζ ))N(ζ ) and the vaccine is perfectlyeffective.It is possible to obtain a system of final size relations. Integration of equationfor S′(t,ζ )S(t,ζ ) in Equation (A.9) giveslnS(0,ζ )S(∞,ζ )=∫ ∞0a(ζ )∫Ωp(ζ ,η)φ(t,η)N(η)dηdt=∫ ∞0a(ζ )∫Ωp(ζ ,η)φ0(t,η)+∫ ∞0 [−S′(t− τ,η)]A(τ,η)dτN(η)dηdt= a(ζ )∫Ωp(ζ ,η)N(η)∫ ∞0(φ0(t,η)+∫ ∞0[−S′(t− τ,η)]A(τ,η)dτ)dtdη= a(ζ )∫Ωp(ζ ,η)N(η)(∫ ∞0φ0(t,η)+∫ ∞0∫ ∞τ[−S′(t− τ,η)]A(τ,η)dtdτ)dτ= a(ζ )∫Ωp(ζ ,η)N(η)(∫ ∞0φ0(t,η)dt +(S(0,η)−S(∞,η))∫ ∞0A(η)dt)dη .If we have only two sub-population, it is straightforward to replace the form of in-tegration by the form of summation. When we perform the analysis of dependenceon initial values (the initial vaccine coverage levels), we will obtain similar results.We now investigate the case which the vaccine are imperfect, and we modelthis by including factors σ(ζ ), where 0≤ σ(ζ )≤ 1. With similar set ups, we havethe multi-group age of infection model with general mixingS′(t,ζ ) =−a(ζ )S(t,ζ )∫Ωp(ζ ,η)φ(t,η)N(η)dη (A.11)V ′(t,ζ ) =−σ(ζ )a(ζ )V (t,ζ )∫Ωp(ζ ,η)φ(t,η)N(η)dη (A.12)132φ(t,ζ ) = φ0(t,ζ )+∫ ∞0[−S′(t− τ,ζ −V ′(t− τ,ζ )]A(τ,ζ )dτ. (A.13)It is easy to use the same way to obtain the relations of final sizes. Integrationsof equations for S′(t,ζ )S(t,ζ ) andV ′(t,ζ )V (t,ζ ) in Equations (A.11) and (A.12) will givelnS(0,ζ )S(∞,ζ )=∫ ∞0a(ζ )∫Ωp(ζ ,η)φ(t,η)N(η)dηdt=∫ ∞0a(ζ )∫Ωp(ζ ,η)φ0(t,η)+∫ ∞0 [−S′(t− τ,η)−V ′(t− τ,η)]A(τ,η)dτN(η)dηdt= a(ζ )∫Ωp(ζ ,η)N(η)∫ ∞0(φ0(t,η)+∫ ∞0[−S′(t− τ,η)−V ′(t− τ,η)]A(τ,η)dτ)dtdη= a(ζ )∫Ωp(ζ ,η)N(η)(∫ ∞0φ0(t,η)+∫ ∞0∫ ∞τ[−S′(t− τ,η)−V ′(t− τ,η)]A(τ,η)dtdτ)dτ= a(ζ )∫Ωp(ζ ,η)N(η)(∫ ∞0φ0(t,η)dt +(N(η)−S(∞,η)−V (∞,η))∫ ∞0A(η)dt)dη ,andlnV (0,ζ )V (∞,ζ )=∫ ∞0σ(ζ )a(ζ )∫Ωp(ζ ,η)φ(t,η)N(η)dηdt=∫ ∞0σ(ζ )a(ζ )∫Ωp(ζ ,η)φ0(t,η)+∫ ∞0 [−S′(t− τ,η)−V ′(t− τ,η)]A(τ,η)dτN(η)dηdt= σ(ζ )a(ζ )∫Ωp(ζ ,η)N(η)∫ ∞0(φ0(t,η)+∫ ∞0[−S′(t− τ,η)−V ′(t− τ,η)]A(τ,η)dτ)dtdη= σ(ζ )a(ζ )∫Ωp(ζ ,η)N(η)(∫ ∞0φ0(t,η)+∫ ∞0∫ ∞τ[−S′(t− τ,η)−V ′(t− τ,η)]A(τ,η)dtdτ)dτ= σ(ζ )a(ζ )∫Ωp(ζ ,η)N(η)(∫ ∞0φ0(t,η)dt +(N(η)−S(∞,η)−V (∞,η))∫ ∞0A(η)dt)dη .133Appendix BAsymptotic stability of thedisease-free equilibriumTo simplify the calculations, we use the following matrix to represent the matrixJ1, a1, a2, a3a4, a5, a6a7, 0, a9 .The entries of this matrix represent the entries having same positions in matrix J1.The characteristic equation corresponding to the matrix has the formula,λ 3 +b1λ 2 +b2λ +b3 = 0.To prove all eigenvalues of this equation are negative, we need to show:1, b1 > 0,2, b3 > 0,3, b1b2 > b3,under the conditionRc < 1.134Proof. In the characteristic equation, we haveb1 =−(a1 +a5 +a9),b2 = a1a9 +a5a9 +a1a5−a2a4−a3a7,b3 = a3a5a7 +a2a4a9−a1a5a9−a2a6a7.Now we prove the three items in order,(1), b1 =−(a1 +a5 +a9) = um +piu +piv +ω > 0;(2),b3 = a3a5a7 +a2a4a9−a1a5a9−a2a6a7= a5(a3a7−a1a9)+a2(a4a9−a6a7)=−piu(σ1σ2b2ThmTmhN?hS?vN?hS?m−µm(piv +ω))+bTmhN?hS?m(− bThmN?hS?u(piv +ω)−σ2bThmN?hS?vω)= piuµm(piv +ω)(1−Rc)> 0,if and only if the control reproduction numberRc < 1.(3),b1b2−b3 =− (a1 +a5 +a9)(a1a9 +a5a9 +a1a5−a2a4−a3a7)− (a3a5a7 +a2a4a9−a1a5a9−a2a6a7)=−a1a1a9−a1a5a9−a1a1a5−a5a5a9−a1a5a5−a1a9a9−a5a9a9−a1a5a9+a1a2a4 +a1a3a7 +a2a4a5 +a3a7a9 +a2a6a7.We observe that only the terms a1a2a4, a1a3a7, a2a4a5 and a3a7a9 are nega-tive. We have:−a1a1a5 +a1a2a4 = µ2mpiu(1−b2TmhThmN?hS?uN?hS?m1µmpiu)> µ2mpiu(1−Rc)> 0,135−a1a5a5 +a2a4a5 = µmpi2u (1−b2TmhThmN?hS?uN?hS?m1µmpiu)> µmpi2u (1−Rc)> 0,−a1a1a9 +a1a3a7 = µ2m(piv +ω)(1−σ1σ2b2TmhThmN?hS?vN?hS?m1µm(piv +ω))> µ2m(piv +ω)(1−Rc)> 0,−a1a9a9 +a3a7a9 = µm(piv +ω)2(1−σ1σ2 b2TmhThmN?hS?vN?hS?m1µm(piv +ω))> µm(piv +ω)2(1−Rc)> 0.All other terms in b1b2−b3 are positive, so we have b1b2 > b3.We now have the conclusion that all eigenvalues of matrix J1 are negative if thecontrol reproduction numberRc < 1.136`

Cite

Citation Scheme:

Citations by CSL (citeproc-js)

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
```                        ```
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
`https://iiif.library.ubc.ca/presentation/dsp.24.1-0305652/manifest`