Qualitative analysis of a stochastic SEITR epidemic model with multiple stages of infection and treatment.

We present a mathematical analysis of the transmission of certain diseases using a stochastic susceptible-exposed-infectious-treated-recovered (SEITR) model with multiple stages of infection and treatment and explore the effects of treatments and external fluctuations in the transmission, treatment and recovery rates. We assume external fluctuations are caused by variability in the number of contacts between infected and susceptible individuals. It is shown that the expected number of secondary infections produced (in the absence of noise) reduces as treatment is introduced into the population. By defining R T , n and R T , n as the basic deterministic and stochastic reproduction numbers, respectively, in stage n of infection and treatment, we show mathematically that as the intensity of the noise in the transmission, treatment and recovery rates increases, the number of secondary cases of infection increases. The global stability of the disease-free and endemic equilibrium for the deterministic and stochastic SEITR models is also presented. The work presented is demonstrated using parameter values relevant to the transmission dynamics of Influenza in the United States from October 1, 2018 through May 4, 2019 influenza seasons.

Introduction

Numerous mathematical models have been developed to study the transmission dynamics of emerging and re-emerging diseases (Diekmann, Heesterbeek, & Metz, 1990; Driessche & Watmough, 2002; Etbaigha, Willms, & Poljak, 2018; Feng, Towers, & Yang, 2011; Hollingsworth, Anderson, & Fraser, 2008; Huo, Chen, & Wang, 2016; Korobeinikov, 2009; LaSalle, 1976; Li, Xiao, Zhang, & Yang, 2012; Melesse & Gumel, 2010; Mendez, Campos, & Horsthemke, 2012; Tornatore, Buccellato, & Vetro, 2005; Otunuga, 2017; Otunuga, 2018; West, Bulsara, Lindenberg, Seshadri, & Shuler, 1979; Yang & Mao, 2013, Mummert & Otunuga, 2019).Without treatment of such diseases, infection advances in stages and infected individuals typically die within certain years. Several authors (Birrell, Presanis, & De Angelis, 2012; Hollingsworth et al., 2008; Korobeinikov, 2009; Melesse & Gumel, 2010; Otunuga, 2018) have studied extensively epidemic models with various stages of infection. Influenza has various stages of infection ranging from the contagious stage before any symptoms appear (period when the flu virus is entering and multiplying in only a few cells in the respiratory tract) to the stage when the flu virus has proliferated enough for the immune system to notice. The general incubation period for Influenza (typically known as the flu) varies for different individuals, usually between one to four days with average incubation period of about two days. This suggests that it is important to study the different stages of flu infection while studying transmission of infectious diseases.

Although it might be impossible to avoid certain infectious diseases, there are different strategies available that protect individuals from infection and treat disease once it has developed. It is of high importance to study how such disease reacts to treatments, and the analysis of treatment stages and treatment effects on infected individuals should be included in models describing the transmission dynamics of treatable diseases. Several programs such as the Biomedical Advanced Research and Development Authority have been developed by the U.S. Department of Health and Human Services to provide an integrated, systematic approach to the development and purchase of vaccines, drugs, therapies, and diagnostic tools necessary for public health medical emergencies.1

According to the work of Hu et al. (Hu, Nigmatulina, & Eckhoff, 2013), contact rates and patterns among individuals in a geographic area drive transmission of directly-transmitted pathogens, making it essential to understand and estimate contacts for simulation of disease dynamics. In their work, Grassly et al. (Grassly & Fraser, 2006) explains different causes of seasonality in infectious diseases of humans. They give different representations of the transmission rate based on the causes of seasonality in the infectious diseases. In this work, we study the global dynamics of a deterministic and stochastic SEITR epidemic model with multiple stages of infection and treatments. We assume the population is completely susceptible at the beginning of the epidemics and derive the measure of the power of an infectious disease to attack a completely susceptible population using the deterministic model. In the absence of noise, we compare mathematically the expected number of secondary cases of infection in the presence and absence of treatments and show that the number decreases as the treatment rate increases. We study the case where the transmission, treatments and recovery rates are assumed to be influenced by external fluctuations caused by variability in the number of contacts between infected and susceptible individuals due to weather patterns, school terms, etc. We assume fluctuations in the treatment rates may be caused by limited availability of drugs or effect of seasonality and this may result in fluctuations in the recovery rates. Such random variations can be modeled by a Gaussian white noise process causing the rate to fluctuate around a mean value. The external noise is able to modify the dynamical behavior of the model by transforming the deterministic SEITR epidemic model to a stochastic epidemic model. We derive the basic reproduction number in the presence of noise and analyze how the presence of noise in the transmission, treatments and recovery rates affects the number of infections produced by an infected individual. The paper is organized as follows. In Section 2, we formulate the deterministic model describing the transmission and spread of certain diseases, as well as its treatments and recovery. In Section 3, the existence of equilibrium points, and derivation of reproduction number using next generation method in the presence and absence of treatments are analyzed. Analysis of the effect of treatments and effect of dropping out of treatment on the number of infection produced by an infected individual are investigated analytically and numerically in Section 4. The local and global stability of the disease-free and endemic equilibriums are discussed in Section 5. By introducing noise in the transmission, treatment and recovery rates, we formulate and derive a stochastic model analogous to the deterministic model in Section 6. The effects of noise on the transmission, treatment and recovery rates, together with the existence and stability of the disease-free equilibrium point in the presence of noise are investigated analytically and numerically.

