Stochastic model of tuberculosis with vaccination of newborns

In this paper we have proposed a stochastic model for studying the dynamics of tuberculosis (TB) by incorporating vaccination of newly born babies. The total population in this model is subdivided in to four compartments, namely susceptible S(t), infected I(t), vaccinated newborns V(t), and recovered R(t). First, the developed model is expressed and analyzed by the deterministic approach. Since this approach neglects the randomness of the dynamics of the process, it has great limitations in the modeling process. To avoid this kind of issues, we transform the deterministic approach into a stochastic one, which is known to play a significant role by providing additional degree of realism compared to the deterministic approach. The analysis of the model is done employing both approaches. The invariant region, positivity of the solution, equilibrium points and their stability are checked. According to the analysis, we came to realize that the basic reproduction number for the stochastic approach is smaller than the deterministic one. We have conducted various numerical experiments and obtained interesting simulation results which indicate that a combination of increased newborn vaccination and appropriate treatment of infected individuals have a great contribution in combating TB. It is worth mentioning that the simulation results confirm the conclusion drawn from the qualitative analysis of the model.


Introduction
Tuberculosis (TB) is a bacterial disease caused by Mycobacterium tuberculosis and is usually acquired through airborne infection from someone who has active TB (smear-positive TB) while coughing, sneezing talking, and singing [5,6]. It usually attacks the lungs and can also affect other parts of the body like central nervous system, lymphatic system, brain and spine or kidney. The symptoms for active TB include fever, weight loss, chest pain, coughing up blood, feeling tired all the time, night sweets, loss of appetite, and anorexia [4].
Today among infectious diseases TB is a major public health concern in the developing countries. It affects all countries and age groups, but overall the best estimates for 2018 were that 89% of cases were adults (57% were male, 32% were adult women) and 11% were children. Furthermore, 8.6% were people living with HIV (72% of them in Africa) [22].
It is also one of the most serious public health challenges in Ethiopia, killing more than 30 thousand people every year [2]. Ethiopia ranks third in Africa and eighth among the 22 highest tuberculosis burdened countries in the world that collectively account for 80% of tuberculosis cases (Federal Ministry of Health [10]). It disproportionately affects young people: 58% of prevalent TB cases in Ethiopia are under 35 years of age, 39% of the estimated 32,000 deaths per year are concentrated among those from 15 to 64 years of age, leading to losses of family wage earners and parents of small children. This is an additional direct and indirect burden on Ethiopia's youth, who are the backbone of the current and future economy [3].
Individuals can prevent TB by proper treatment, that is, taking all doctor prescribed medicines for the specified time period, keeping all doctor's appointments, always covering the mouth with a tissue when coughing or sneezing, placing the used paper tissues in a plastic bag, then throwing them away, washing hands after coughing or sneezing, not visiting other people and not inviting them for a visit, as well as staying home from work, school or other public places, using a fan or open windows to move around fresh air, not using public transportation. Vaccination is another preventer [18]. Vaccine which is used for TB disease is Bacillus Calmette Guerin (BCG). In 1973, BCG vaccine was part of the immunization program [14]. BCG's efficacy on preventing pulmonary TB in adults is highly variable [21]. The BCG vaccination of newborns significantly reduces the risk of TB by over 50% on average [8]. Newborn usually refers to a baby from birth to about two months of age [9].
Mathematical modeling has been an important tool in analyzing the spread and control of infectious diseases and also in making decision as regards the intervention mechanisms for the control of infectious diseases [1]. Understanding the transmission characteristics of the infectious diseases in communities, regions, and countries can lead to better approaches to decrease the transmission of these diseases [17].
Modeling is necessary for the infectious TB for a number of reasons, for example, since TB has a complex and poorly understood natural history, the fact that it is not easy to conduct interventional research due to the lag between infection and the disease, the behavior of the susceptible population needs to be studied more, economic challenges in conducting interventions in low and middle income countries, and due to many unanswered questions about the impact of interventions [23]. Since transmission of the M. tuberculosis and infection of TB are influenced by various complex biological processes, the existence of randomness in the transmission dynamics of the disease could be believed [12]. Taking into consideration all these factors, today a stochastic model of TB is common by introducing perturbation in the parameters. Moreover, the authors of [16] constructed a deterministic SIR type model of TB epidemic and conducted a simulation by using fourth order Runge-Kutta method. The result shows that the TB spread can be restrained from an epidemic by sending down the spreading rate and increasing the recovery rate. A way to decrease the spreading rate is by keeping at a distance a TB-infected individual from the susceptible population, while to increase the recovery rate, maximal treatment needs to be conducted. But they did not analyze the model qualitatively. The work done by [15] developed two deterministic models of TB spread, which are SIR and SEIR. By implementing Lyapunov function method, it was obtained that the disease would disappear if the basic reproduction number was less than one, and it would persist if the basic reproduction number was greater than one. Additionally, the authors of [7] formulated a series of mathematical models to study the dynamics of TB. They formulated a distributed delay model to study the effect of long and variable periods of latency on the disease dynamics, which showed qualitative behavioral change at a critical value R 0 = 1. When R 0 < 1, a stable disease-free equilibrium existed and, for R 0 > 1, the disease-free equilibrium became unstable, and a unique endemic equilibrium existed. But they did not discuss sensitivity analysis of the basic reproduction number. Another model they formulated considered the role of nonadherence to drug taking by TB patients on the development and maintenance of antibiotic-resistant TB strain. Moreover, the authors of [13] studied a deterministic tuberculosis transmission dynamics model with vaccination and treatment for both high-risk latent and active TB infected classes. The reproduction number was calculated and the equilibrium points were described. They showed that the disease-free equilibrium point P * 0 is globally asymptotically stable when R 0 < 1, so that the disease dies out. Finally, they showed that an increase in the treatment and vaccination coverage gives rise to a decrease in the number of infected TB patients. Since the efficiency of the BCG vaccine is not complete, they assumed that some portion of vaccinated individuals will be susceptible to bacteria, but they did not consider the immunity waning in the vaccinated host. Moreover, the authors of [20] also developed a stochastic SVIR epidemic model of TB. Their investigations covered two important aspects: exponential stability of the diseasefree equilibrium and optimal control of vaccination. Regarding stability, the main result of their paper has a particularly simple formulation. Essentially, they said that they have almost sure exponential stability whenever the basic reproduction number of the underlying deterministic model is below unity. It will be good to know how an increase in vaccination rate would lead to better stability of the stochastic model. In their study numerical simulation enabled them to assess the feasibility of the option they followed for specific examples. We know that the transmission of different diseases is not the same, so that the developed model may not work for some diseases.
In this paper we developed an SVIRS stochastic model for the spread of tuberculosis considering vaccination of newborns. The paper is organized as follows. Section 2 introduces the formulation and description of the proposed tuberculosis model. In Sect. 3, the analysis of the model is discussed. Section 4 discusses a numerical simulation for the model. Finally, Sect. 5 contains discussion and conclusions.

