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

Abstract

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.

1. 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 Kermack and McKendrick. There are more sophisticated models that allow for an incubation period for the pathogen after entering the body of a host. One of the ways of dealing with this phenomenon is by way of delay differential equations, for instance, in the papers [6] of De la Sen et al. and [7] of Li et al. Another way of handling an incubation period is by introducing another compartment. A comparison of these two approaches can be found in the work [8] of Kaddar et al. Other models allow for certain entities such as force of infection or incidence rate to be nonconstant. Such a model, in both a deterministic and a stochastic version, is considered in [9] by Lahrouz et al.

In this paper, we study a disease of the SIR type, prevailing in a population that can be regarded as consisting of two subpopulations. We compare it with similar models existing in the literature. We study stability of equilibrium solutions and optimal roll out of the vaccination. Such a study, in the case of a homogeneous population, was done in [10] by Zaman et al. For more complex population structures, there is a study by Piccolo and Billings [11]. A model similar to that of Piccolo and Billings has been studied in a stochastic setting in the work [12] of Yu et al. In [12], such a population is being referred to as a two-group population. A model of SEIR type for such a diversified population was proposed in [4] by Jia et al. In the latter paper, they analyse stability of solutions, but they do not consider vaccination. Our paper aims to follow the approach of [4], but for the SIR case and to include vaccination. Some papers have addressed epidemic models with pulse vaccination, for instance, an SIR model with pulse vaccination strategy to eradicate measles is presented in [13] by Agur et al. A model of an SIR epidemic in a two-group population, separated by age, is presented in the paper [14] of Acedo et al. They present a vaccination strategy similar to that in [13]. Much work in pulse vaccination has been done following on and inspired by [13]. However, two different diseases of SIR type may require completely different strategies for effective control of the disease. In this paper, we cater for those diseases for which pulse vaccination is not the best solution. We will assume the so-called proportional vaccination. A very interesting control problem is solved in the paper [15] of Tchuenche et al. In [15], the control vector is 3-dimensional, providing for a two-dimensional control on vaccination and a control on treatment. In the current paper, the control problem and its solution follow more closely along the lines of [10]. We obtain a simplification over [10] by observing that some of the pivotal costate variables vanish.

This paper is organized as follows. In Section 2, we formulate the model by way of a system of six ordinary differential equations. Then, we analyse the disease-free equilibrium and derive the threshold parameters in Section 3. In Section 4, we consider the optimal control problem, controlling vaccination on both the locals and the migrants. The percentages of susceptibles being vaccinated are taken as the control variables. We include a simulation. Finally, Section 5 has concluding remarks and offers a brief outlook on further research possibilities.

2. Model Formulation

