On the reproduction number and the optimal control of infectious diseases in a heterogenous population

The effect of infectious diseases cannot be overemphasised. The continuing surfacing of the infectious diseases gives the stakeholders a great concern. In this paper, the nature of the spread of Ebola virus outbreak in West Africa in 2014 is studied. We develop a model that analyses the spread of infectious diseases, and the reproduction number is determined by using the next generation matrix method. Finally, the effects of treatment of the infected individuals and vaccination of the susceptible population as the control strategies are looked into. The optimal control system showed that the combination of the two strategies proved more effective.


Introduction
Infectious diseases are deadly diseases of humans and animals which account for about one-tenth of the causes of deaths worldwide. These deadly diseases cause major change to the size of the population of any community during their outbreaks. These are diseases like HIV, tuberculosis, malaria, including Ebola virus disease which West African countries have just witnessed. The dynamics of the population change due to the high mortality rate caused it to be a global problem to the practitioners of public health.
However, some diseases incubate in their hosts for a period of time before the hosts become infectious [1], and in a heterogenous population, this is a cause for concern.
Several modelers have done justice to the various models of infectious diseases to show their transmission and spread. The authors in [2] modified the S-E-I-R model and developed an S-E-I-H-F-R model to study the spread and transmission of Ebola virus to include those at hospital and funeral but assumed homogeneous population, which according to the authors was too simple, which may not have effects in countries where the structure of the community favours infection in households, and recovered individuals were not returned to the population. Still on the model formulation of infectious disease, a model with nonlinear incidence to analyse the transmission of Zika virus was formulated by [3]. The model considered three forms of nonlinear forces of infection. In 2013, [4] studied the global stability conditions of two models having relapse and bilinear standard incidence rate. He discovered that the disease produced deaths in the infectious class. He used the construction of a Lyapunov function for the systems. The analysis of epidemiological data of the 2014 Ebola virus disease outbreak in Nigeria was done by dividing the population into five categories: susceptible; exposed; infectious and symptomatic individuals; hospitalised individuals; and those individuals separated from isolation centre after recovery or those that died as a result of the disease. The model developed was S-E-I-H-P from S-E-I-R [5]. In the work of [6] in (2018), they studied the stability analysis of infectious disease using an S-E-I-R model incorporating recruitment rate. It was discovered that the model was stable both locally and globally. Many other researchers have equally worked on the stability and control analysis of infectious disease. The authors in [7] examined the significance of short-term human travellers on transmission dynamics of Malaria. The chance of an individual developing breast cancer depends on the level of the immune system, the efficacy of the anti-cancer drug and the rate at which the ketogenic diet is taken [8,9]. [10] in 2015 used an SIR epidemic model to model the 2014 Ebola outbreak in West Africa and introduced vaccination to the susceptible as part of the control system. They concluded that vaccination is a very efficient factor in reducing the number of infected individuals. On the account of controlling the spread of these deadly diseases, [11] examined the possible impact of using combination of three control strategies to combat the spread of Malaria. The effective treatment and prevention of Schistomiasis infection assists in eradication of Malaria [12]. Vaccination combined with effective treatment reduces the spread of malaria disease [13]. The control strategies of using treated bed-nets, medication and insecticides spray have great impact on the control of malaria disease [14]. To this end, the main aim of this work is to model the spread of Ebola virus disease in West Africa with the incubation period by applying mathematical models of the systems of ordinary differential equations that describe the dynamics of the outbreak of 2014 Ebola virus in some parts of West Africa, to determine the reproduction number and to better understand the dynamics and find optimal control strategies to reduce the effect of the disease spread. We considered two control strategies: treatment of the infected individuals and vaccination of the susceptible.

Model formulation
The dynamics of the spread is described by a system of ordinary differential equations. Because Ebola virus incubates in the body of the host for some time before manifesting, the population is divided into four (4) compartmental groups: the susceptible population (S), the exposed population (E), the infected population (I) and the recovered population (R). We introduce a recruitment rate into the population which is dynamic, i.e. migration and birth are possible. To describe the rate of change of the susceptible, the rate at which the susceptible is being populated by the recruitment rate is reduced as the infected individual comes into contact with them with a constant rate of "a", which means the rate of change of the susceptible population with respect to time is now (2.1) The equation of the exposed population will be those leaving the susceptible population to populate it and will later reduce by those that will die natural death during the outbreak, which we represented as constant of μ, and will also be reduced by those that will not be affected upon exposure and those that will be infected and progress into infected population with their constant rates being represented as b and δ respectively. Thus the equation of the exposed population expressing the rate of change of the exposed with respect to time is now The equations of the infected and recovered groups will relate in the same way as above, i.e. and That is, the rate of change that describes the equation of the infected population is those leaving the exposed group which will be reduced by those that will recover from the infection, those that will die as a result of the disease and those that will die naturally with the infection. The equation of the recovered will be those that recovered from the infection, added to those that recovered upon exposure to the infection and reduced by those that will recover and return to the susceptible population at a constant rate of c. So our model is now represented by system (2.5) of Eqs. (2.1)-(2.4) as follows: where a is the rate of infection, b is the rate of recovery, μ is the natural death rate, c is the rate at which the recovered individuals move back to the susceptible to populate it and δ is the rate at which the exposed individuals progress to the infected class.

Positivity and boundedness of solutions
From system (2.5), we have the state variables to be dS dt | S=0 = + cR, dE dt | E=0 = aSI, dI dt | I=0 = δE and dR dt | R=0 = bI + bE. This implies that all the above rates are nonnegative on bounding planes R 4 + . Furthermore, we next show that the solutions of system (2.5) are positively invariant and attracting. We note that Then the derivative of (3.1) taking into account (2.5) iṡ Since the term eI in (3.2) is nonnegative, Hence, the following positively invariant bound set Thus, it is sufficient to consider the dynamics of system (2.5) in with initial conditions