Model formulation and description
The proposed model divides the entire population into four compartments or classes according to their disease status: susceptible S(t), infected I(t), vaccinated V (t), and recovered R(t).
The susceptible class, S(t), consists of individuals of all age groups in the population who have not come into effective contact with the Mycobacterium. The infected class, I(t), consists of individuals of all age groups infected with TB in the active stage; from the infected class an individual gets treatment and moves to the recovered class R(t). The vaccinated class, V (t), consists of individuals who had been vaccinated when they were newborns and still possess partial immunity against TB.
The vaccinated class is increased by birth with the recruitment π of probability p (0 ≤ p ≤ 1).
The susceptible class is increased by birth with recruitment rate π of probability (1p) and from vaccinated class V with rate b (b ≥ 0), as well as from the recovered class R with rate α.
In all subclasses, μ 1 is the natural death rate, μ 2 the disease-induced death rate for individuals in compartment I, β is the probability that susceptible (S) and vaccinated (V ) individuals become infected by one infectious individual per contact per unit of time.
And then, the rate at which the susceptible individuals are infected is βc; the rate at which the vaccinated individuals are infected is γβc. Here 0 ≤ γ ≤ 1. If γ = 0 then the vaccination protection efficacy is 100%, if γ = 1 then the vaccination protection effectiveness is 0, and 1-γ represents the reduction in infection risk due to vaccination effectiveness. Also r is the rate at which an infected individual leaves the infectious compartment I and joins the class R.
The model is based on the SVIRS transmission model. The total population size at time t is denoted by N(t), and therefore we have N(t) = S(t) + V (t) + I(t) + R(t).