To study the transmission of a disease in two interacting populations, we consider the total population with size 𝑁, as being divided into two subpopulations, the migrant subpopulation of size 𝑀, and the local subpopulation of size 𝐿. 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. The sizes of these classes change with time and will be denoted by 𝑆₀(𝑡), 𝐼₀(𝑡), 𝑅₀(𝑡), 𝑆₁(𝑡), 𝐼₁(𝑡), and 𝑅₁(𝑡). Let us agree henceforth to suppress the subscript (₀) for local population, writing simply 𝑆(𝑡) instead of 𝑆₀(𝑡), 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 𝑆𝐼𝑅 model as, for instance, in the paper [10] by Zaman et al. Let us normalize the variables, using the new variables 𝑠₁=𝑆₁/𝑀,𝑖₁=𝐼₁/𝑀,𝑟₁=𝑅₁/𝑀,𝑠=𝑆/𝐿,𝑖=𝐼/𝐿 and 𝑟=𝑅/𝐿. After normalizing our model, which we shall refer to as model (2.1) and (2.2), becomes as follows: 𝑑𝑠₁(𝑡)𝑑𝑡=𝑣₁−𝑣₁+𝑢₁(𝑡)𝑠₁(𝑡)−𝛽₁𝑖₁(𝑡)𝑠₁(𝑡),𝑑𝑖₁(𝑡)𝑑𝑡=𝛽₁𝑖₁(𝑡)𝑠₁(𝑡)−𝛾₁+𝑣₁𝑖₁(𝑡),𝑑𝑟₁(𝑡)𝑑𝑡=𝛾₁𝑖₁(𝑡)−𝑣₁𝑟₁(𝑡)+𝑢₁(𝑡)𝑠₁(𝑡),𝑑𝑠(𝑡)((𝑡𝑠(𝑡)(𝑡)𝑠(𝑡)𝑑𝑡=𝑣−𝑣+𝑢))−𝛽𝑖−𝛽₂𝑖₁(𝑡)𝑠(𝑡),𝑑𝑖(𝑡)(𝑡)𝑠(𝑡)𝑑𝑡=𝛽𝑖+𝛽₂𝑖₁(𝑡)𝑠(𝑡)−𝑖(𝑡),(𝑡)𝛾+𝑣𝑑𝑟(𝑡)(𝑡)(𝑡)𝑠(𝑡).𝑑𝑡=𝛾𝑖−𝑣𝑟+𝑢 (2.1)(2.2) Here 𝑣₁ and 𝑣 are the mortality rate (equal to the birth rate) in the migrant subpopulation, and the local subpopulation, respectively. The functions 𝑢₁(𝑡) and 𝑢(𝑡) are the percentages of susceptible individuals being vaccinated in the respective subpopulations per unit time. Individuals enter the recovered compartment at rates 𝛾₁ and 𝛾 for the respective subpopulations. Also 𝛽₁ 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 𝛽₂. The term 𝛽₂𝑖₁𝑠 models the influence of the migrant subpopulation onto the locals as in the paper [4] of Jia et al.

Figure 1 

Flow chart of two interacting populations.

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 𝑐₀ must be in step with the ratio 𝑀/𝐿. The feasible region for the system is the following set: Ω=𝑋∈ℝ⁶₊∶𝑋₁+𝑋₂+𝑋₃=1,𝑋₄+𝑋₅+𝑋₆=1. (2.3)

3. Equilibria and Their Stability

Equilibrium points are time-independent solutions to the given system of equations. Therefore, in this subsection, we assume 𝑢₁(𝑡) and 𝑢(𝑡) to be constant functions, 𝑢₁(𝑡)≡𝑢₁ and 𝑢(𝑡)≡𝑢. Stability properties of the equilibria are closely linked with the numbers 𝐾₁=𝛽₁𝑣₁𝑣₁+𝑢₁𝛾₁+𝑣₁,𝐾=𝛽𝑣(.𝑣+𝑢)𝛾+𝑣 (3.1) We shall prove that the basic reproduction ratio is the number 𝑅_𝑢,𝑢1=max{𝐾,𝐾₁}. This will follow from Proposition 3.3.

Notation 1. If 𝑢 and 𝑢₁ are both identically zero, then 𝑅_𝑢,𝑢1 will be written as 𝑅₀. For an equilibrium point 𝐸, the coordinates will be denoted by 𝐸_𝑠, 𝐸_𝑠1, and so on.

Remark 3.1. Suppose that in the model (1a, 1b) in [4] of Jia et al., we make the following modifications, transforming the model into SIR: replace the compartments 𝐸_𝑀 and 𝐼_𝑀 by a single compartment 𝐽_𝑀, and similarly replace 𝐸_𝐿 and 𝐼_𝐿 by a single 𝐽_𝑀,
Then the model takes the same form as our model (2.1) and (2.2), if in (2.1) and (2.2) we put 𝑢(𝑡)≡0, 𝑢₁(𝑡)≡0 and 𝑣₁=𝑣.
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 𝑢(𝑡)≡0 and 𝑢₁(𝑡)≡0, and let us further assume that 𝑣₁=𝑣.
If 𝑅₀<1, then the disease-free equilibrium 𝐹 with 𝐹_𝑠=1 and 𝐹_𝑠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 𝑣₁=𝑣, we can identify three possible equilibrium points.