Deterministic model formulation

By assuming the human population is completely susceptible at the beginning of an epidemics and sub-dividing the total population, , into susceptible humans , exposed humans , infected untreated humans in stage j of infection, infected humans under treatment and in stage j of infection , and the recovered population , at time t, we investigate the transmission and treatment of certain infectious diseases. We assume the total human population satisfies and humans are recruited into the susceptible population at a rate . The general population is reduced by natural death at a rate μ. The population of susceptible humans is reduced by infection due to contact with infectious (untreated or treated) individual at a full rate . It is well known (Godoy et al., 2018) that influenza vaccination may not prevent infection but reduces the severity of the disease. The Center for Disease and Control 2 claimed that in randomized clinical trials, there was evidence that some influenza viruses developed resistance or reduced susceptibility to one or more influenza antiviral CDC recommended FDA-approved drugs like oseltamivir (Tamiflu), zanamivir (Relenza), peramivir (Rapivab), and baloxavir (Xofluza) . Several authors (Feng et al., 2011; Gani et al., 2005; Kretzschmar, Schim van der Loeff, Birrell, Angelis, & Coutinho, 2013; Liu and Zhang, 2011; Otunuga, 2018; Qiu & Feng, 2010) have considered introducing parameter that accounts for the reduction in infectiousness due to treatments among individuals in their model. In our model, we let be the reduced infectiousness due to treatment in stage j of infection and include the reduced rate due to treatment. Infected (but not yet infectious) individuals become untreated infectious individuals in stage 1 of infection at a rate π. Untreated infected individuals in stage k of infection migrate into stage of untreated infection at a rate and die of infection at a rate . These individuals receive treatment (and migrate to stage k of treated infected compartment) at a rate . Treated infected individuals in stage k of infection migrate to stage of treated infection at a rate and die of infection at a rate . Individuals that stop receiving treatment migrate to stage k of untreated infected compartment at a rate . Untreated and treated infected individuals in stage k of infection recover and migrate to the recovered compartment at a rate of and , respectively. The schematics describing the transmission described above is given in Fig. 1.

Fig. 1

Schematic diagram for the SEITR model. The circle compartments represent group of individuals.

The deterministic model governing S, E, , , R for , is described as follows:where the parameters in the model are described in Table 2, with . Since the limit , we consider the solution of the model (2.1) in the feasible regionwhere denotes set of nonnegative real numbers. For the rest of this work, we define . It can be shown that is positively invariant with respect to (2.1). We set the sizes of S, E, , , R, for as percentages by setting .

Table 2

Description of parameters for the epidemic model.

Parameter	Description
Λ	Recruitment rate into the population
β	Transmission rate of infection
h_k	Infectivity of untreated individuals in stage k of infection
ε_k	Reduced infectiousness due to treatment in stage k of infection
μ	Natural death rate
π	Infectious rate for exposed individuals
δ_k	Death rate associated with untreated infection in stage k of infection
δ¯_k	Death rate associated with treated infection in stage k of infection
τ_k	Treatment rate of infected individuals in stage k of infection
φ_k	Rate of dropping out of treatment in stage k
ρ_k	Transition rate from stage k to k+1 for untreated individuals
γ_k	Transition rate from stage k to k+1 for treated individuals
ψ_k	Recovery rate for untreated individuals in stage k of infection
η_k	Recovery rate for treated individuals in stage k of infection

Existence of equilibrium points in the presence and absence of treatments

We discuss the existence and stability of the equilibrium points of (2.1) in the presence and absence of treatment. Under certain conditions (which are discussed in (3.14) and Section 5), system (2.1) has two unique equilibrium points namely, the disease-free (denoted ) and endemic (denoted ) equilibrium points described as

The equilibrium points and are derived in 3.1, 3.2, respectively.

Disease-free equilibrium

The disease-free equilibrium of (2.1) has entries