Assumptions
1. Here we assume a homogeneous mixing of individuals in the population which means that every uninfected individual has an equal likelihood of being infected when coming into adequate contact with infectious individuals and that transmission of the infection occurs with a standard incidence rate. 2. We also assume that some recruits will emerge in the susceptible class, S, at a rate (1p)π and the vaccinated class V at a rate pπ . 3. The efficacy of the BCG vaccine is not complete so that some portion of vaccinated individuals will be infected by bacteria. 4. Due to the immunity waning in the vaccinated host, some portion of vaccinated individuals will be susceptible to bacteria with rate b. 5. Due to loss of immunity, we assume that recovered individuals move to the susceptible class at rate α. 6. We further assume that all parameters used in this model are positive. Considering the definitions, assumptions, and interrelations between the variables and the parameters, the basic dynamics of TB with vaccination for newborns is illustrated as a flow diagram in Fig. 1.
From the above flow diagram, the following system of differential equations is obtained: with initial conditions The deterministic approach has some limitations in the mathematical modeling of transmission of an infectious disease because biological processes involved in the dynamics of TB are stochastic rather than deterministic, so neglecting their built-in randomness may lead to misleading and erroneous results [12].
Stochastic differential equation (SDE) models play a significant role in various branches as they provide some additional degree of realism compared to their deterministic counterparts. Recently, many authors have introduced parameter perturbation into epidemic models and have studied their dynamics.
Therefore, in this study taking account of the effect of randomly fluctuating environment, we incorporate white noise in each equation of model (1). Suppose that some stochastic environmental factor acts simultaneously on each individual in the population. We use W i (t) for the mutually independent standard Brownian motions with W i (0) = 0 and σ i (i = 1, 2, 3, 4), the intensities of white noise. By introducing stochastic perturbation, the stochastic version corresponding to the deterministic model (1) takes the following form: (2)

Qualitative analysis of the model
In this section, the invariant region, positivity of solution, disease-free equilibrium point, endemic equilibrium point, basic reproduction number, as well as local and global stability of disease-free equilibrium point are discussed.

Positivity of solution
In this subsection, we show all solutions of the models (1) and (2) remain positive for future time if their respective initial values are positive.
Proof From the third equation of model system (1), Hence, this proves that I(t) > 0 for all t ≥ 0. From the first equation of model system (1), This proves that S(t) > 0 for all t ≥ 0. Similarly it can be shown that the remaining variables of the system are positive for all t ≥ 0.

Invariant region
In this subsection, we obtain a region in which the solutions of (1) and (2) are bounded.

Theorem 2 The feasible solution set {S, V , I, R} of the model remains bounded in the region
Proof For this model the total human population is N = S + V + I + R. Then, after differentiating N with respect to time and substituting the expression for dS dt , dV dt , dI dt , dR dt and from Eq. (1), we obtain If there is no infected individual in the population (i.e., I = 0), Eq. (3) becomes Solving equation (4) gives (1) and (2).

Disease-free equilibrium point
In this subsection we determine the equilibrium point at which the disease is eradicated from the population.
We consider the state where there is no infection, i.e., I = R = 0. By taking the second equation of model (1) and making it equal to zero, we have Since I = 0 in our consideration, Eq. (4) becomes Similarly, by taking the first equation of the model (1) and making it equal to zero, we have Since I = R = 0 in our consideration, Eq. (6) becomes (1p)π + bVμ 1 S = 0, and solving for S, we obtain Hence, Therefore, our model has DEF given by E 0 = (S 0 , V 0 , I 0 , R 0 ) = ( , pπ b+μ 1 , 0, 0).

Endemic equilibrium point
In this subsection we obtain the equilibrium point at which the disease persists in the population.
The endemic equilibrium (EE) of the system (1) is obtained by equating all equations of the model to zero: Adding all the equations of system (7) gives πμ 1 Nμ 2 I = 0, and solving for I, we have From the fourth equation in (6), we have, from (8). Hence, From the second equation in (6), we have Hence, From the first equation in (7) we have from (8), (9), and (10). Hence, Therefore, by (8), (9), (10), and (11) our model has EE given by

Basic reproduction number
The basic reproduction number (R 0 ) is defined as the average number of secondary infections caused by a typical infected individual during his entire period of infectiousness.
There are two types of basic reproduction, the deterministic and stochastic R 0 .