Proposition 3.3. (a) If 𝑅_𝑢,𝑢1<1, then the disease-free equilibrium 𝐹 is locally asymptotically stable and its coordinates are 𝐹_𝑠1=𝑣₁𝑣₁+𝑢₁,𝐹_𝑖1=0,𝐹_𝑟1=𝑢₁𝑣₁+𝑢₁,𝐹_𝑠=𝑣𝑣+𝑢,𝐹_𝑖=0,𝐹_𝑟=𝑢.𝑣+𝑢 (3.2)
(b) If 𝐾₁<1 and 𝐾>1, then there is a unique feasible equilibrium 𝐵 with 𝐵_𝑠1=𝑣₁𝑣₁+𝑢₁,𝐵_𝑖1=0,𝐵_𝑟1=𝑢₁𝑣₁+𝑢₁,𝐵_𝑠=𝛾+𝑣𝛽,𝐵_𝑖=𝑣1𝛾+𝑣1−𝐾,𝐵_𝑟=1−𝐵_𝑠−𝐵_𝑟. (3.3)
(c) The endemic equilibrium 𝐷 has coordinates as follows: 𝐷_𝑠1=𝛾₁+𝑣₁𝛽₁,𝐷_𝑖1=𝑣₁𝛾₁+𝑣₁11−𝐾₁,𝐷_𝑟1=1−𝐷_𝑠1−𝐷_𝑟1, (3.4)𝐷_𝑠 is a root 𝑥 of the quadratic polynomial 𝑃(𝑥)=𝐶₂𝑥²+𝐶₁𝑥+𝐶₀ with 𝐶₀=𝛽₁𝑣²𝑣₁+𝛽₁𝑣²𝛾₁+𝛽₁𝛾𝛾₁𝑣+𝛽₁𝛾𝑣₁𝐶𝑣,₁=𝐶₀+𝛽₂𝑣₁𝛾₁𝛾+𝛽₂𝑣²₁𝛾−𝛽₂𝑣₁𝛽₁𝑣−𝛽₂𝑣₁𝛽₁𝛾−𝛽₁𝛾𝛾₁𝑢−𝛽₁𝛽𝛾₁𝑣+𝛽₂𝑢₁𝑣₁𝑣+𝛽₂𝑣₁𝛾₁𝑣+𝛽₂𝑣²₁𝑣−𝛽₁𝛽𝑣₁𝑣−𝛽₁𝑣𝑣₁𝑢−𝛽₁𝛾𝑣₁𝑢−𝛽₁𝑣𝛾₁𝑢+𝛽₂𝑢₁𝑣₁𝛾+𝛾₁𝛾+𝛾₁𝑣,𝐶₂=𝛽₁𝛽𝛾₁+𝑣₁(),𝐷𝑣+𝑢_𝑖=𝑣−(𝑢+𝑣)𝑠𝑣+𝛾,𝐷_𝑟=1−𝐷_𝑠−𝐷_𝑖. (3.5)