In the following, we derive the measure of the power of an infectious disease to attack a completely susceptible population. It is the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual. This number, called the basic reproduction number and denoted by , is calculated explicitly considering n stages of infection and treatment. The endemic equilibrium, , is expressed in terms of . We also discuss a case where no treatment is received in the population and denote the corresponding reproduction number by . We show that in order for the number of infection to diminish to zero on the long run, appropriate parameters in the model must be controlled so that the number is at most one. That is, as long as the number of secondary infection produced by an infected individual is not more than one, the number of infections diminish to zero on the long run. Above the number , disease endemic presist.

Elimination threshold quantity, , in the presence of treatments

Define

In the presence of treatments, we write (2.1) in the formusing the next-generation matrix (Driessche & Watmough, 2002), where.

The derivatives and of and , respectively, are evaluated at and partitioned so that , , , , , and

The spectral radius of the matrix is given bywhere and satisfyand ; ; . We note here that for .where and are defined in (3.7) and the matrices and are coefficient matrices of the differential equationgoverning and in (2.1) for .

The reproduction number (3.6) can be re-written in matrix form as

Description of the derivation of

For a model with one stage of infection, if represent compartments and , respectively, then the entry of the inverse of the matrix defined in (3.5), and obtained asis the average time an individual introduced into compartment j spent in compartment i. It follows directly from (3.9) that the average time an individual introduced into the exposed compartment spent in the untreated infected compartment is , while the average time an individual introduced into the exposed compartment spent in the treated infected compartment is . An infected individual in the untreated and treated infected compartments and produces new infection in the exposed compartment E at a rate and , respectively. Thus, the number is the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual in compartment 1. In general, the average time an individual introduced into the exposed compartment spent in the untreated infected compartment is , while the average time an individual introduced into the exposed compartment spent in the treated infected compartment is . Hence, .

Reproduction number in the absence of treatment

Define

In the absence of treatment (that is, for ) we have , for and the reproduction number simplifies to the treatment free reproduction number given by

This is the reproduction number associated with the model without treatment

In a completely susceptible population receiving no treatment, we describe the quantity as the expected number of secondary infection produced by a typical untreated infected individual in a completely susceptible population.

The disease-free equilibrium point of (3.12) reduces to

Endemic equilibrium point, , in the presence of treatment

The endemic equilibrium of system (2.1) described in (3.1) is obtained asprovided , where and are defined in (3.7).

Endemic equilibrium in the absence of treatment.

In the absence of treatment, the endemic equilibrium reduces towhere is derived from (3.14) by setting and obtained asprovided .

Effect of treatment and dropping out treatment in the system

In this section, we study how receiving treatment and dropping out of treatment affect the system.

Effect of treatment of infection in the system

Consider the reproduction number corresponding to model (2.1) with j stage(s) of infection (derived by setting in (3.6)). Write as a function of for . We define the quantities and as the expected number of secondary infection produced by a typical infected individual (in a completely susceptible population with stages of infection) as treatment capacity goes to infinity and as no treatment is administered in stage i of infection, respectively.

We can show, after rigorous calculations, thatwhere , are defined in (3.7) for , , , and

Furthermore,

It follows from (4.4) that the derivative and the graph of concaves up for all if and only if . Likewise, and the graph of concaves down for all if and only if . By definition, we expect . This shows that in a population with j stages of infection, the number of secondary infection, , produced by an infected individual in a completely susceptible population decreases as the treatment rate increases.

Case where for all

Definewhere is defined in (3.10). For fixed , , we write (defined in (3.6)) as a function of τ. The number of secondary infection, , has the property:

The functionis a rational function of τ referred to as the relative elimination threshold. The graph of the function has y-intercept (following directly from Remark 3.1.1) and negative zeros. The vertical asymptotes are the negative vertical lines , for . Define

The function as . The value is the horizontal asymptote of . It measures the infection transmission potential when treatment capacity goes to infinity relative to the transmission potential when no treatment is administered. It follows from property of rational functions that if and if . This is represented in Fig. 2 below.

Fig. 2

Graphs of against τ for the cases where and .

Fig. 2 (a) and (b) show the trajectory of for the cases where and , respectively.

The quantity can be described as the expected number of secondary infection produced by a typical infected individual as the treatment capacity goes to infinity. From the description of in Remark 3.1.1, we expect , that is, we expect the expected number of secondary infection produced when the treatment capacity goes to infinity to be smaller than the expected number of secondary infection produced when no treatment is administered. This implies , so that . This shows that as the treatment rate increases, the expected number of infection decreases. The highest expected number of infection produced by an infected individual in a completely susceptible population is (which is attained when ) while the lowest expected number of infection is (attained as ).

Effect of dropping out of treatment