Basic reproduction number for deterministic model
By the principle of next generation matrix [19], we consider the following equation from the model (1): The next generation matrices are given by The next step is obtaining the Jacobian matrices of f and v by differentiating with respect to I and at the disease-free equilibrium E 0 .
The Jacobian matrices of f and v is obtained as F and V , respectively, where The Jacobian matrices F and V at the disease-free equilibrium pointE 0 are Then the product of F and V -1 becomes The eigenvalue of FV -1 can be obtained from Hence by the principle of next generation matrix, the dominant eigenvalue is the basic reproduction number.
Since in the limit dt → 0, (dt) 2 and dt dW 2 (t) →0 faster than (dW 2 (t)) 2 , which is dtfor(dW 2 (t)) 2 (due to the variance of Wiener process), The next generation matrices are F = [βCS + γβCV -1 2 σ 2 3 ]and V = [μ 1 + μ 2 + r]. Now F and V at the disease-free equilibrium become: Thus The product of F and V -1 is obtained as The eigenvalue of FV -1 can be obtained from By the principle of next generation matrix, the dominant eigenvalue is the basic reproduction number. Hence, Therefore, R S 0 < R D 0 because the stochastic version approximate reality closer than the deterministic one.

Theorem 3 The disease-free equilibrium point for the deterministic model is locally
asymptotically stable if R 0 < 1 and unstable if R 0 > 1.
Proof To prove the stability of the disease-free equilibrium pint for the deterministic model, we obtain the Jacobian matrix of the system at the disease-free equilibrium (E 0 ).
The Jacobian matrix of the system (1) is The Jacobian matrix of the system at the disease-free equilibrium is The eigenvalues of J (E 0 ) can be obtained from |J(E 0 ) -λI| = 0, that is, Then, From the first equation in (14), we have -(μ 1 + λ) = 0 ⇒ λ 1 = -μ 1 < 0. From the second equation in (14), we have-(b + μ 1 + λ) = 0 ⇒ λ 2 = -(b + μ 1 ) < 0. From the third equation in (14), we have According to Routh-Hurwitz stability criterion formula for the second-degree polynomial, one needs both roots in the open left half-plane (and then the system of the char-acteristic equation is stable) which is true if and only if both coefficients are greater than zero [20].
Since B > 0, we want to show C > 0. Now Therefore, the disease-free equilibrium (E 0 )is locally asymptotically stable if R D 0 < 1.
Proof By using Eq. (12), By integrating both sides, we have where for the martingale G(t) = t 0 σ 3 dW 3 (t),by the strong law of large numbers for martingales (X. Mao, 1997), we have lim t→∞ sup G(t) t = 0 almost surely. Dividing both sides of (15) by t and letting t → ∞, Hence, R S 0 < 1. Therefore, for the extinction of disease from a community R S 0 should be less than 1.

Global stability of disease-free equilibrium point
Therefore, the DFE is globally asymptotically stable in the feasible region if R D 0 < 1.

Sensitivity analysis of basic reproduction number
Sensitivity analysis for the basic reproduction number R 0 is being investigated to identify the parameters that have high impact on expansion or control of the disease in the community.
To perform this, we use normalized sensitivity index, ϒ For m = β, For m = c, For m = π , For m = γ , For m = p, For m = μ 2 , For m = r, For m = μ 1 , .
For m = d, From Table 1, we see that parameter β is with positive indices, which shows us that when we increase the value of β keeping other parameters constant then it increases the value of R 0 implying the disease expands in the community. Since the parameter c (contact rate) is positive and if we increase the value of c keeping other parameters constant, then the value of R 0 increases and as a result the disease expands in the population. Since the parameter π (the recruitment rate) is positive, increasing the value of π keeping other parameters constant, the value of R 0 is increased implying the disease expands in the community. Since the parameter γ (efficacy of vaccination) is positive, increasing the value of γ while keeping other parameters constant, the value of R 0 is increased and as a result the disease expands in the population. And the parameter b (the immunity waning rate in the vaccinated host) is positive, so when we increase the value of b keeping other parameters constant, then it increases the value of R 0 implying the disease expands in the community. We see that the parameter p (the fraction of the newborn vaccinated) is negative and when we increase the value of p keeping other parameters constant, it decreases the value of R 0 and then disease is eliminated from the community. While increasing the value of μ 2 (disease-induced death rate) keeping other parameters constant, the value of R 0 is decreased and as a result disease is eliminated from the community as it has a negative index. The parameter r (treatment rate for infectious individual) is negative and when we increase the value of r by keeping other parameters constant, the value of R 0 is increased and the disease is eliminated from the population. And μ 1 (natural death rate) has negative index which shows us that when we increase the value of μ 1 keeping other parameters constant, then it decreases the value of R 0 implying that the disease is eliminated from the community.

