Modeling the Dynamics of an Epidemic under Vaccination in Two Interacting Populations

We present a model for an SIR epidemic in a population consisting of two components—locals and migrants. We identify three equilibrium points and we analyse the stability of the disease free equilibrium. Then we apply optimal control theory to find an optimal vaccination strategy for this 2-group population in a very simple form. Finally we support our analysis by numerical simulation using the fourth order Runge-Kutta method.


Introduction
Mathematical modeling of the numerical evolution of infectious diseases has become an important tool for disease control and eradication when possible. Much work has been done on the problem of how a given population is affected by another population when there is mutual interaction. The mere presence of migrant people poses a challenge to whatever health systems are in place in a particular region. Such epidemiological phenomena have been studied extensively, described by mathematical models with suggestions for intervention strategies. The epidemiological effect of migration within the population itself was modeled for sleeping sickness in a paper 1 by Chalvet-Monfray et al. In the case of malaria, there is for instance a study 2 by Tumwiine et al. on the effect of migrating people on a fixed population. The latter two diseases are vector borne. Diseases that propagate without a vector spread perhaps more easily when introduced into a new region. Various studies of models with immigration of infectives have been undertaken for tuberculosis, see for instance 3 by Zhou et al., or the work 4 of Jia et al., and for HIV, see the paper 5 of Naresh et al.
A very simple compartmental model of an epidemic would be an autonomous system comprising a system of two or three differential equations, such as, for instance, the model of

Model Formulation
To study the transmission of a disease in two interacting populations, we consider the total population with size N, as being divided into two subpopulations, the migrant subpopulation of size M, and the local subpopulation of size L. We assume that each subpopulation size is constant the rate of birth equals the mortality rate and that the population is uniform and homogeneously mixing. Divide each subpopulation into disjoint classes called the susceptible class S , the infectious class I , and the class of the removed R . Thus, there will be three such classes for the local population and also three classes for the migrant population. Figure 1: Flow chart of two interacting populations.
The sizes of these classes change with time and will be denoted by S 0 t , I 0 t , R 0 t , S 1 t , I 1 t , and R 1 t . Let us agree henceforth to suppress the subscript 0 for local population, writing simply S t instead of S 0 t , and so on. The model is described by a system of six differential equations as follows. The schematic diagram depicted in Figure 1 illustrates the model and informs the differential equations. We note that the first three equations in 2.1 constitute an SIR model as, for instance, in the paper 10 by Zaman et al. Let us normalize the variables, using the new variables s 1 S 1 /M, i 1 I 1 /M, r 1 R 1 /M, s S/L, i I/L and r R/L. After normalizing our model, which we shall refer to as model 2.1 and 2.2 , becomes as follows: dr t dt γi t − vr t u t s t .

4 Journal of Applied Mathematics
Here v 1 and v are the mortality rate equal to the birth rate in the migrant subpopulation, and the local subpopulation, respectively. The functions u 1 t and u t are the percentages of susceptible individuals being vaccinated in the respective subpopulations per unit time. Individuals enter the recovered compartment at rates γ 1 and γ for the respective subpopulations. Also β 1 and β are the transmission coefficients from the susceptible compartment into the infectious, for the migrant subpopulation and the local subpopulation, respectively. The transmission coefficient from migrants to locals is denoted by β 2 . The term β 2 i 1 s models the influence of the migrant subpopulation onto the locals as in the paper 4 of Jia et al. In the normalized system above, the sizes of the two groups in the population are not visible. At least we should be aware of their relative sizes. In particular, the weighting constant c 0 must be in step with the ratio M/L. The feasible region for the system is the following set:

Equilibria and Their Stability
Equilibrium points are time-independent solutions to the given system of equations. Therefore, in this subsection, we assume u 1 t and u t to be constant functions, u 1 t ≡ u 1 and u t ≡ u. Stability properties of the equilibria are closely linked with the numbers We shall prove that the basic reproduction ratio is the number R u,u 1 max{K, K 1 }. This will follow from Proposition 3.3.