Write as a function of for . Using similar definition in Subsection 4.1, we define the quantities and as the expected number of secondary infection produced by a typical infected individual (in a completely susceptible population with stages of infection) as drop out treatment rate goes to infinity and as no one drops out of treatment in stage i of infection, respectively.

We obtain, after rigorous calculationswhere , are defined in (3.7) for , , , and

Furthermore,

It follows from (4.11) that the derivative and the graph of concaves down for all if and only if . Likewise, and the graph of concaves up for all if and only if . By definition, we expect . This shows that in a population with j stages of infection, the number of secondary infection, , produced by an infected individual in a completely susceptible population increases as the treatment dropout rate increases.

Assume for , and write . We see thatandwhere is defined in (3.7) for . The vertical asymptotes of the rational function are the negative vertical lines , for . Since is a rational function of φ whose numerator and denominator have the same degree, it follows that is an increasing function of φ if and only if , for . By definition, we expect . This shows that as the rate of dropping out of treatment increases, the expected number of secondary infection produced by an infected individual increases to .

Numerical results verifying the effects of treatment and dropping out of treatment on the number of infections

Here, we use relevant parameters to the transmission dynamics of influenza disease in the United States for the numerical simulations of the reproduction number as a function of the treatment and dropout rates. We set the life expectancy of the United States population to 80 years3 and the total population to be as of July 2018.4 Using the parameters collected from the Center for Disease Control and Prevention (CDC), the time from when a person is exposed and infected with flu to when symptoms begin is about 2 days, but can range from about 1 to 4 days5 and uncomplicated influenza signs and symptoms typically resolve after 3–7 days for the majority of people.6 Antiviral drugs, when used for treatment, can reduce symptoms and shorten sick time by 1 or 2 .

CDC7 estimates that, from October 1, 2018, through May 4, 2019, there have been million flu illness, million flu medical visits, thousand flu hospitalizations and thousand flu death. We define as a reduction factor in infectiousness (in stage j of infection) due to flu treatment and it reduces the infectious period to . For more information about the parameter , we refer readers to the work of Lipsitch et al. (Liu & Zhang, 2011), Feng et al. (Feng et al., 2011), Kretzschmar et al. (Kretzschmar et al., 2013) and . In their work, Lipsitch (Liu and Zhang, 2011) introduced a parameter which is the reduction in hazard of infection for an individual on prophylaxis. They claimed with probability , transmission is blocked and of those blocked infections, a proportion are only partially blocked. Using two infectious stages, we set , , , , , , , , , , , , , , , , , , , , . The value for the number of individuals under treatment is close to the number reported by Biggerstaff et al. (Biggerstaff, Jhung, Kamimoto, Balluz, & Finelli, 2012). According to the paper published by Tokars at al. (Tokars, Olsen, & Reed, 2018), between and of the U.S. population gets infected and develops flu symptoms each year. The value is approximately in this reported range. See Table 1, Table 2, Table 3 for parameter values and descriptions.

Table 1

Description of variables for the epidemic model.

Variable	Description
S	Population of susceptible individuals
E	Population of exposed individuals
I_k	Population of untreated infected individuals in stage k of infection
T_k	Population of treated infected individuals in stage k of infection
R	Population of individuals who recovered from disease

Table 3

Parameter values for the epidemic model: Case study Influenza.

Parameter	Description	Default Value	References
Λ	Recruitment rate into the population	180×365day^−1	CIA^3
β	Transmission rate of infection	∑j=1nβh_j=0.5	Feng et al. (2011)
h_k	Infectivity of untreated individuals in stage k of infection	0.5	(Feng et al., 2011; Roosa & Chowell, 2019)
ε_k	Reduced infectiousness due to treatment in stage k of infection	0.2	Feng et al. (2011)
π	Infectious rate for exposed individuals	1π=2 (days)	CDC^5
μ	Natural death rate	Λ	CIA^3
δ_k	Death rate associated with untreated infection	1.43×10^−4	Murphy, Xu, Kochanek, and Arias (2018)
δ¯_k	Death rate associated with treated infection		Assumed
τ_k	Treatment rate of individuals in stage k of infection	∑j=1nτ_j∈[0.05,0.2] (day^−1)	CDC^6
φ_k	Rate of dropping out of treatment in stage k	∑j=1n1φ_j=7 (days)	Assumed
ρ_k	Average duration of untreated infection	∑j=1n1ρ_j∈[3,7] (days)	CDC^6
γ_k	Average duration of treated infection	∑j=1n1γ_j∈[1,6] (days)	CDC^6
ψ_k	Recovery rate for untreated individuals in stage k of infection	∑j=1n1ψ_j∈[3,15] (days)	(Feng et al., 2011; Roosa & Chowell, 2019)& Assumed
η_k	Recovery rate for treated individuals in stage k of infection	∑j=1n1ηj∈[2,14] (days)	(Feng et al., 2011; Roosa & Chowell, 2019)& Assumed
S(0)	Initial susceptible Population		Assumed
E(0)	Initial Exposed Population		Assumed
∑j=1nI_j(0)	Initial Untreated Infected Population	37.4329.27	CIA^3, CDC^7
∑j=1nT_j(0)	Initial Treated Infected Population	20.1329.27	CIA^3, CDC^7
R(0)	Initial Recovered Population		Assumed