Numerical results and discussion
In this section, we present and analyze some simulation results for the dynamics of TB infection represented by models (1) and (2).
For this purpose, we used parameters of the model, namely p, c, r, and β, whose values are given in Table 2. We vary the values of these parameters and investigate their impact on the models. In addition, we consider the initial values of the variables of the model at time t = 0, i.e., S(1), V (0), I(0), and R(0).   due to the fact that stochastic models better represent physical situations compared to deterministic approaches. It can also be observed from these figures that the number of infected people keeps decreasing as time goes on, smoothly in the deterministic case and with some ups and downs in the case of the stochastic approach.

Effect of probability of contact rate on infected population
In this section we try to investigate the impact of the contact rate parameter β on the number of infected individuals I(t). Figure 3 presents the numerical results obtained by varying the value of β while keeping other parameters fixed. In the left side of Fig. 3 (i.e., deterministic model), when the value of β increased from β = 0.07 to β = 0.1, there is a significant and regular increase in the number of infected individuals. Moreover, when β = 0.2, the number of infected people quickly increases to 20 and start to go down a bit, but still manages to be higher than in the previous two cases. On the other hand, the results from the stochastic model, depicted in the right-hand side of Fig. 2, are also increasing, maintaining their zigzagging pattern (mimicking real-life case better) due to the randomness behavior. However, the overall outcome is that the number of infected people still increases significantly with increasing the value of β. Therefore, we can conclude that,

Effect of contact rate on TB infected population
In Fig. 4, the numerical results obtained by varying the value of contact rate c are plotted. We keep the other parameters fixed and vary c to facilitate investigation of the impact of c on I(t), the number of infected people. It is easy to see that, in both deterministic (left) and stochastic (right) cases, the curve is higher for larger values of c, implying the expansion of the disease in the community. That is, the number of infected people increases with increasing contact rate.

Effect of recovery rate on TB infected population
In this section, we investigate the impact of recovery rate r on the size of infected population I(t). The experimental results obtained by varying the value of r and keeping the other parameters constant are depicted in Fig. 5. The numerical results reveal that the number of infected individuals decreases with increasing the value of recovery rate r. It is also observed in the aforementioned figure that the decrement of infected individuals is smooth in the case of deterministic model (left-hand side of Fig. 5) whereas it is irregular in the stochastic case (right-hand side of Fig. 5). Hence, we can conclude that the increment in the value of the recovery rate r plays a great role in eliminating the disease from the community.  Figure 6 presents numerical results obtained by varying the value of parameter p (vaccination rate) and keeping the other parameters fixed. In the deterministic case, left-hand side of Fig. 6, the difference between the number of infected individuals becomes more obvious as time passes by, obviously the curve I(t) is slightly higher for smaller values of the vaccination rate p. Therefore, vaccinating the target population has a significant contribution in eliminating the disease. On the contrary, the stochastic case (right-hand side of Fig. 6) exhibits a very different behavior. For the first few moments, the number of infected people I(t) is seen to go upwards for all values of vaccination rate (but it is pronounced for p = 0.5). However, the value of I(t) turns to go down as time goes by. It is worth noticing that the increment in the value of the vaccination rate plays a vital role in eliminating the disease in the stochastic case as well.

Summary and conclusion
In this paper we have developed deterministic and stochastic models of tuberculosis with vaccination of newborns. In Sect. 2, we briefly described and formulated the tuberculosis infection model. In Sect. 3, we analyzed the models by obtaining the feasible region, positivity of the solution set, basic reproduction number, equilibrium points, and described their stability. The sensitivity of basic parameters, interpretation of the sensitivity index and global stability were also analyzed in Sect. 3. In Sect. 4, we presented and analyzed some simulation results. We started Sect. 4 by comparing the deterministic and stochastic TB infection trends in the community. According to the numerical results, we came to realize that the number of infected people keeps decreasing if one carefully combines vaccination with appropriate treatment. We have also observed that the stochastic model mimics and better represents real-life phenomena compared to the deterministic approach. Moreover, we have investigated the impact of the parameters β, c, r and p on both the deterministic and stochastic model. It was observed that increasing the probability of contact rate β contributes to the expansion of the disease whereas decreasing contact rate c, increasing vaccination rate r, and increasing the value of the recovery rate p all play vital roles in eliminating the disease from the community. Therefore, we recommend a combination of a decrease in contact between infected and susceptible individuals, increasing vaccination coverage, creating awareness, and proper treatment to effectively control TB infection.