Modeling and Stability Analysis for Measles Metapopulation Model with Vaccination
Leopard C. Mpande1, Damian Kajunguri1, Emmanuel A. Mpolya2
1Nelson Mandela African Institution of Science and Technology, School of CoCSE, Arusha, Tanzania
2Nelson Mandela African Institution of Science and Technology, School of LiSBE, Arusha, Tanzania
To cite this article:
Leopard C. Mpande, Damian Kajunguri, Emmanuel A. Mpolya. Modeling and Stability Analysis for Measles Metapopulation Model with Vaccination.Applied and Computational Mathematics.Vol.4, No. 6, 2015, pp. 431-444. doi: 10.11648/j.acm.20150406.16
Abstract: In this paper, a metapopulation model is formulated as a system of ordinary differential equations to study the impact of vaccination on the spread of measles. The disease-free equilibrium is computed and proved to be locally and globally asymptotically stable if and unstable if . We show that when there are no movements between the two patches, there exists at least one endemic equilibrium for all and bifurcation analysis of endemic equilibrium point proves that forward (supercritical) bifurcation occurs in each patch. Numerical simulation results are also presented to validate analytical results and to show the impact of vaccination on the incidence and prevalence of measles in a metapopulation.
Keywords: Vaccination, Metapopulation, Measles, Bifurcation Analysis
Measles is a contagious disease and is due to infection of Paramyxovirus of the genus Morbillivirus [34,35]. An incubation period for measles is found somewhere between 9 and 12 days and its infectivity period is found between 4 and 9 days . Globally, the disease is said to be one of the most prominent causes of death among young children, despite the presence of an effective vaccine . Measles is easily transmitted by coughing and sneezing, especially when someone stays in direct contact with an infected nasal secretions . It has been emphasized that in the year 2013 there were 145,700 measles induced deaths globally, which is equivalent to 400 deaths every day or 16 deaths every hour .
Measles cases occur if there is no high coverage of vaccination . The high number of cases occurs in places where there is an aggregation of individuals who have not been vaccinated or infected by the disease . Measles has a basic reproduction number of the range 6 to 45 which implies that the mean number of secondary infections caused by a single infected individual in a susceptible population is found somewhere between 6 and 45 .
Several studies have been conducted on the use of mathematical models to control infectious diseases such as measles [1,2,27,29,34,39]. These studies respectively, studied the effect of vaccination  and area  on transmission dynamics of measles, estimated basic reproduction number for measles , studied control of measles by vaccination incorporating two phases of infectiousness , used bifurcation theory on the mathematical model to study measles dynamics , and predicted an optimal vaccine coverage level needed to control measles [18,39]. There are also other studies which use metapopulation models to control infectious diseases such as measles [4,5,14,37,43]. These models play an important role in studying disease epidemics because they can describe the dynamics of individuals between patches which may be cities, towns, and so forth. These studies respectively, presented a system of 4p ordinary differential equations to describe disease spread in an environment divided into p patches and extended their system to include cross infection between several patches and keeping track of both the current patch and the patch in which an individual usually resides [4,5], presented a fractional SEIR metapopulation system modeling the spread of measles by considering 4 distinct patches which are cities , proposed a metapopulation model for regional measles dynamics on the basis of a gravity coupling model and a time series susceptible-infected-recovered (TSIR) model for local dynamics , formulated a disease transmission model as a system of ordinary differential equations for a population with individuals traveling between discrete geographic patches .
In this study, we propose a metapopulation mathematical model as a system of ordinary differential equations to study the impact of vaccination on the spread of measles. Our metapopulation model consists of two regions one with high measles infection (patch 1) and the other region with a low measles infection (patch 2) and movement of individuals between the patches in all direction at constant rates is considered.
2. Model Formulation
In this section, we formulate a measles metapopulation model incorporating vaccination as a control strategy. Our model consists of two patches, where each patch is divided into the following epidemiological classes (for : Susceptible, Vaccinated , Exposed , Infected , and Recovered. We assume that individuals mix homogeneously. Recruitment is assumed to be through birth at constant rate.
Natural mortality rate is constant for all patches. We assume one dose of vaccination for susceptible individuals at a rate . Once an individual is vaccinated, he or she goes to recovered class with permanent immunity at a constant rate . The average number of effective contacts of an infectious individual per unit time is, and standard incidence is assumed. The exposed individuals move from exposed class to infectious class at a rate. The infectious individuals recover permanently after treatment at the rate . Our metapopulation model represents two regions, patch 1 with high measles infection and patch 2 with a low measles infection with an assumption of individual movements between patches in both directions at equal rates as shown in figure 1. The forces of infections for each patch are given by and respectively.
Per capita birth rate in patch .
Contact rate (the average number of adequate contacts per person per unit time) in patch .
The rate of progression from latent class to infectious class in patch .
Vaccine coverage rate in patch .
Recovery rate of treated infectious individuals in patch .
Per capita natural mortality rate in patch .
Disease induced death rate in patch .
Recovery rate of vaccinated individuals in patch .
From the description of the dynamics of measles and with the aid of the compartmental diagram in Figure 1, we have the following set of differential equations.
with initial conditions, , , , and for [4, 5, 37].
Here, is the total population in each patch and satisfies .
The total population size in all patches is .
The following two lemmas show that the model is well posed and that all variables lie in the interval where .
Lemma 1: The solution for the model system 1 is positively invariant in the positive orthant .
Consider the first equation. Assume there exist a time such that , and for .
But we have which is a contradiction to the assumption . This implies that remains positive for all. Similarly, it can be shown that for all , the variables and remain positive for all . Hence solutions remain non-negative for nonnegative initial conditions. Therefore the model is considered to be mathematically and epidemiologically well-posed. Basing on biological considerations, model system (1) will be studied in the region
Lemma 2. Consider the system (1) with nonnegative initial conditions. Assume that for all , the variables , , , and remain non-negative, then remain positive, and the total populationis bounded above for .
Proof. Assume non-negative initial conditions.
For all , we have.
Thus for which shows thatprovided. Thusprovided that.
By summing all the equations we have
If at a certain time, , then at , so is non-increasing at . Thus is bounded above by .
The right hand sides of (1) are continuously differentiable, hence basic theorems  can be used to show that there is a unique solution to the system with given non-negative initial conditions and that this solution exists for all . Therefore the model is considered to be mathematically and epidemiologically well-posed.
3. Model Analysis
The model system (1) is analyzed qualitatively to give better understanding of the impact of vaccination on the epidemiology of measles.
3.1. Disease Free Equilibrium (DFE),
The metapopulation model is at equilibrium if the time derivatives are zero. In the case of system 1, the metapopualtion model is at DFE if for all . Thus, at a disease free equilibrium we have .
Solving the system (1), we get a disease-free equilibrium point where
3.2. The Effective Reproduction Number,
Stability of equilibrium can be analyzed using the basic reproduction number [3,11,20]. The basic reproduction number is the expected number of secondary cases produced by a typical infective individual introduced into a completely susceptible population, in the absence of any control measure. A general method for computing is the next generation method [13,41]. Mathematically, is the spectral radius of the so-called next generation matrix. Here, we compute the control reproduction number, denoted by , to describe the average number of secondary cases generated by primary cases under specified controls such as vaccination [3,20]. Using the method described by , we use to denote the rates of appearance of new infections in each compartment; , being the vector of individual transfer rates into the particular compartment, and the vector of individual transfer rates out of the particular compartment. The two vectors are given by
, and .
The next generation matrix is defined as , where and are both the Jacobian matrices of and evaluated at disease free equilibrium with respect to exposed and infectious classes.
After some calculations we found
The next generation matrix , has a nonzero eigenvalue corresponding to the spectral radius which represents the control reproduction number of the model as given in (2).
If , the disease cannot invade the metapopulation and the infection will die out over a period of time, and also, if , then an invasion is possible and infection can spread through the mepopulation. Generally, the larger the value of, the more severe, and possibly widespread the epidemic will be.
where , , ,,, and .
When there is no vaccination in all patches, we set the parameters and to zero and we get and . Thus we get the basic reproduction number as shown in (3)
where , , , and .
We now consider the case when there are no movements between the given two patches. This means that the parameters and become zero. Hence the control reproduction numbers for patch 1 and patch 2 are given in the form (for )
When there are no vaccination strategies, we set the parameter and equal to zero and hence the reproduction numbers for the two patches when there are no movements between them are given in the form (for )
3.3. Local Stability of the Disease-Free Equilibrium
Here, we investigate the local stability of the disease free equilibrium point, by employing the method described in [17, 23, 33, 39] to linearize the model system (1) by computing its Jacobian matrix . The Jacobian matrix is computed at disease free equilibrium point by differentiating each equation in the system with respect to the state variables and . We get
where , , , , , , , , , , , and .
An equilibrium point is locally asymptotically stable if the Jacobian matrix has a negative trace and a positive determinant or if all of its eigenvalues have negative real parts [15, 18, 23, 30]. Using the idea of [17, 23] we write the jacobian matrix in the form
,, and .
The disease-free equilibrium is locally asymptotically stable if and only if all the eigenvalues of the matrices and have negative real parts. The eigenvalues of are
It can also be shown that all eigenvalues of have negative real parts. Thus, it is clear that for , the DFE is locally asymptotically stable, so that the infection does not persist in the metapopulation and under this condition the endemic equilibrium point does not exist. The DFE is unstable for , and then the endemic equilibrium point exists and the infection persists in the mepopulation. Therefore we established the following Lemma.
Lemma 3. With nonnegative initial conditions the disease-free equilibrium of the system (1) is locally asymptotically stable if and unstable if .
3.4. Global Stability of Disease Free Equilibrium Point (DFE)
In this section, we use the method developed in [11,32,34] to analyze the global stability of disease free equilibrium point. We state two conditions which guarantee the global stability of disease free equilibrium point. The model system (1) can be written in the form
where denotes (its components) the number of uninfected individuals and denotes (its components) the number of infected individuals including latent, infectious, etc.
We use as a disease free equilibrium of this system. According to  the conditions and below must be met to guarantee local asymptotic stability.
: For , is globally asymptotically stable (g.a.s).
: ,for ,
whereis an M-matrix (the off-diagonal elements of are non-negative) and is the region where the model makes biological sense. Considering our model system (1), we have
Which clearly shows that is globally asymptotically stable (g.a.s). Therefore, the condition is satisfied.
For the second condition we have
where , , , and .
Since and , it is clear that .
Now consider the right hand side of
So the condition is also satisfied. Thus, is globally asymptotically stable (g.a.s). Therefore, we have the following important Lemma.
Lemma 4. With non-negative initial conditions, the DFE of the model system (1) is globally asymptotically stable if and unstable if .
3.5. Existence and Local Stability of Endemic Equilibrium (EE) Point,
In the presence of infection the model system (1) has a non-trivial equilibrium point, known as endemic equilibrium point given by
. The endemic equilibrium is an equilibrium where at least one of the components or is nonzero [12, 33]. We compute the endemic equilibrium point by setting the equations of the model system (1) to zero. Since the endemic equilibrium cannot be clearly expressed in closed form, we find the conditions for its existence as done in [24, 40]. We can reduce the model by eliminating,, and to obtain the system
For the existence of an endemic equilibrium the following condition must be satisfied
or or or i.e.
Or or or or or .
Adding equations in the system (6) above at an endemic equilibrium we have
which is equivalent to
Since and we can observe that
and which implies or or or or or .
Therefore endemic equilibrium pointof the model exists.
The reduced model system given in (6) can be studied as means of attacking the model system (1). Thus, we will use this reduced model system for checking global stability of endemic equilibrium in section 3.7. Since the disease free equilibrium is locally asymptotically stable as we have proved in section 3.3, this will imply local stability of endemic equilibrium point for the model system (1). In the next section we are going to investigate the existence and local stability of endemic equilibrium point for patch 1and patch 2 when there are no individual movements between them using bifurcation analysis theory.
3.6. Stability Analysis Using Bifurcation Analysis
Bifurcation analysis plays an important role in disease control and eradication. In this section we study the existence and stability of endemic equilibrium point of the two patches when there exists no individual movements between them and determine the existence of either forward (supercritical) or backward (subcritical) bifurcation. When a forward bifurcation occurs then we guarantee that reducing basic reproduction number to a value less than one is a sufficient condition for disease eradication. On the other hand when a backward bifurcation occurs, an endemic equilibrium may also occur for. This means that must be reduced further so as to avoid endemic states and ensure the eradication . We apply theorem 1 as done in [7, 10, 15, 26] which is based on the use of center manifold theory , to establish local stability of endemic equilibrium point corresponding to patch 1 and patch 2.
Considering patch 1 and patch 2 in isolation, we have the following model system (for )
It can be shown that for existence of endemic equilibrium in patch , the system (7) must satisfy the equation where and .
It follows that
From (8) it can be proved that a positive endemic equilibrium exists in patch if .
The model system (7) has effective reproduction number and a basic reproduction number as defined in (4) and (5) respectively.
For studying the direction of bifurcation we transform the system (7) by setting , , , and .
The model system (7) can be written in the form as follows
We choose as a bifurcation parameter. Solving for when we get
Theorem 1. Consider the general system of ordinary differential equations with a parametersuch that 
Without loss of generality we assume that is an equilibrium point of the system. Thus for all .
1. is Jacobian (linearization) matrix of the system around the equilibrium with evaluated at . Zero is a simple eigenvalue of and all other eigenvalues of A have negative real parts.
2. Matrix has a (nonnegative) right eigenvector and a left eigenvector corresponding to the zero eigenvalue.
Let denote the component of and,.
Then the local dynamics of the system around are totally determined by the sign of and . In particular, if then a backward bifurcation occurs at .
i. . When with x=0 is locally asymptotically stable and there exists a positive unstable equilibrium; when is unstable and there exists a negative and locally asymptotically stable equilibrium.
ii. . When with x=0 is unstable; when x=0 is locally asymptotically stable and there exists a positive unstable equilibrium.
iii. . When with is unstable and there exists a locally asymptotically negative stable equilibrium; when is stable and a positive unstable equilibrium appears.
iv. . Whenchanges from negative to positive, changes its stability from stable to unstable. Correspondently, a negative unstable equilibrium becomes positive and locally asymptotically stable.
Remark. The requirement that is non-negative is unnecessary .
Clearly, at a transcritical bifurcation takes place: more precisely, when such a bifurcation is forward; when the bifurcation is backward.
Now applying theorem 1, the Jacobian matrix of the system (9) at disease free equilibrium evaluated at is given by
The eigenvalues of are , , , and .
Since is a simple eigenvalue of and all other eigenvalues have negative real parts, then assumption 1 of theorem 1 is verified.
The right eigenvector of corresponding to zero eigenvalue is given by
The left eigenvector of satisfying is given by where
Considering system (7) and only nonzero components of the left eigenvector, we compute the values of and at disease free equilibrium as defined in theorem 1 as follows.
The disease free equilibrium in patch is given by
We consider the functions and as defined in (9).
Associated nonzero partial derivatives at the disease fee equilibrium and are given by
, , , ,
, and .
It follows that,
The sign of in (10) depends on the sign of . If then , and if then . Thus, we have established the following theorem.
Figure 2. a and shows the forward bifurcation diagrams for patch 1 and patch 2 respectively obtained from numerical simulations. DFE stands for disease free equilibrium and EE stands for endemic equilibrium. The two diagrams show that the disease free and endemic equilibria exchange stability when for . This means that the disease free equilibrium is locally asymptotically stable when and unstable when . Furthermore, a unique endemic equilibrium exists for and it is locally asymptotically stable. So, the total number of infectious individual in each patch goes to a unique endemic equilibrium.
(i) If then patch exhibit forward bifurcation at . When changes from negative to positive, the disease free equilibrium changes its stability from stable to unstable. Correspondently, a negative unstable endemic equilibrium becomes positive and locally asymptotically stable when.
(ii) If then patch exhibits backward bifurcation at. When with, the disease free equilibrium is locally asymptotically stable and there exists a positive unstable endemic equilibrium; when , the disease free equilibrium is unstable and there exists a negative and locally asymptotically stable endemic equilibrium.
The bifurcation diagrams for patch 1 and patch 2 are shown in figure 2 and respectively.
An implication of EE point being locally asymptotically stable is that the disease can still invade in the metapopulation and transmission dynamics can persist if control measures for the disease are not highly considered in each patch. Therefore, our study agrees that reducing the reproduction number to value less than one is a sufficient condition to eliminate the disease.
3.7. Global Stability of Endemic Equilibrium Point
In this section we analyse the global stability of the endemic equilibrium point by constructing a suitable Lyapunov function. For simplicity, we consider the reduced model system (6) to prove for global stability. We employ the approach of  as it is used for many complicated epidemiological models. We consider the Lyapunov function of the form
where (for is a properly chosen positive constant in the given region. is a population of compartment and is the equilibrium level. So we define the Lyapunov function as
The time derivative of L is
At an endemic equilibrium pointwe have
is non-positive following the modified version of Barbalat’s Lemma  or by following the approach of [25, 32]. Thus, for . Hence and is zero when , , , , , . Therefore, the largest invariant set in Ω such that is the singleton which is our endemic equilibrium point. By LaSalle’s invariant principle  we conclude that is globally asymptotically stable (g.a.s). Thus, we establish the following theory.
Theorem 3. When the endemic equilibrium point is globally asymptotically stable in.
4. Simulation and Discussion
The main objective of this study was to study the impact of vaccination on the spread of measles in a metapopulation. In order to support the analytical results, graphical representations showing the variations in parameters with respect to different state variables have been presented in this section. This is done by using a set of parameter values whose sources are mainly from literature as well as estimation in order to have more realistic simulation results. We will vary key parameters to investigate the impact of vaccination on the transmission dynamics of measles. The parameter values are shown in table 2.
Figure 4. and b shows variations of susceptible, vaccinated, exposed and recovered individuals in patch1and patch 2 respectively when individual movements between them are not allowed.
In figures 3 and 4 above we can see that the susceptible population in both patches decrease rapidly to lower levels with time due to high number of individuals who become vaccinated or exposed due to high contact rates. Exposed population increases more rapidly in patch 1 than in patch 2 due to high contact rates in patch 1. The exposed population later starts to decrease to lower levels due to large number of individuals who become infected or vaccinated. In both patches, the infected population decreases to lower levels with time due to high vaccination and treatment rates. On the other hand, due to treatment and vaccination, recovered population increases to higher levels in both patches as shown in the figures 3 and 4 above.
It can be observed that as vaccination rates increase, the measles prevalence and incidence decrease to lower levels. Thus figures 5 and 6 depict positive impact of vaccination on measles prevalence and incidence in metapopulation. Therefore, our study suggests higher vaccination coverage in all patches in order to eradicate the disease in metapopulation.
In this paper we presented a mathematical model for the control of measles in a metapopulation by considering two regions (patches). We used estimated data and data from literature in numerical simulation. We started by showing nonnegativity of solutions to the metapopulation model, thereby addressing the problem of its well posedness. We proved the disease equilibrium points of the model to be locally and globally asymptotically stable ifand unstable if. We performed bifurcation analysis of endemic equilibrium points of the two patches when there exist no movement of individuals between them and found that forward (supercritical) bifurcation occurs in both cases, which agrees with an intuition that reducing reproduction number to values less than one is a necessary and sufficient condition for disease eradication in the community [8, 41]. Simulation results of different epidemiological classes revealed that most of the individuals undergoing treatment or vaccination join the recovered class. Through simulations, we also showed that vaccination has a positive impact on measles incidence and prevalence in a metapopulation.
The authors would like to thank the Nelson Mandela African Institution of Science and Technology for financial assistance.