Fig. 3 (a) shows the graph of against . Fig. 3 (b) shows the graph of against . The graphs show that with no treatment, the reproduction number is , and as more treatment is introduced into the population the number of secondary infection reduces until it approaches , which is the least number of secondary infection that can be produced by an infected individuals when introduced into susceptible population. This is explained in Subsection 4.1.

Fig. 3

Effect of treatment on the reproduction number for and .

Fig. 4 (a) shows the graph of against . Fig. 4 (b) shows the graph of against . The graphs show that the number of secondary infection increases to as individuals drop out of treatment. This is explained in Subsection 4.2.

Fig. 4

Effect of dropping out of treatment on the reproduction number for cases and .

Fig. 5 (a) shows the graph of against and . Fig. 5 (b) shows the graph of against and .

Fig. 5

Effect of treatment and dropping out of treatment on the reproduction number for the cases and , and .

Existence and stability of equilibrium points

In this section, we discuss the endpoint behavior of the solution of (2.1). We give conditions under which the solution converges on the long run to the disease-free or endemic equilibrium.

Existence and stability of disease-free equilibrium in the presence of treatment

The following theorems show the condition for the local and global stability of the disease-free equilibrium, . We study condition(s) under which disease elimination exists on the long run. The idea presented here is similar to the work in Otunuga (Otunuga, 2018). To analyze the local asymptotic stability of , we linearize (2.1) about and show that the real part of all eigenvalues of the coefficient matrix of the linear associated system is negative.

Define . The linearization of (2.1) along the disease-free equilibrium is obtained aswhere with , , , , , , , , , , , , , , , and are defined in (3.5). We can express the characteristic polynomial of A in the formwhere is the square matrix formed by deleting the first row and column of A in (5.1) and r is the eigenvalue of A.

The real part of all eigenvalues of is negative if . One of the eigenvalues of is zero if and at least one of the eigenvalues is positive real if .

Proof. It suffices to show that the maximum real part of all eigenvalues of , denoted, , is less than zero if . To do this, we use relations and in (Plemmons, 1977) to show that the real part of each eigenvalues of the matrix is positive. The matrix can be written in the formwhere and are lower and upper diagonal matrices, respectively, with positive diagonals. The matrices and are computed rigorously as follows:where , and for , and can be simplified aswhere is the determinant operator, and { } are defined in (3.3) and (3.6), respectively. Since and for , it follows that and so for . Therefore, if , it follows from (5.4) that for , and the diagonal entries for . Since is a Z-matrix and the diagonal entries for , it follows from relations and in (Plemmons, 1977) that the real part of each eigenvalues of matrix is positive, which is in turn equivalent to . The determinant of the matrix is , which is the product of all -eigenvalues of . If , then , which means at least one of the eigenvalues of is zero. If , then , which means at least one of the eigenvalues of is positive.

The disease-free equilibrium of (2.1) is locally asymptotically stable if and unstable if .

Proof. The proof follows from (5.2) and Theorem 5.1.

The above theorem shows that if , the system approaches the equilibrium point whenever it starts somewhere near it in . The local stability of the disease-free equilibrium of system (3.12) without treatment follows immediately from Theorem 5.2 by setting for all . We state the theorem below without proof.

The disease-free equilibrium of (3.12) is locally asymptotically stable if and unstable if .

The following theorem gives the threshold under which disease elimination (considered independent of the initial conditions in ) exists.

The disease-free equilibrium of (2.1) is globally stable in the feasible region if .

Proof. Define the Lyapunov function bywhere is the set of positive real numbers, ϖ, and satisfyand reduces towhereand and are recurssive sequences defined by

The coefficients ϖ, and satisfy , for , and . It follows from (5.5) and (5.6) that the derivative of L computed along solution of (2.1) is

If , then . Thus, it follows from (5.6) and (5.7) that and are positive for andusing the fact that and . If , then if and only if , , and for all . Substituting this into the equation for in (2.1) shows that as . If , then if and only if . The largest invariant set of (2.1) contained in is the set . The global stability of follows from the LaSalle invariance principle (LaSalle, 1976).