0
It is known or described as the anticipated cases of newly created or fresh cases of infections from one individual that is infectious in a population that is wholly susceptible through the entire length of the infection period denoted by R 0 . It defines the dynamical behaviour of the model, whether the disease dies out or persists in the system. If R 0 < 1, the infection in one individual cannot reinstate itself so the pathogen dies out (stable diseasefree population). If R 0 > 1, the number of infectious persons increases and the disease persists, and if R 0 = 1, there is an equilibrium: the endemic and the disease-free being equal. By the use of next generation matrix, where F is the matrix of the newly created infection, V is the matrix of transferred infection / the movement of individuals from one compartment to another.
V -1 is the inverse of matrix V. So, Let V cf be the cofactor The dominant eigenvalue of G is the basic reproduction number denoted by R 0 , i.e. |G -λI| = 0, where I is the identity matrix. Therefore, , The dominant eigenvalue is .
Therefore, the basic reproduction number of our model is . (3.5)

Local stability of disease-free equilibrium
The stability analysis of system (2.5) is discussed in this section. The Jacobian matrix of system (2.5) at P o is Two eigenvalues of (3.6) are λ 1 = -μ, λ 2 = -(c + μ), while the other two eigenvalues are obtained from the 2 by 2 matrix We next apply the Routh-Hurwitz condition, namely (i) Tr(A) < 0 (ii) Det(A) > 0 to show that the remaining two eigenvalues are real and negative. Consequently, Tr(A) = -(μ + b + δ) -(b + e + μ) = -(2μ + 2b + δ + e) < 0, The above result is summarized as follows.

Theorem 1
The DFE Po of system (2.5) is locally asymptotically stable if R 0 < 1, otherwise unstable.

Local stability of the endemic equilibrium
We observed that the endemic equilibrium of system (2.5) can be re-written in terms of the basic reproduction number R 0 as follows: , Therefore, I can be written as so also, .
Similarly, R can also be written as The following result follows immediately.

Theorem 2
The endemic equilibrium P * of system (2.5) is unique and positive if R 0 > 1.
Next, we discuss the stability of the endemic equilibrium point P * by obtaining the Jacobian matrix of system (2.5) at P * , thus where , , , and r 0 = (c + μ).
It then follows that all the eigenvalues of (3.9) are real and negative if the following condition holds: (i) Tr(J(P * )) < 0 and (ii) Det(J(P * )) > 0. Now, Tr(J(P * )) = -am 2 ( Thus all eigenvalues of (3.9) are real and negative. The foregoing can be stated as follows.

Theorem 3
The endemic equilibrium P * of system (2.5) is locally asymptotically stable if R 0 > 1, otherwise unstable.

Optimal control of the model
This is to find an optimal control strategy U throughout the length of 0 ≤ t ≤ t f such that the number of infected individuals I and that of exposed individuals E are minimised while minimising the cost of control U. Let U 1 be the cost of the treatment of infected and exposed individuals and U 2 be the cost of the vaccine for them. From system (2.5), the model now becomes The objective function now is

Analysis of the optimal control problem of infectious diseases
Let H be the Hamiltonian function, and using the maximum principle given by Pontryagin to derive necessary conditions for the optimal control problem, we have where λ i , i = 1 . . . 4, are the adjoint or the co-state variables.

The adjoint conditions
In this section, the adjoint condition is obtained as follows: with the boundary conditions at the final time

The optimality conditions
The optimality condition in this section is hereby found as follows: At the absolute minimum, the slope of the function is zero. Then at ∂H ∂U 2 = 0, So, for the control system, the conditions have been set for U 1 and U 2 .

Numerical simulation and discussion of results
In this section, we present the numerical simulation of the model using the parameters of [15]. These are presented in Figs. 1-5.

Results
This paper has studied the dynamics of the spread of infectious diseases using the outbreak of 2014 Ebola virus disease as a case study and presented a model which allowed move- The optimal control strategy 1 using treatment only Figure 4 The optimal control strategy 2 using vaccination only ment in and out of the population during the outbreak of the epidemic. The reproduction number R 0 was found to be .
It was noted that when R 0 < 1, the disease died out, and when R 0 > 1, the disease persisted in the system which is endemic. The results in this study were in perfect agreement with the parameters estimated in [16]. In Fig. 1, the variation of the susceptible population was presented. The susceptible population dropped drastically at the onset of the epidemic and reduced further as the rate of infection increased.
In the study we also tested the effect of treatment in Fig. 2. It showed that when the treatment of infected humans was intensified from 0.0 to 0.8, the number of infected humans dropped drastically. Figures 3  and 4 showed the optimal control strategies of treatment and vaccination of susceptible population. We found out that using a combination of treatment and vaccination in Fig. 5

Figure 5
The optimal control system using the combination of the two strategies proved to be more effective than using only one control, as we found out in Figs. 3 and 4, in that the disease can die out in about 100 days of the outbreak.

Conclusion
The reproduction number which determines disease eradication or persistency was found using the next generation matrix. The optimal control of the model was studied using Pontryagin's maximum principle. Control strategies for the model were considered using time-dependent control: treatment U 1 (t) and vaccination U 2 (t). The simulation result revealed that the combination of the two controls proved to be more effective. In the control of infectious diseases such as Ebola, it is advised that the medical practitioners should start preparing for vaccines which can curb the effect of Ebola virus and act as an immunity to the susceptible population. Also during the outbreak, the treatment of the infected individuals should be combined with the vaccination to have a more robust strategy to curb the effect.