Notation 1.
If u and u 1 are both identically zero, then R u,u 1 will be written as R 0 . For an equilibrium point E, the coordinates will be denoted by E s , E s 1 , and so on.
We take advantage of the aforementioned equivalence in presenting our next theorem.

Theorem 3.2.
Let one consider the unvaccinated version of model 2.1 and 2.2 , that is, with u t ≡ 0 and u 1 t ≡ 0, and let us further assume that v 1 v. If R 0 < 1, then the disease-free equilibrium F with F s 1 and F s 1 1 exists and is globally stable.
Proof. In view of Remark 3.1, this theorem is a direct consequence of 4, Theorem 1 .
Turning to the more general model 2.1 and 2.2 , with vaccination and without the assumption v 1 v, we can identify three possible equilibrium points. Proposition 3.3. a If R u,u 1 < 1, then the disease-free equilibrium F is locally asymptotically stable and its coordinates are 3.3 c The endemic equilibrium D has coordinates as follows:

3.5
Proof. The given points F, D, B ∈ R 6 clearly are equilibrium solutions, which may or may not be feasible. a The Jacobian associated with the system 2.1 and 2.2 at point F is We set out to find the eigenvalues of W. This amounts to solving for λ in the equation, where q 1 λ and q 2 λ are the quadratic expressions below: Now from 3.9 we can write q 1 in the form where A 1 and A 2 are the constants:

3.12
Substituting the equilibrium values at the point F of s 1 , i 1 , s and i, we can rewrite

3.13
The roots of q 1 have negative real parts if both A 1 and A 2 are positive. Now we note that A 2 is positive if and only if 14 that is, Journal of Applied Mathematics 7 If K 1 < 1, then also A 1 > 0. From 3.10 we have q 2 as follows: Now let us define the coefficients Q 1 and Q 2 as

3.17
By applying a similar analysis as for q 1 , we find that the roots of 3.16 have negative real parts if and only if both Q 1 and Q 2 are positive, which is equivalent to the condition Therefore, the disease-free equilibrium is locally asymptotically stable if K 1 < 1 and K < 1, that is, when R u,u 1 < 1. b and c : The points are obtained by direct computation. Feasibility of B is clear when K and K 1 are as given in b .
We include a computational example of an endemic equilibrium point D. We note that P x also has a root x 0.50866, but this is not a feasible value for D s since the corresponding D i value −1.00141 is negative.
In line with the terminology of 4 , we shall refer to the point B as a boundary equilibrium. Stability analysis of the points B and D would take more effort than in the case of F, and could distract from the main purpose of this paper.

Optimal Vaccination Strategy
We wish to design optimal vaccination strategies u * t and u * 1 t , respectively, for the local population and the migrant population. We have six state variables s 1 t , s t , . . . , r t . The variable u t denotes the percentage of susceptible individuals being vaccinated per unit of time in the local population, and u t is assumed to be bounded, 0 ≤ u t ≤ α ≤ 1. A similar 8 Journal of Applied Mathematics interpretation holds for u 1 t , and we assume that for some constant α 1 , 0 ≤ u 1 t ≤ α 1 ≤ 1. Our optimal control problem amounts to minimizing the objective function below where c 0 , c, and c 1 are positive weighting constants. The integral in the objective function can be regarded as follows. The first two terms in the integrand represent the suffering, the lost working hours, the cost of hospitalization, and so on, due to infections. The other two terms represent the cost of vaccination. Similar objective functions are considered in the book 16 of Lenart and Workman and in, for instance, the paper 10 of Zaman et al. Our problem is then as follows.
Problem 1. Minimize J u t , u 1 t subject to the system 2.1 and 2.2 of differential equations, together with the initial conditions and terminal conditions, s 1 T , i 1 T , r 1 T , s T , i T , and r T are free, while the control variables are assumed to be measurable functions that are bounded above The Hamiltonian for this problem is as follows: H t, s 1 , i 1 , r 1 , s, i, r, u, λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 i t c 0 i 1 t cu t 2 c 1 u 1 t 2 In the theorem below, the controls, the state variables, and the costate variables are functions of time. However, notationally this dependence will be suppressed except when required explicitly. The upper dot denotes the time derivative.