The above theorem shows that disease can be eliminated on the long run from the population if parameters are controlled so that the elimination threshold is at most 1. This elimination is independent of the initial number of infection. The global stability of the disease-free equilibrium of system (3.12) without treatment follows immediately from Theorem 5.4 by setting for all . We state the theorem below without proof.

The disease-free equilibrium of (3.12) is globally asymptotically stable in the feasible region if .

Numerical results verifying global stability of disease-free equilibrium

Here, we use relevant parameters (given in Table 3) to the transmission dynamics of influenza disease in the United States for the numerical simulations of the number of susceptible, untreated infected, treated infected and recovered individuals satisfying the SEITR models (2.1) and (3.12).

Fig. 6 (a) shows the comparison of the trajectories of the number (in percentages) of exposed (), untreated infected () population in stage 1 of infection for model (3.12) (no treatment) with the trajectories of the number of exposed (E), untreated infected () and treated infected () population in stage 1 of infection for model (2.1) (with treatment) for the case where . Fig. 6 (b) shows the comparison of the trajectories of the number of exposed (), untreated infected () and () population in stages 1 and 2 of infection, respectively, for model (3.12) with the trajectories of the number of exposed (E), untreated infected (), () and treated infected (), () populations in stages 1 and 2 of infection, respectively, for model (2.1) with the case . It is clear from the graph that the introduction of treatment in the system reduces the number of exposed and infected individuals (that is, , and ) after some days. The number of exposed and infected individuals tends to zero on the long run and the number of susceptible individuals tends to 1. In this case, , , . and . The graph of the solution of system (3.12) converges to as . This confirms Corollary 5.5. Likewise, the graph of the solution of system (2.1) converges to as . This confirms Theorem 5.4.

Fig. 6

Graphs of comparison of deterministic trajectories of solution of system (2.1) and (3.12) for the cases where and , respectively.

Existence and stability of endemic equilibrium in the presence of treatment

The endemic equilibrium (given in (3.14)) of (2.1) exists if and only if and does not exist if . It becomes disease-free (that is, ) if .

Proof. It follows directly from (3.14) that , , , and for , if . The result for the case where follows from (3.14).

The following theorem gives the threshold for persistence of endemic (considered independent of the initial number of infection).

The endemic equilibrium of the system (2.1) is globally stable in the feasible region if

Proof. The existence of the endemic equilibrium follows from Theorem 5.6 if . Assume . Define the Lyapunov function bywhere , and , , are positive constants defined byand reduces towhere is given in (5.7). It follows from (5.9) and (3.14) that , , for , and .

The derivative of computed along solution of (2.1) is

Define

We havewherehence, from (5.10)–(5.11) and the fact that the arithmetic mean of a list of non-negative real numbers is greater than or equal to the geometric mean of the same list (Steele, 2004), it follows that ; ; ; ; ; for , and , for , and

Equality holds if and only if , for . Using (3.14) and the fact that satisfies (2.1), it follows that as . The largest invariant set of (2.1) contained in is the singleton . By the LaSalle’s Invariance Principle (LaSalle, 1976), it follows that is globally stable in the feasible region if .

The global stability of the endemic equilibrium of system (3.12) without treatment follows immediately from Theorem 5.7 by setting for all . We state the theorem below without proof.

The endemic equilibrium (given in (3.16)) of (3.12) is globally asymptotically stable if .

Numerical results verifying the global stability of and effect of treatment

Using two infectious stages, we use the same values of parameters given in Table 3 except that we set .

Fig. 7 (a) shows the comparison of the trajectories of the number of exposed (), untreated infected () individuals for model (3.12) with trajectories of the number of exposed (E), untreated infected () and treated infected () individuals for model (2.1) for the case where and . Fig. 7 (b) shows the comparison of the trajectories of the number of exposed (), untreated infected (), () individuals for model (3.12) with trajectories of the number of exposed (E), untreated infected (), (), and treated infected (), () individuals for model (2.1) for the case where and . It is clear from the graph that the introduction of treatment in the system reduces the number of exposed and infected individuals (that is, , and ) after some days. In this case, , , . and . The endemic equilibrium point for system (3.12) is () for the case and () for the case . Likewise, the endemic equilibrium points for system (2.1) for cases and are () and (), respectively. The graph of the solution of system (3.12) converges to as . This confirms Corollary 5.8. Likewise, the graph of the solution of system (2.1) converges to as . This confirms Theorem 5.7.

Fig. 7

Graphs of comparison of deterministic trajectories of solution of system (2.1) and (3.12) for the cases where and , with .