Proof. The given points 𝐹,𝐷,𝐵∈ℝ⁶ 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 𝐹 is 𝑎𝑊=¹𝛽00000₁𝑖₁𝑏₁0000−𝑢₁𝛾₁−𝑣₁0000𝛽₂𝑠0𝑎−𝛽₂𝑖₁−𝛽𝑠00𝛽₂𝑠0𝛽𝑖+𝛽₂𝑖₁𝑐0000𝑢𝛾−𝑣, (3.6) where 𝑎₁=−𝑣₁−𝑢₁−𝛽₁𝑖₁,𝑏₁=𝛽₁𝑠₁−𝛾₁−𝑣₁,𝑐=𝛽𝑠−𝛾−𝑣. (3.7) We set out to find the eigenvalues of 𝑊. This amounts to solving for 𝜆 in the equation, 𝑞₁⋅𝜆+𝑣₁⋅𝑞₂⋅(𝜆+𝑣)=0, (3.8) where 𝑞₁(𝜆) and 𝑞₂(𝜆) are the quadratic expressions below: 𝑞₁=𝜆+𝑣₁+𝑢₁+𝛽₁𝑖₁𝜆−𝛽₁𝑠₁+𝛾₁+𝑣₁+𝛽²₁𝑖₁𝑠₁,𝑞₂=𝜆−𝛽𝑠+𝛾+𝑣𝜆+𝑣+𝑢+𝛽𝑖+𝛽₂𝑖₁+𝛽²𝑠𝑖+𝛽𝛽₂𝑖𝑠. (3.9)(3.10) Now from (3.9) we can write 𝑞₁ in the form 𝑞₁=𝜆²+𝐴₁𝜆+𝐴₂, (3.11) where 𝐴₁ and 𝐴₂ are the constants: 𝐴₁=𝑣₁+𝑢₁+𝛽₁𝑖₁−𝛽₁𝑠₁+𝛾₁+𝑣₁,𝐴₂=𝑣₁+𝑢₁+𝛽₁𝑖₁𝛾₁+𝑣₁−𝛽₁𝑠₁+𝛽²₁𝑖₁𝑠₁. (3.12) Substituting the equilibrium values (at the point 𝐹) of 𝑠₁,𝑖₁,𝑠 and 𝑖, we can rewrite 𝐴₁=𝑣₁+𝑢₁−𝛽₁𝑣₁𝑣₁+𝑢₁+𝛾₁+𝑣₁,𝐴₂=𝑣₁+𝑢₁𝛾₁+𝑣₁−𝛽₁𝑣₁𝑣₁+𝑢₁. (3.13) The roots of 𝑞₁ have negative real parts if both 𝐴₁ and 𝐴₂ are positive. Now we note that 𝐴₂ is positive if and only if 𝛾₁+𝑣₁−𝛽₁𝑣₁𝑣₁+𝑢₁>0, (3.14) that is, 𝐾₁=𝛽₁𝑣₁𝑣₁+𝑢₁𝑣₁+𝛾₁<1. (3.15) If 𝐾₁<1, then also 𝐴₁>0. From (3.10) we have 𝑞₂ as follows: 𝑞₂=𝜆²+𝑣+𝑢+𝛽𝑖+𝛽₂𝑖₁−𝛽𝑠+𝛾+𝑣𝜆+𝛾+𝑣−𝛽𝑠𝑣+𝑢+𝛽𝑖+𝛽₂𝑖₁+𝛽²𝑠𝑖+𝛽𝛽₂𝑠𝑖. (3.16) Now let us define the coefficients 𝑄₁ and 𝑄₂ as 𝑄₁=𝑣+𝑢+𝛽𝑖+𝛽₂𝑖₁,𝑄−𝛽𝑠+𝛾+𝑣₂=𝛾+𝑣−𝛽𝑠𝑣+𝑢+𝛽𝑖+𝛽₂𝑖₁+𝛽²𝑠𝑖+𝛽𝛽₂𝑠𝑖. (3.17) By applying a similar analysis as for 𝑞₁, we find that the roots of (3.16) have negative real parts if and only if both 𝑄₁ and 𝑄₂ are positive, which is equivalent to the condition 𝐾=𝛽𝑣(𝑣+𝑢)𝛾+𝑣<1. (3.18) Therefore, the disease-free equilibrium is locally asymptotically stable if 𝐾₁<1 and 𝐾<1, that is, when 𝑅_𝑢,𝑢1<1.
((b) and (c)): The points are obtained by direct computation. Feasibility of 𝐵 is clear when 𝐾 and 𝐾₁ are as given in (b).

We include a computational example of an endemic equilibrium point 𝐷.