Theorem 4.1. An optimal solution for Problem 1 exists. An optimal solution satisfies the identity
Journal of Applied Mathematics 9 and also satisfies the following system of differential equations:

4.6
with transversality conditions: Furthermore, the optimal vaccination strategy is given by

4.8
Proof. Existence of a solution follows since the Hamiltonian is convex with respect to u t and u 1 t . We check the first-order conditions for this optimization problem. We calculate the partial derivatives of the Hamiltonian with respect to the different state variables, in order to obtain the time derivativesλ i t of the costate variables. Due to s 1 T , i 1 T , r 1 T , s T , i T and r T being free, the following terminal conditions hold: We start off by observing that, This implies that λ 3 t and λ 6 t are of the form for some constants A and B, respectively. The terminal conditions λ 3 T 0 and λ 6 T 0, forces A and B to vanish. Therefore, λ 3 and λ 6 are identically zero, that is, λ 3 t ≡ 0 and λ 6 t ≡ 0 as claimed in the theorem. Now we calculatė and we obtain the equations as asserted in the theorem.

Journal of Applied Mathematics
We now turn to the final part of the proof, which is about the form of the controls, u * t and u * 1 t . The function u * t must optimize H. So we calculate Consider a fixed value of t. Now if 2cu t − λ 4 t s t is zero for some value of u t in 0, α , then the given value of u t is optimal. If for every number u ∈ 0, α , we have 2cu − λ 4 t s t ≥ 0 resp., 2cu − λ 4 t s t ≤ 0 , 4.14 then we must choose u t 0 resp., u i α . Thus, we must have The function u * 1 t also must optimize H, and by a similar argument we obtain the stated expression for u * 1 t .

Numerical Simulation
We  We note that if both groups have the infection on a significant scale, then the optimal strategy is to vaccinate in both groups on a comparable scale. The optimal vaccination rollouts for the two groups are similar in form Figures 2, 3, and 4 . Naively, one would expect to see that in such a case the optimal strategy should be to vaccinate the migrants at much higher ratios than the locals. Our simulation reveals that although the initial infection on locals is zero, it is optimal to immediately start on vaccination of the locals whenever there are infected migrants Figure 5 . We conclude this section with a comparison of the values of the objective functional, comparing cases of contant vaccination with optimal and variable vaccination strategies.  Table 1   u 1 0.7 and u 0.7 of u 1 t , and u t , respectively, and in the last column we use the values u 1 0.7 and u 0.3. The computations in the table show that the application of just a constant vaccination strategy would render excessively large values of the objective functional in comparison with the optimal vaccination strategy. The resulting benefit therefore makes it worth the effort of calculating the optimal control.

Concluding Remarks
We observe the influence of the migrant subpopulation onto a given local population, and then determine an optimal vaccination strategy for the two-group population. In the model of Jia et al., the emphasis is on the high impact of migrants. The paper 4 of Jia et al., and also the work 2 of Tumwiine et al., together with some other papers, support the theory that migrants have considerable influence in transmission of most communicable diseases. Our model facilitates numerical illustration of these phenomena. Example 4.3, in particular, shows how optimal control theory informs the best strategy, eliminating the risk of naive decision making. With the results of this paper, we are now able to efficiently plan the rollout of the appropriate vaccination strategy on such a two-group population. Future research may include similar investigations on stochastic two-group SIR models such as the model in 12 of Yu et al. In particular, the approach of Lahrouz et al. 9 to the stochastic model, permits an SIR model in which the total population stays constant, while stochasticity prevails in the propagation of the disease. This method shows much promise in epidemiological modeling.