Fig. 8 (a) shows the graph of against and . Fig. 8 (b) shows the graph of against and . The graphs show that for fixed φ, as more (less) treatment is introduced into the population, the number of secondary infection reduces (increases) until it approaches (), which is the least (highest) number of secondary infection that can be produced by an infected individuals when introduced into susceptible population. This is explained in Subsection 4.1. Also, the number of secondary infection increases to as individuals drop out of treatment. This is explained in 4.1, 4.2.

Fig. 8

Effect of treatment and dropping out of treatment on the reproduction number for the cases and , with .

Derivation of stochastic model: effect of fluctuations and stability of disease-free equilibrium

In this section, we study the effect of noise on the transmission rates and infectivities, ; the treatment rates ; the recovery rates and in stage k of untreated and treated individuals, respectively, for . We assume the noise/external fluctiations in the system is caused by variability in the number of contacts between infected and susceptible individuals and such random variations can be modeled by a Gaussian white noise (Mendez et al., 2012). We also assume that fluctuations in the treatment rates may be caused by limited availability of drugs or effect of seasonality. This, in turn, causes fluctuations in the recovery rates. By allowing these rates to fluctuate about a mean value, we introduce external fluctuations in the model as follows:where and are independent Gaussian noise terms with zero mean, and , , and are the noise intensities, a measure of the amplitude of fluctuations, for . By substituting (6.1) into (2.1), we get a Langevin equation. The resulting equation is a stochastic differential equation. It is important to be able to interprete and evaluate the noise structure of this equation. The Itô approach on stochastic differential equation depends on Markovian and Martingale properties. These properties do not obey the traditional chain rule. Whereas, the Stratonovich approach obeys the traditional chain rule and allows white noise to be treated as a regular derivative of a Brownian or Wiener process. It has been suggested by several authors like West et al., Wong et al. (West et al., 1979; Wong & Zakai, 1965) that Stratonovich calculus is appropriate for Langevin equations with both internal and external noise. For this reason, by substituting (6.1) into (2.1), we extend the resulting equation to a Stratonovich stochastic model of the formwhere denotes the Stratonovich integral (Arnold, 1974); , , , , are standard Wiener process on a filtered probability space (); the initial process is measurable and independent of , , and , .

The Stratonovich dynamic model (6.2) is converted to its It’s equivalent (stated below) using the Stratonovich-It conversion theorem given in Bernardi et al. (Bernardi, Madday, Blowey, Coleman, & Craig, 2001) and Kloeden et al. (Kloeden & Platen, 1995).

The Itô stochastic differential equation having the same solution as the -dimensional Stratonovich stochastic differential equation (6.2) is given by

Proof. The proof follows using the Stratonovich-It conversion theorem given in Bernardi et al. (Bernardi et al., 2001) and Kloeden et al. (Kloeden & Platen, 1995).

Following similar approach presented in Otunuga (Otunuga, 2018), we can show, using the function, that and , where is a differential operator called the - operator defined bywhere and . It follows from Theorem 3.5 of Khasminskii (Rafail, 2012) that there exists a solution of (6.3) which is an almost surely continuous stochastic process and is unique up to equivalence if is independent of the processes , , , , . The solution described above can be shown to be nonnegative and in the feasible region using a similar idea presented in (Yang & Mao, 2013).

Equilibrium points and basic reproduction number in the presence of noise

The point defined in (3.1)–(3.2) is also the disease-free equilibrium of system (6.3). We calculate an equivalent of RT,n in (3.6), denoted by and derive threshold under which system (6.3) becomes disease-free on the long run. We first linearize the non-linear stochastic system about the disease-free equilibrium and study the stability of the solution of the linear system.

Define . The linearization of (6.3) about the disease-free equilibrium results inwhere with , , , , , , , , and defined in (5.1),, , , where.

, , and are defined in (3.3), , , and are matrices with entries , , , , , , , , , , and zero otherwise for . Define . The function satisfies the differential equation

The characteristic polynomial of can be expressed aswhere is the matrix obtained be deleting the first row and column of in (6.5), and is the eigenvalue.

Using the idea presented in Mendez et al. (Mendez et al., 2012) and in Section 3.1.1, we calculate the reproduction number with respect to the deterministic model (6.6) in the presence of treatment aswhere with , . We note here that the threshold is nonnegative provided

For the rest of this work, we assume condition (6.9) is satisfied.

We note here that the number reduces to if for all .

Condition (6.9) indicates that the noise intensities , and must not exceed the rates , and , respectively, for the model to be well defined.

Effect of noise in the treatment, and recovery rates

In this section, we study the effect of fluctuations in the treatment and recovery rates.

Effect of noise in the treatment rates