Example 3.4. Let us choose parameter values: 𝑣,𝛽,𝛾,𝑢,𝑣₁,𝛽₁,𝛾₁,𝑢₁,𝛽₂=(0.11,0.40,0.09,0.5,0.15,0.55,0.05,0.25,0.3). (3.19) Then we obtain 𝐾₁=1.03125 and 𝐷 has coordinates: 𝐷_𝑠1=0.36364,𝐷_𝑖1=0.02273,𝐷_𝑠=0.17726,𝐷_𝑖=0.00936. (3.20) We note that 𝑃(𝑥) also has a root 𝑥=0.50866, but this is not a feasible value for 𝐷_𝑠 since the corresponding 𝐷_𝑖 value (=−1.00141) is negative.

In line with the terminology of [4], we shall refer to the point 𝐵 as a boundary equilibrium. Stability analysis of the points 𝐵 and 𝐷 would take more effort than in the case of 𝐹, and could distract from the main purpose of this paper.

4. Optimal Vaccination Strategy

We wish to design optimal vaccination strategies 𝑢^∗(𝑡) and 𝑢^∗₁(𝑡), respectively, for the local population and the migrant population. We have six state variables 𝑠₁(𝑡),𝑠(𝑡),…,𝑟(𝑡). The variable 𝑢(𝑡) denotes the percentage of susceptible individuals being vaccinated per unit of time in the local population, and 𝑢(𝑡) is assumed to be bounded, 0≤𝑢(𝑡)≤𝛼≤1. A similar interpretation holds for 𝑢₁(𝑡), and we assume that for some constant 𝛼₁, 0≤𝑢₁(𝑡)≤𝛼₁≤1. Our optimal control problem amounts to minimizing the objective function below 𝐽𝑢(𝑡),𝑢₁(𝑡)=^𝑇₀𝑖(𝑡)+𝑐₀𝑖₁(𝑡)+𝑐𝑢²(𝑡)+𝑐₁𝑢²₁(𝑡)𝑑𝑡, (4.1) where 𝑐₀, 𝑐, and 𝑐₁ 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 𝐽(𝑢(𝑡),𝑢₁(𝑡)) subject to the system (2.1) and (2.2) of differential equations, together with the initial conditions 𝑠₁(0)=𝑠₁₀≥0,𝑖₁(0)=𝑖₁₀≥0,𝑟₁(0)=𝑟₁₀𝑠(0)≥0,=𝑠₀≥0,𝑖(0)=𝑖₀≥0,𝑟(0)=𝑟₀≥0, (4.2) and terminal conditions, 𝑠₁(𝑇),𝑖₁(𝑇),𝑟₁(𝑇),𝑠(𝑇),𝑖(𝑇), and 𝑟(𝑇) are free, while the control variables are assumed to be measurable functions that are bounded above 0≤𝑢(𝑡)≤𝛼≤1,0≤𝑢₁(𝑡)≤𝛼₁≤1. (4.3) The Hamiltonian for this problem is as follows: 𝐻𝑡,𝑠₁,𝑖₁,𝑟₁,𝑠,𝑖,𝑟,𝑢,𝜆₁,𝜆₂,𝜆₃,𝜆₄,𝜆₅,𝜆₆=𝑖(𝑡)+𝑐₀𝑖₁(𝑡)+𝑐𝑢(𝑡)²+𝑐₁𝑢₁(𝑡)²+𝜆₁(𝑡)𝑣₁−𝑣₁+𝑢₁(𝑡)𝑠₁(𝑡)−𝛽₁𝑖₁(𝑡)𝑠₁(𝑡)+𝜆₂(𝑡)𝛽₁𝑖₁(𝑡)𝑠₁(𝑡)−𝛾₁+𝑣₁𝑖₁(𝑡)+𝜆₃(𝑡)𝛾₁𝑖₁(𝑡)−𝑣₁𝑟₁(𝑡)+𝑢₁(𝑡)𝑠₁(𝑡)+𝜆₄(𝑡)𝑣−(𝑣+𝑢(𝑡))𝑠(𝑡)−𝛽𝑖(𝑡)𝑠(𝑡)−𝛽₂𝑖₁(𝑡)𝑠(𝑡)+𝜆₅(𝑡)𝛽𝑖(𝑡)𝑠(𝑡)+𝛽₂𝑖₁(𝑡)𝑠(𝑡)−𝛾+𝑣𝑖(𝑡)+𝜆₆(𝑡)𝛾𝑖(𝑡)−𝑣𝑟(𝑡)+𝑢(𝑡)𝑠(𝑡). (4.4) 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 𝜆₃(𝑡)=0=𝜆₆(𝑡)∀0≤𝑡≤𝑇, (4.5) and also satisfies the following system of differential equations: 𝜆̇₁=𝜆₁𝑣₁+𝑢₁+𝛽₁𝑖₁−𝜆₂𝛽₁𝑖₁,𝜆̇₂=−𝑐₀+𝜆₁𝛽₁𝑠₁−𝜆₂𝛽₁𝑠₁−𝛾₁−𝑣₁+𝜆₄𝛽₂𝑠−𝜆₅𝛽₂𝜆̇𝑠,₄=𝜆₄𝑣+𝑢+𝛽𝑖+𝛽₂𝑖₁−𝜆₅𝛽𝑖+𝛽₂𝑖₁,𝜆̇₅=−1+𝜆₄𝛽𝑠−𝜆₅𝛽𝑠−𝛾−𝑣, (4.6) with transversality conditions: 𝜆₁(𝑇)=0,𝜆₂(𝑇)=0,𝜆₄(𝑇)=0,𝜆₅(𝑇)=0. (4.7) Furthermore, the optimal vaccination strategy is given by 𝑢^∗(𝑡)𝜆=minmax^∗4(𝑡)𝑠^∗(𝑡),𝑢2𝑐,0,𝛼_∗1(𝑡)𝜆=minmax^∗1(𝑡)𝑠^∗₁(𝑡)2𝑐₁,0,𝛼₁. (4.8)

Proof. Existence of a solution follows since the Hamiltonian is convex with respect to 𝑢(𝑡) and 𝑢₁(𝑡). 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 𝜆̇_𝑖(𝑡) of the costate variables. Due to 𝑠₁(𝑇),𝑖₁(𝑇),𝑟₁(𝑇),𝑠(𝑇),𝑖(𝑇) and 𝑟(𝑇) being free, the following terminal conditions hold: 𝜆₁(𝑇)=0,𝜆₂(𝑇)=0,𝜆₃(𝑇)=0,𝜆₄(𝑇)=0,𝜆₅(𝑇)=0,𝜆₆(𝑇)=0. (4.9) We start off by observing that, 𝜆̇₃(𝑡)=−𝜕𝐻𝜕𝑟₁=−𝑣₁𝜆₃(𝑡),𝜆̇₆(𝑡)=−𝜕𝐻𝜕𝑟=−𝑣𝜆₆(𝑡). (4.10) This implies that 𝜆₃(𝑡) and 𝜆₆(𝑡) are of the form 𝜆₃(𝑡)=𝐴𝑒^−𝑣1𝑡,𝜆₆(𝑡)=𝐵𝑒^−𝑣𝑡, (4.11) for some constants 𝐴 and 𝐵, respectively. The terminal conditions 𝜆₃(𝑇)=0 and 𝜆₆(𝑇)=0, forces 𝐴 and 𝐵 to vanish. Therefore, 𝜆₃ and 𝜆₆ are identically zero, that is, 𝜆₃(𝑡)≡0 and 𝜆₆(𝑡)≡0 as claimed in the theorem.
Now we calculate 𝜆̇₁(𝑡)=−𝜕𝐻^∗𝜕𝑠₁,𝜆̇₂(𝑡)=−𝜕𝐻^∗𝜕𝑖₁,𝜆̇₄(𝑡)=−𝜕𝐻^∗,𝜆̇𝜕𝑠⁵(𝑡)=−𝜕𝐻^∗,𝜕𝑖 (4.12) and we obtain the equations as asserted in the theorem.
We now turn to the final part of the proof, which is about the form of the controls, 𝑢^∗(𝑡) and 𝑢^∗₁(𝑡). The function 𝑢^∗(𝑡) must optimize 𝐻. So we calculate 𝜕𝐻𝜕𝑢=2𝑐𝑢−𝜆₄𝑠. (4.13) Consider a fixed value of 𝑡. Now if 2𝑐𝑢(𝑡)−𝜆₄(𝑡)𝑠(𝑡) is zero for some value of 𝑢(𝑡) in [0,𝛼], then the given value of 𝑢(𝑡) is optimal. If for every number 𝑢∈[0,𝛼], we have 2𝑐𝑢−𝜆₄(𝑡)𝑠(𝑡)≥0resp.,2𝑐𝑢−𝜆₄(𝑡)𝑠(𝑡)≤0, (4.14) then we must choose 𝑢(𝑡)=0 (resp., 𝑢(𝑖)=𝛼). Thus, we must have 𝑢^∗(𝑡)𝜆=minmax^∗4(𝑡)𝑠^∗(𝑡).2𝑐,0,𝛼 (4.15) The function 𝑢^∗₁(𝑡) also must optimize 𝐻, and by a similar argument we obtain the stated expression for 𝑢^∗₁(𝑡).