Assuming condition (6.9) is satisfied, and for , we wish to study how the number of infection changes due to changes in the treatment intensity rates . Define (given in (3.6)) and . It is easy to show that . As discussed in Subsection 4.1, the derivative if and only if , that is, is a decreasing function of if and only if . It follows that provided . The same result follows for the case where for all , that is, provided . An increase in the noise intensity in the treatment rate increases the number of secondary infection cases produced by a typical infective individual.

Effect of noise in the recovery rates of untreated infected individual

Assuming condition (6.9) is satisfied, and for . We wish to study how the number of infection changes due to changes in the untreated recovery intensity rates of infected individual. Write as a function of . Since the functions and are increasing function of t for , and can be expressed in terms of and , it follows from the increasing property of that . The higher the noise intensity in the untreated infected recovery rates, the higher the number of secondary infection cases produced by a typical infective individual.

Effect of noise in the recovery rates of treated infected individual

Assuming condition (6.9) is satisfied and for . By writing as a function of , we wish to show that . Since the functions and are increasing function of for , it follows that , that is, as the noise intensity in the recovery rate of treated infected individuals increases, the number of secondary infection cases produced by a typical infective individual increases.

Numerical analysis

We use the parameters presented in Table 3 to verify the results claimed in Sub1, 2, 3, 4, 5, 6.

Fig. 9 (a), (b) and (c) show the graphs of , and against (fixing ), (fixing ) and (fixing ), respectively. Fig. 9 (d) shows the graph of against and . The trajectories of these graphs suggest that the higher the intensity of noise in the treatment rate, recovery rates of untreated and treated infected individuals, the higher the number of secondary infections produced by an infected individuals when introduced into a susceptible population.

Fig. 9

Effect of noise on treatment rates and recovery rates of untreated and treated infected individuals for the case .

Fig. 10 (a) and (b) show the graphs of against and and against and . The trajectories of these graphs suggests that the higher the intensity of noise in the treatment rate, recovery rates of untreated and treated infected individuals, the higher the number of secondary infections produced by an infected individuals when introduced into a susceptible population.

Fig. 10

Effect of noise on treatment rates and recovery rates of untreated and treated infected individuals for the case .

Stability of infection-free equilibrium of (6.3)

In this section, we discuss conditions for stability of the infection-free equilibrium of (6.3) in the presence of noise. We study the conditions for stochastic stability of the disease-free equilibrium of the linear associated system (6.5) and later use Theorem A.2 in (Tornatore et al., 2005) to extend the result to that of the nonlinear system (6.3).

Assume condition (6.9) is satisfied. The real part of all eigenvalues of is negative if .

Proof. The proof follows from (6.9) and Theorem 5.1 by setting , , , , and into matrix A in (5.1).

Writing the system of non-linear stochastic differential equation (6.3) in terms of reduces tofor , where and are defined in (3.3).

Let F and be the drift and diffusion coefficients of the linear system (6.5), respectively, and f and g the drift and diffusion coefficients of the non-linear system (6.10), respectively. We give a theorem concerning the global stability of the disease-free equilibrium point by showing that Theorems A.1 and A.2 of Tornatore et al., (2005) is satisfied with respect to systems (6.5) and (6.10).

The disease-free equilibrium of the system (6.3) is globally asymptotically stable in the feasible region if .

To prove this, we first show that if , the trivial solution of the linear stochastic differential equation (6.5) is assymptotically stable and later show that the drift and diffusion coefficients and , respectively, of the nonlinear system (6.10) satisfy the inequalityin a sufficiently small neighbourhood of , with a sufficiently small constant ξ.

Proof. If , it follows from Theorem 6.2 that the real part of all eigenvalues of is negative. Hence, there exist a diagonal matrix (with positive diagonal entries, say, ) and a real number such that for every nonzero vector (see relation of (Plemmons, 1977)) . Let be a vector satisfying the linear system (6.5) and define by

Let such that , , for . Using (6.4), the -operator defined in (6.4) satisfies

Let and be and , respectively. Then . It follows from Theorem A.1 of Tornatore et al., (2005) that the trivial solution of (6.5) is asymptotically stable. We deduce from this result that if the initial condition (in ) of system (6.5) is near 0, then the solution approaches on the long run if . To prove the global stability of the solution of (6.10) (equivalent to the disease-free equilibrium of (6.3)), we choose sufficiently small in a neighbourhood of so that and reduces towhere The global stability result follows from Theorem A.2 of (Tornatore et al., Vetro).

Numerical verification of global stability of infection-free equilibrium points for the stochastic model

Fig. 11 (a) shows the trajectories of E, and satisfying model (6.3) for the case where and . Fig. 11 (b) shows the trajectory of E, , , , satisfying model (6.3) for the case where and . In this case, and .

Fig. 11

Graphs of stochastic trajectories of solution of system (6.3) for the cases where and , respectively, and .