Numerical Simulation
We present two simulations in the examples below, and we use the Runge-Kutta fourth-order method. For both of these examples, we use the same parameter values, but the initial conditions on the state variables will be different. The parameter values are as follows: 𝑐₀=1;𝑐=0.3;𝑐₁=0.2;𝑔=0.4;𝑑=0.0222;𝛽=0.09;𝛽₁=0.12;𝛽₂=0.02;𝛼=0.7;𝛼₁=0.8;𝜇₁=0.0222;𝛾₁=0.3;𝑇=300. (4.16) The time horizon of a control problem in epidemiology is usually dependent on economic factors such as budgeting, biological, and medical considerations, or even maybe influenced by political dynamics. For the purpose of our illustrative examples, the chosen value of 𝑇 is nominal.

Example 4.2. Consider the initial conditions 𝑠(0)=0.7;𝑖(0)=0.28;𝑟(0)=0.02;𝑠₁(0)=0.7;𝑖₁(0)=0.25;𝑟₁(0)=0.05. (4.17) 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).

Figure 2 

This plot shows the proportions of susceptible and infected individuals in the local subpopulation, according to Example 4.2.

Figure 3 

For the case of Example 4.2, we show the proportions of susceptible and infected individuals in the migrant group.

Figure 4 

The vaccination strategies for the two groups, Example 4.2.

Example 4.3. In this case we assume at time 𝑡=0 to have the local population to be infection-free. We consider the initial conditions 𝑠(0)=1;𝑖(0)=0;𝑟(0)=0;𝑠₁(0)=0.7;𝑖₁(0)=0.25;𝑟₁(0)=0.05. (4.18) 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).

Figure 5 

The vaccination strategies for the two groups, Example 4.3.

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.

Example 4.4. In Table 1 we use model parameters as for the simulations of Examples 4.2 and 4.3. Only the state variables and the controls are different. We take the values 𝑠(0)=1 and 𝑖(0)=0 as fixed (with 𝑟=1−𝑠−𝑖), and we select values for 𝑠₁(0) and 𝑖₁(0) as indicated in the table. For comparison with the case of optimal control, we consider the constant values 𝑢₁=0.7 and 𝑢=0.7 of 𝑢₁(𝑡), and 𝑢(𝑡), respectively, and in the last column we use the values 𝑢₁=0.7 and 𝑢=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.

Table 1 

5. 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.

Acknowledgments

The authors gratefully acknowledge the anonymous referees whose constructive suggestions led to a much improved paper. P. J. Witbooi acknowledges a research grant from the South African National Research Foundation (NRF) as well as funding from the South African Herbal Science and Medicines Institute (SAHSMI). K. Patidar acknowledges a research grant from the NRF.