A nonlinear multi-population behavioral model to assess the roles of education campaigns, random supply of aids, and delayed ART treatment in HIV/AIDS epidemics

1.   Introduction

HIV (human immunodeficiency virus) is the agent that causes AIDS (acquired immunodeficiency syndrome). HIV weakens the human immune system by attacking the body's special defensive cells, CD4 cells (T cells), against infection. With the depletion of the T-cells in the body, the human being becomes vulnerable to secondary infections or infection related cancers ^[1]. There is no cure for HIV, and infected persons live with the virus for life. According to the WHO ^[2], nearly $37.9$ million people live with HIV by the end of 2018.

While over two thirds of the infected population lives in sub-Saharan Africa, the most affected sub-populations globally are men who have sex with men (MSM), drug users, prisoners and people living in closed settings, sex workers and their clients, and transgender people.

There is biomedical treatment available against HIV. The main treatment used to control HIV infection is called antiretroviral therapy or ART ^[1]. With proper ART, the viral load is reduced and may become undetectable. People treated properly with ART, and with undetectable viral loads, live healthy long lives, and exhibit effectively no risks of transmitting the virus to HIV-negative persons ^[1]. Moreover, if HIV is diagnosed early and properly treated, the individuals live nearly natural lives as uninfected individuals. However, without proper treatment, or in the absence of treatment, the infected persons progress to AIDS, and this can occur in $2$ to $15$ years ^[2]. Thus, there is a critical time delay $\tau_{2}$ ^[3] to diagnosis and the onset of proper treatment necessary for a healthy longer lifespan; there is also the natural time delay of $2$ to $15$ years, $\tau_{1}$, until the onset of AIDS.

As remarked by WHO ^[2], between $2000$ and $2018$, $13.6$ million lives were saved due to ART. Moreover, new HIV infections decreased by $37\%$, and HIV related deaths also decreased by $45\%$. And these achievements were the result of national HIV programs supported by civil society and international organizations and partners. National and multinational campaigns against HIV have a history of success in many communities. For example, the success of Uganda in reducing the prevalence of HIV since the late 1980s is attributed to governmental information, education and treatment campaigns against HIV ^[4,5,6]. Several other studies have investigated the role of information and education campaigns (IEC) on the prevalence of HIV/AIDS ^[7,8].

Other forms of national, multinational and multilateral assistance in combating HIV/AIDS have been in the form of aids. Top multilateral organizations fighting against HIV/AIDS such as Global Fund ^[9], and PEPFAR ^[10] provide official development assistance (ODA)^{* [11]}, funding and supporting large and small scale projects globally designed to prevent and treat HIV/AIDS. Such assistance in treatment or prevention has saved many lives against the disease ^[9,10,11,12]. On the downside, the supply of ODAs are sometimes sporadic and random, and this can upset disease control. For instance, recently, the American government discussed options to cut back funds and spending on ODAs against HIV/AIDS ^[13]. This announcement led to some studies ^[14] to forecast long-term effects of such policy change on global HIV/AIDS prevalence. These studies have projected insignificant monetary savings by the US government, and rather devastating clinical and epidemiological impacts on the global spread of HIV/AIDS. Thus, it may be helpful to mathematically explore the effects of either cutting back funds on ODAs, or randomly supplying ODAs such as national, multinational and multilateral assistantships, to fight against the prevalence of HIV/AIDS.

^* All forms of domestic or foreign monetary support and other aids provided to combat HIV are referred to as ODAs in this paper.

HIV/AIDS prevalence is also strongly associated with poverty. As clearly remarked by ILOAIDS ^[15], HIV/AIDS is at the same time the cause and an outcome of poverty, and poverty is both a cause and an outcome of HIV/AIDS. Indeed, HIV/AIDS impoverishes poor households that expend savings and assets to afford medical care; HIV/AIDS retards economic growth, and saps out the vitality of healthy economies through a weakened labor force, as the skilled, and experienced are killed by the disease.

On the other hand, poverty drives and exposes the poor and unemployed into unhealthy employment situations, such as, migratory manual laborers who are in pursuit of temporary or seasonal migratory work, and because of challenges in daily living and basic accommodations, they may secure accommodations with sex workers, who are HIV positive ^[16]. Poverty also drives and exposes the poor and unemployed into risky sexual behaviors and practices such as prostitution, especially among women who may exchange sex for food ^[15,16], or women who are economically dependent on their partners and are more susceptible to risky sexual favors in order to please their partners ^[15,17]. Thus, a suitable HIV/AIDS model should account for the effects of poverty and living standards either explicitly or implicitly.

Mathematical models have certainly advanced understanding about the dynamics of HIV/AIDS epidemics ^{[4,7,8,18,19,20]}. Mathematical epidemic models with information intervention also advance understanding about the role of information and education in changing attitudes and behavior that lead to disease control ^{[4,7,8,21,22]}. Joshi et. al. ^[4] studied a SIRE model for HIV/AIDS epidemic in Uganda. Two susceptible classes with change of behavior namely- practicing abstinence or condom use, emerge from a general uneducated susceptible class via interaction with information in the community, from information and education campaigns (IEC). The three susceptible classes nevertheless experience the disease from interactions with infected persons, but with different per capital disease transmission rates. In their model, the density of education in the community has a separate dynamics. Huo et.al. ^[18] also studied a HIV/AIDS epidemic model with a treatment class with ART. The treated class does not transmit the disease, and it can either relapse to the active infectious class, if treatment is discontinued, or progress into full blown AIDS.

Employing similar reasoning in the studies ^[4,18], a more generalized HIV/AIDS epidemic model is studied in this paper. It is assumed that the IECs in the community educate the adult population with multiple preventive measures (more than the two in ^[4]) against HIV. Moreover, the response to the education results in multiple distinct behaviorial changes that can be visibly characterized into distinct sub-susceptible classes $S_j, j = 0, 1, \dots, n, \quad n\in \mathbb{N}$, exhibiting lowered disease vulnerabilities. Furthermore, the density of information or education, $Z(t)$, at any time $t$ has a separate dynamics. In this model, sporadic random supply of ODA and varying poverty levels in the population are also studied. Adding randomness in the supply of ODAs and poverty rates leads to a stochastic differential equation model for HIV/AIDS with treatment, multiple behaviorial changes, and time delays to onset of treatment and full blown AIDS.

In this paper, the impacts of IECs, the random supply of ODAs and delayed ART treatment on HIV/AIDS control are investigated via conducting a stochastic stability analysis of the system of differential equations for the HIV/AIDS epidemic. Moreover, a theoretical exploration of these HIV/AIDS epidemic factors is conducted via applying the epidemic model to Uganda HIV/AIDS data, wherein a multiple linear regression model for the incidence rate of the disease is derived, and utilized to simulate the disease dynamics.

This paper is organized as follows. In Sections 2–4, the stochastic epidemic model is derived. In Section 5 model validation results are presented. In section 6, stochastic stability of the model is conducted, sensitivity analysis of the BRN $\mathfrak{R}_{0}$ and discussion of the disease control conditions are given. Numerical simulation results for Uganda HIV/AIDS epidemic are presented in Section 7.

2.   Description of model

The HIV/AIDS epidemic model is based on the following assumptions summarized into definitions.

Model-Assumption 2.1. Notations and abbreviations

(A). The notation $I(k, n), k\geq 0$ represents the set of consecutive natural numbers between $k$ and $n$. E.g. $I(0, n) = \{0, 1, 2, \ldots, n\}$.

(B). The following abbreviations are used: information and education campaigns: IECs; official development assistance: ODA; antiretroviral therapy: ART; disease free equilibrium: DFE; stochastic solution process: SSP; basic reproduction number: BRN $\mathfrak{R}_{0}$; Pre-exposure prophylaxis: PrEP.

Model-Assumption 2.2. Population structures and human behavioral categories:

(A) Sexually active adults in a community are considered. The primary means of HIV transmission is sexual contact. Vertical transmission is not considered, and alternative means of transmission such as contact with infected needles are not considered. The IECs are designed to change adult sexual behaviors, especially of the susceptible population.

It is assumed that the IECs lead to $n\in \mathbb{N}_{ > 1}$ distinct behavioral categories based on safe-sex measures that are taught and learnt via the IECs ^{[6,23,24,25,26]}. Examples of the safe-sex behavioral categories include preventive measures such as: abstinence, mutually monogamous relationships (i.e. be faithful to partners), condom use, use of lubricants, voluntary medical male circumcision, counseling, harm reduction interventions for people who use drugs and all other distinct preventive measures that reduce vulnerability and transmission rates of HIV. For a comprehensive WHO recommended HIV prevention package and other HIV preventive measures, see ^{[6,23,24,25,26]}. It is assumed that nearly all adults learn and actively practice at most one distinct measure at a time. That is, everyone practices a $j^{th}$ measure, $j\in I(0, n)$, where the category $j = 0$ represents the state of "naivety", where no safe-sex measure is practiced. All persons who practice at least two measures at a time are collectively grouped into one of the $(n+1)$ behaviorial categories $j\in I(0, n)$.

The total sexually active human population $N(t)$ is decomposed into five major states namely: the susceptible state $S(t)$, which is vulnerable to HIV infection; the HIV infected individuals not receiving ART $I(t)$; the treatment state $T(t)$, representing all HIV infected individuals receiving ART treatment; the AIDS state $A(t)$, representing all HIV infected persons in the advanced stages of their HIV infection, and experiencing full symptoms of AIDS; the removed state $R(t)$, representing all those who practice safer sexual behavior, and fully protected from HIV infection. The removed state can consist of individuals who are adhering to HIV prevention measures such as PrEP (see Model-Assumption 2.3). The susceptible state $S(t)$ is further decomposed into $(n+1)$ distinct susceptible states $S_{j}(t), j\in I(0, n)$, based on the $n+1$ IECs behavioral categories above. Hence,

where

Indeed, the state $S_{0}(t)$ represents all susceptible individuals at time $t$, who do not practice any HIV preventive measures, either due to negligence or limited knowledge about HIV preventive measures. The states $S_{j}(t)$, $j\in I(1, n)$ represent all susceptible individuals who through the IECs, learn and actively practice the $j^{th}$ major HIV/AIDS preventive measure, where $j\in I(1, n)$. Note that for obvious reasons, it is not necessary to decompose the other states $I, T, A, R$ into the $(n+1)$ behaviorial categories.

It is also assumed that there is a constant influx $B$ per unit time of susceptible adults in the population. Moreover, all new individuals into the population are susceptible of type $S_{0}$.

Model-Assumption 2.3. Information density and interaction rates:

(B) The density of information at anytime $t$ is denoted $Z(t)$. The "naive" susceptible individuals $S_{0}(t)$ modify their behavior to become at most one of $S_{j}(t)$, $j\in I(1, n)$ after receiving education, $Z(t)$, about the disease at time $t$, at the effective response rate $\gamma_{j}\equiv \gamma_{S_{0}S_{j}}, j\in I(1, n)$.

The rate per unit time at which the susceptible individuals in class $S_{0}(t)$ change their behavior into class $S_{j}(t)$ is given by the expression $\gamma_{j}S_{0}(t)H_{j}(Z(t))$, where $H_{j}, j \in(1, n)$ is a nonlinear function describing the response of the susceptible class $S_{0}(t)$ to the density of information $Z(t)$ in the population.

It is also assumed that some individuals in the susceptible class $S_{0}(t)$ experience the highest impacts of the IECs after interacting with information $Z(t)$ at effective contact rate $\gamma_{0}\equiv\gamma_{S_{0}R}$. The impacts of the education obtained from the IECs result to reform their sexual behavior, and produce actions all through their lives that never result in contracting the disease. In other words, these individuals are considered to be immune to the disease and removed, $R(t)$, at time $t$. Indeed, according to WHO and CDC ^[24,25,26], proper training in the use of HIV preventive medications such as Pre-exposure prophylaxis (PrEP), leads to reduction of HIV infection risk by 99%. Thus, individuals in the population who are properly trained and adhere to correct daily use of PrEP can be considered into the removal class $R(t)$.

Thus, $\gamma_{0}S_{0}(t)H_{0}(Z(t))$ is the rate per unit time at which $S_{0}(t)$ individuals are removed ($R$), by fully reforming their sexual behaviors and attitudes via interacting with information $Z(t)$, and practicing safe-sex measures that will never lead to HIV infection.

Indeed, it is assumed that the effects of the content of the information related to HIV prevention measures from the IECs, initially rises in susceptible individuals who have not heard it, due to excitement about knowledge of new preventive measures, then saturates due to familiarity with the content of the information. The properties of $H_{j}$ are given in Assumption 2.1. Using ideas in ^[27,28,29] we adopt assumptions for the nonlinear function $H_{j}, j\in I(0, n)$.

Assumption 2.1. $A1:$ $H_{j}(0) = 0$; $A2$: $H_{j}(Z)$ is strictly monotonic on $[0, \infty)$; $A3:$ $H_{j}\in \mathcal{C}^{2}(\mathbb{R}_{+}, \mathbb{R}_{+})$ and $H^{''}_{j}(Z) < 0$; $A4:$ $\lim_{Z\rightarrow \infty}H_{j}(Z) = C_{1}$, $0\leq C_{1} < \infty$; and $A5:$ $H_{j}(Z)\leq Z, \forall Z > 0$.

In addition, using ideas from ^[22], the rate $\gamma_{j}$ can be expressed as $\gamma_{j} = \nu_{j}a_{j}$, where $a_{j}$ is the interaction rate by which individuals of type $S_{j}$ change their behavior, and $\nu_{j}\in [0, 1]$ is the response intensity.

The density of information in the population $Z (t)$ at time $t$ from the IECs is assumed to grow at a rate that is proportional to the number of infected individuals $I(t), T(t), A(t)$ in the population. Furthermore, the growth rate exhibits nonlinear character in response to the number of infected individuals of all states $I(t), T(t), A(t)$ in the system. This growth rate per unit time is represented by the function $F_{Z}(I(t), T(t), A(t))$.

Apparently, the content of information in the IECs relates to the infected classes $I(t), T(t), A(t)$, as preventive measures are taught against these states. Moreover, the rate of supply of information, $F_{Z}(I(t), T(t), A(t))$, saturates over time with increase in $I(t), T(t), A(t)$ (see B). Also, it is assumed that the effectiveness or strength of the content of the information in the IECs degrades at the rate $\mu_{Z}$. A special form for the function $F_{Z}(I(t), T(t), A(t))$ is given in (3.10).

Note that there are several different ways to quantity HIV information density $Z(t)$ in the population over time. For instance, in ^[4], the amount of information present at anytime is a function of the number of national, international, governmental and non-governmental organizations involved in IECs against HIV in Uganda at any time.

Model-Assumption 2.4. Disease transmission and generalized standard incidence rates of the disease:

(C) The active HIV infectious class $I(t)$ passes infection to all susceptible states $S_{j}, j\in I(0, n)$. However, because of preventive measures learnt from the IECs, the sub-classes $S_{j}, j\in I(1, n)$ experience a reduced rate of transmission from $I$, than the class $S_{0}$. That is, at the rate $\beta_{j}\equiv \beta_{S_{j}}, j \in I(0, n)$, the interaction between the susceptible state $S_{j}$ and infectious state $I(t)$ results in HIV transmission. The rate $\beta_{j}\equiv \beta_{S_{j}}, j \in I(0, n)$ represents the average number of effective contacts (i.e. sufficient contacts to transmit disease) per person per unit time. Moreover, $\beta_{0}\geq \beta_{j}, j\in I(1, n)$. In Section 7, multiple linear regression is applied to model the disease transmission rates $\beta_{j}\equiv \beta_{S_{j}}, j \in I(0, n)$ over time, for a given HIV/AIDS data for Uganda.

Recognizing the viewpoints regarding the incidence rates of human epidemics ^[30], a nonlinear generalization of the standard incidence rate $\frac{\beta_{j} S_{j}(t)I(t)}{N(t)} = \beta_{j} S_{j}(t)i(t), i(t) = \frac{I(t)}{N(t)}$ is considered, where $N(t)$ is the total population at time $t$.

Indeed, $i(t) = \frac{I(t)}{N(t)}$ is the fraction of the HIV infectious persons in the population at time $t$. Observe that when $N(t)$ is a constant, then $\frac{\beta_{j} S_{j}(t)I(t)}{N(t)} = \beta_{j} S_{j}(t)i(t) = \theta (t)\beta_{j} S_{j}(t)I(t)$, where the fraction $0 < \theta(t) = \frac{1}{N(t)} < 1$, reflects that the incidence rate rises linearly with respect to the infectious state $I$, and this pattern is unsuitable for most human epidemics ^[31]. Also, when the total population size $N(t)$ grows or declines proportionately with a rise or a drop in disease transmission in the population, respectively, i.e. $N(t)$ and the infectious state $I(t)$ both change (increase or decrease) proportionately, then the fraction $i(t) = \frac{I(t)}{N(t)}$ is constant overtime, and the standard incidence rate $\frac{\beta_{j} S_{j}(t)I(t)}{N(t)} = \beta_{j} S_{j}(t)i(t)$ no longer reflects the true incidence rate of the disease in the population.

Thus, to increase flexibility in the standard incidence rate to represent more real life scenarios, a nonlinear incidence function $G_{j}$ is introduced with assumptions in Assumption 2.2. The properties of the nonlinear function $G_{j}$ in Assumption 2.2 signify a psychological response from the susceptible classes $S_{j}, j\in I(0, n)$, where more susceptibles apply more appropriate preventive measures and actions that limit contacts with infectious persons, as the HIV infectious state $I(t)$ increases in the community over time.

The modified nonlinear incidence rate of HIV in the state $S_{j}$, is given by the expression $\beta_{j} S_{j}(t) G_{j}(i(t)), j\in I(0, n)$. The function $G_{j}$ satisfies the assumptions in Assumption 2.2.

Using ideas in ^[27,28,29] we adopt assumptions for the nonlinear function $G_{j}, j\in I(0, n)$.

Assumption 2.2. $A1:$ $G_{j}([0, 1])\subseteq \mathbb{R}^{+}$, $G_{j}(0) = 0$; $A2$: $G_{j}(I)$ is strictly monotonic on $[0, 1]$; $A3$: $G_{j}\in \mathcal{C}^{2}([0, 1], \mathbb{R}_{+})$ and $G_{j}^{''}(i)\leq 0$; $A4:$ $\lim_{i\rightarrow \infty}G_{j}(i) = C_{2}, 0\leq C_{2} < \infty$; and $A5$: $G_{j}(i)\leq i, \forall i > 0$.

An example of a function in the $G$-family in Assumption 2.2 is $G(x) = \frac{ax}{1+bx}, x\geq 0$. Indeed, to illustrate further, Figure 1 depicts the behaviors of the standard incidence function $i(t) = I(t)/N(t)$, and a modified standard incidence function $G(i(t)) = \frac{ai(t)}{1+bi(t)}$, as the number of infectives increase, where $I(t)\in [0, 1000000]$ and $a = 0.05, b = 10$.

Note that disease transmission from the AIDS state $A(t)$ and from treatment state $T(t)$ are not considered ^[4,18]. Indeed, it is assumed that individuals in the AIDS stage of their HIV infection exhibit the typical visible symptoms of the disease ^[32], and as a result they are either aware of their disease status and take precautionary measures via abstinence to not infect others, or they are unwell due to symptoms of the disease.

Model-Assumption 2.5. Random supply of ODAs and delays in the disease dynamics:

(D) The effects of randomness in the supply of ODAs will be assessed through the parameters $\varepsilon_{j} \in(0, 1)$, and $\bar{\varepsilon}_{j}\equiv 1-\varepsilon_{j} \in(0, 1)$, where $\varepsilon_{j} \in(0, 1), j\in I(0, n)$ represents the proportion per unit time of newly infected individuals from the class $S_{j}, j\in I(0, n)$ who do not receive ART treatment, and consequently progress to the AIDS state $A(t)$ after the time delay $\tau_{1}$, and $\bar{\varepsilon}_{j}\equiv 1-\varepsilon_{j} \in(0, 1)$, $j\in I(0, n)$ is the other proportion per unit time of the newly infected from the class $S_{j}, j\in I(0, n)$, who after a time delay $\tau_{2}$ following exposure to HIV, proceed to be tested, begins ART treatment and become the state $T(t)$, respectively.

In the absence of noise in the proportions $\varepsilon_{j} \in(0, 1)$, and $\bar{\varepsilon}_{j}\equiv 1-\varepsilon_{j} \in(0, 1)$, it is expected that a constant significant supply of ODAs in a community enables more people to be easily tested and to afford ART. And as a result, the proportion $\bar{\varepsilon}_{j}\equiv 1-\varepsilon_{j} \in(0, 1)$ increases, and the delay until ART begins, $\tau_{2}$, decreases on average. Also, if there is little or no supply of ODAs in the community, then more people cannot afford testing and ART, and consequently the proportion $\varepsilon_{j} \in(0, 1)$ increases.

The question of how random and sporadic supplies of ODAs in the community affect the HIV/AIDS epidemic dynamics will be answered by introducing white noise into the parameters $\varepsilon_{j} \in(0, 1)$, and $\bar{\varepsilon}_{j}\equiv 1-\varepsilon_{j} \in(0, 1)$ (see (4.2)).

Clearly, the time delay until the onset of ART, $\tau_{2}$, varies with available resources from ODAs, and also varies with the attitudes of newly exposed individuals in the population. Indeed, some communities attach huge social stigmas to HIV/AIDS, and as a result many people are dissuaded by such stigmas from testing and beginning early ART. Sometimes late commencing of ART may be the result of simple ignorance about the benefits of early testing and ART. These combinations of attitudes towards HIV can delay testing and diagnosis of HIV, and consequently, lead to delayed ART. Also, as remarked earlier, progression from HIV without treatment to full-blown AIDS can occur after time delay $\tau_{1}$ varying between 2 to 15 years. Therefore, distributed time delays $\tau_{1}$ and $\tau_{2}$ are considered to represent the variabilities in the delays, with probability density functions $f_{\tau_1}, t_0 \leq \tau_1 \leq h_1$ and $f_{\tau_2}, t_0 \leq \tau_2 \leq h_2.$

Model-Assumption 2.6. Withdrawal from treatment and developing full-blown AIDS:

(E) The only form of treatment considered is ART. In the advanced stages of HIV without treatment, the infectious individual $I(t)$ develops full-blown AIDS $A(t)$ after the natural incubation period $\tau_{1}$. Individuals who begin treatment $T(t)$ in the advanced stages of HIV, where considerable damage to T-cells has occurred, can still progress to full-blown AIDS $A(t)$ at the per capita rate $\alpha_{TA}$, as treatment fails. Individuals diagnosed early with HIV, and receiving treatment can discontinue ART due to a random sporadic limited supply of ODAs in the community, or due to personal self-limiting factors against the ART. Such individuals in the state $T(t)$ relapse to the active infectious state $I(t)$ at the rate $\alpha_{TI}$.

Model-Assumption 2.7. The per capita death rates:

(F) All individuals of all states in the population die naturally at the per capita rate $\mu_{k}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$, while the infected classes-$I, T, A$ die from the disease at the rate $d_{k}, k\in\{I, T, A\}$. Since in most natural settings the natural death causes are uniform across all disease states, then $\mu_{k} = \mu, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$, where $\mu$ is a constant. However, the distinct notation for the natural deathrates $\mu_{k}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$ will be used.

Model-Assumption 2.8. poverty indicators in the disease dynamics:

(G) The effects of poverty in the HIV/AIDS epidemic dynamics are implicitly assessed through the per capita natural deathrates $\mu_{k}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$. Indeed, it is apparent that a major indicator of the living standards of a community is the life expectancy, and assuming that natura deaths occur independently with homogenous Poisson rates $\mu_{k}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$, then the average life expectancy is related to natural death rate by $\frac{1}{\mu_{k}}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$. Furthermore, higher living standards correlate with richer economies, and small values for $\mu_{k}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$, and vice versa. Moreover, even within a given geographical region, the life expectancy varies. In this paper, the random poverty rates within the community over time are assessed by introducing independent white noises into the natural death rates $\mu_{k}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$.

3.   Derivation of the epidemic model in the absence of noise

It follows from the assumptions (A)–(G) in Model-Assumptions 1-8 above that a compartmental framework depicting the transitions between the different states of the population is given in Figure 2, and the deterministic HIV/AIDS epidemic dynamic model with treatment and information intervention follows immediately.

and

where the initial conditions are given in the following: define $h = \max{(h_{1}, h_{2})}$,

where $\mathcal{C}((-\infty, t_{0}], \mathbb{R}_{+})$ is a Banach space of continuous functions endowed with the uniform norm

The reader is directed to 9, for the detailed derivation and interpretation of the distributed delay terms of the epidemic model (3.1)–(3.7).

To complete the model formulation, applying some ideas in ^[22] to Model-Assumption 2.3, we take the function

where $\phi_{i}$ is the growth rate of the information and $\hat{\phi}_{i}$ is the saturation constant owing to the $i^{th}$ class $i\in \{I, T, A\}$.

Observe from Eq (3.10) that for $I(t), T(t), A(t)\in \mathbb{R}_{+}$,

where $\hat{\phi}_{min} = \min{(\phi_{k})}, k\in \{I, T, A\}$, and $\phi_{max} = \max{(\phi_{k})}, k\in \{I, T, A\}$.

The form for $F_{Z}$ in Eqs (3.10) and (3.11) signifies that the growth rate of the information in the population saturates as the infected population increases. Furthermore, it is assumed that there is relatively more preventive information related to getting infected (becoming ($I$) state) and getting treatment (becoming $T$ state) than progressing from HIV infectious ($I$) state into full-blown AIDS ($A$). That is,

Note that the family of epidemic models in Eqs (3.1)–(3.7) contains an interesting sub-family of models describing the dynamics of HIV/AIDS in the population, with information intervention, but in the absence of treatment (ART). This sub-family of models is obtained from Eqs (3.1)–(3.7) easily by deleting all treatment related parameters, leading to the system

with initial conditions in Eqs (3.8) and (3.9). Observe that the model (3.13) is a generalization of the model by Joshi et. al ^[4], that would be used to investigate the impacts of information intervention in changing human behavior in HIV/AIDS epidemics, whenever ART treatment is not available.

4.   Derivation of the stochastic HIV/AIDS model

From Assumptions (D) and (G) in Model-Assumptions 2.1–2.8, it is assumed there are noises in the HIV/AIDS dynamics from the random supply of ODAs and random poverty levels/living standards in the community reflected by life expectancy. That is, the proportions per unit time $\bar{\varepsilon}_{j} = 1-\varepsilon_{j}$ and $\varepsilon_{j}, j \in I(0, n)$ are random variables, and likewise the natural death rates $\mu_{k}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$ are random variables per unit time. Denote these random variables, respectively, by $\tilde{\bar{\varepsilon}}_{j}$, $\tilde{\varepsilon_{j}}, j \in I(0, n)$ and $\tilde{\mu_{k}}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$. We represent the environmental variabilities as independent white noise processes applying similar techniques in the earlier studies ^[27,28].

For $t\geq t_{0}$, let $(\Omega, \mathfrak{F}, P)$ be a complete probability space, and $\mathfrak{F}_{t}$ be a filtration (that is, sub $\sigma$- algebra $\mathfrak{F}_{t}$ that satisfies the following: given $t_{1}\leq t_{2} \Rightarrow \mathfrak{F}_{t_{1}}\subset \mathfrak{F}_{t_{2}}; E\in \mathfrak{F}_{t}$ and $P(E) = 0 \Rightarrow E\in \mathfrak{F}_{0}$). The variabilities in $\tilde{\bar{\varepsilon}}_{j}$, $\tilde{\varepsilon_{j}}, j \in I(0, n)$ and $\tilde{\mu_{k}}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$ in any small time interval of length $dt$ are represented by independent white noise processes as follows:

and the $w_k(t)$'s are the normalized Wiener processes for the $k^{th}$ state at time $t$ ($k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$), with the following properties: $w_{k}(0) = 0, E(w_{k}(t)) = 0, Var(w_{k}(t)) = t$. {Note from Eqs (4.1) and (4.2) that the random variables $\tilde{\bar{\varepsilon}}_{j}\in \mathbb{R}$, $\tilde{\varepsilon}_{j}\in \mathbb{R}, j \in I(0, n)$ and their means $\mathbb{E}[\tilde{\bar{\varepsilon}}_{j}] = \bar{\varepsilon}_{j} = 1-\varepsilon_{j}\in (0, 1)$ and $\mathbb{E}[\tilde{\varepsilon}_{j}] = \varepsilon_{j}\in (0, 1)$.

Also, for each $k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$, it is easy to see from Eqs (4.1) and (4.2) that $\mathbb{V}ar[\tilde{\mu_{k}}dt] = \sigma^{2}_{k}dt$, where $\sigma_k^{2}$ is the intensity of the noise in the natural death rate of the $k^{th}$ state; $\mathbb{V}ar[(1-\tilde{\varepsilon_{j}})dt] = \mathbb{V}ar[\tilde{\varepsilon_{j}}dt] = \sigma^2_{\varepsilon_j}dt\equiv \sigma^2_{\bar{\varepsilon}_j}dt$, where $\sigma^2_{\varepsilon_j}\equiv \sigma^2_{\bar{\varepsilon}_j}$ is the intensity of the noise in the random variable $\tilde{\varepsilon_j}, j\in I(0, n)$. Note, $\sigma_{\varepsilon_j}$ and $\sigma_{\bar{\varepsilon}_j}$ are identical, however, the distinct notations are utilized to emphasize the origins of the noises. The reader should see 10 for further discussion of the sub-models and white noise processes in Eqs (4.1) and (4.2).

Substituting Eqs (4.1), (4.2) and (10.4) into the deterministic systems (3.1)–(3.7) leads to the following generalized system of Ito-Doob stochastic differential equations.

and

where the initial conditions are given in the following: define $h = \max{(h_{1}, h_{2})}$,

and $\mathcal{C}((-\infty, t_{0}], \mathbb{R}_{+})$ is a Banach space of continuous functions endowed with the uniform norm

Furthermore, the random continuous functions $\vartheta_{j}(t), \varphi_{k}(t) > 0, \forall j\in I(0, n), k\in I(1, 5)$, are $\mathfrak{F}_{0}-measurable$, or independent of $w(t)$ for all $t\geq t_{0}$.

Observe that for $F_{Z}$ in Eqs (3.10) with (3.12), the Eqs (4.7) and (4.8) decouple from the stochastic systems (4.3)–(4.9). The following vectors defined below, will be used.

5.   Model validation results

Applying standard techniques in ^[27,28,33], it is proven that there exists a unique positive stochastic solution process $\{Y(t), t\geq t_{0}\}$ that satisfies Eqs (4.3)–(4.9).

It becomes apparent that the existence and behavior of the paths of the solution process for Eqs (4.3)–(4.9) depend on the sources and magnitudes of the intensities $\sigma_{k}, k \in\left\{S_{j}, I, T, A, R\right\}, j\in I(0, n)$ of the noises in the system. The results are given below.

Theorem 5.1. Define a closed ball in $\mathbb{R}_{+}^{n+6}$, centered at the origin and radius $r > 0$ in Eq (5.2) as follows:

where

Given the initial conditions (4.10) and (4.11), there exists a unique positive stochastic solution process (SSP) $\{Y(t, w), t\geq t_{0}, w\in \Omega\}$ for the systems (4.3)–(4.9). That is, $S_{j}(t, w) > 0, I(t, w) > 0, T(t, w) > 0, A(t, w) > 0, R(t, w) > 0$, and $Z(t, w) > 0, t\geq t_{0}, w\in \Omega$, almost surely. Moreover, the following hold:

a. If the intensities in Eqs (4.3)–(4.9) satisfy $\sigma_{\bar{\varepsilon}_{j}}\equiv\sigma_{\varepsilon_{j}} = 0$, $\sigma_{S_{j}} = 0, \forall j\in I(0, n)$ and $\sigma_{k} = 0, \forall k\in \{I, T, A, R\}$, then the solution of the ensuing deterministic systems (3.1)–(3.7) lies in a positive self-invariant space $\mathbb{\bar{B}}_{\mathbb{R}^{n+6}}(0, r)$ in (5.1).

b. If the intensities satisfy $\sigma_{\bar{\varepsilon}_{j}}\equiv\sigma_{\varepsilon_{j}} > 0$, $\sigma_{S_{j}} = 0, \forall j\in I(0, n)$ and $\sigma_{k} = 0, \forall k\in \{I, T, A, R\}$, then the solution of the ensuing stochastic systems (4.3)–(4.9) almost surely lies in a positive self-invariant space $\mathbb{\bar{B}}_{\mathbb{R}^{n+6}}(0, r)$ in Eq (5.1).

c. If at least one of $\sigma_{k} > 0$, $k\in \{S_{j}, I, T, A, R\}$, and $j\in I(0, n)$, then the solution of the ensuing stochastic systems (4.3)–(4.9), almost surely lies in the unbounded space $\mathbb{R}^{n+6}_{+}$, for all time $t\geq t_{0}$, regardless whether $\sigma_{\bar{\varepsilon}_{j}}\equiv\sigma_{\varepsilon_{j}}\geq 0$.

Proof. See 11 for the complete proof of Theorem 5.1. In addition, an interpretation of the results is given in Remark 11.1.

6.   Stochastic Stability of the disease-free steady state

Denote by $E = X^{*} = (S^{*}_{0}, \ldots, S^{*}_{n}, I^{*}, T^{*}, Z^{*})^{T}$ the steady states of the systems (3.1)–(3.7) and (4.3)–(4.9). Similarly to Theorem 5.1, the existence of a disease-free equilibrium (DFE) of Eqs (4.3)–(4.9) depends on the intensities of the noises in the system.

Theorem 6.1. Let Theorem 5.1 hold. The following are true about the DFE $E_{0}$:

(a.) When Theorem 5.1(a.) holds, the systems (4.3)–(4.9) and (3.1)–(3.7) are equivalent and the DFE is given by

(b.) When Theorem 5.1(b.) holds, the stochastic systems (4.3)–(4.9) has a DFE given by (6.1).

(c.) When part of Theorem 5.1(c.) holds, i.e. $\sigma_{S_{0}} = 0$, $\sigma_{S_{j}} > 0, \forall j\in I(1, n)$ and $\sigma_{k} > 0, \forall k\in \{I, T, A, R\}$, and $\sigma_{\bar{\varepsilon}_{j}}\equiv\sigma_{\varepsilon_{j}} > 0$, then the DFE exists and is given by Eq (6.1).

(d.) When $\sigma_{S_{0}} > 0$ in Eqs (4.3)–(4.9), the systems (4.3)–(4.9) no longer has a DFE.

Proof. The proof of Theorem 6.1 is given in 12. In addition, the interpretation of the theorem in relation to the disease dynamics is given in Remark 12.1.

The next result presents the basic reproduction numbers (BRN $\mathfrak{R}_{0}$) of the two families of deterministic models (3.1)–(3.7) and (3.13).

Theorem 6.2. Define

(a.) The BRN of the deterministic models (3.1)–(3.7) is given by

(b.) The BRN of the deterministic model (3.13) is given by

Proof. The proof is easily obtained by applying the method of next generation matrix.

Remark 6.1.

(1.) The BRN $\mathfrak{R}_{0}$ in (6.4) is interpreted in the following. Indeed, given one infectious individual placed in the disease free population $E_{0}$, the term $\varepsilon_{0}R_{1}E(e^{-\mu_{I}\tau_{1}})$ represents all newly infected individuals who fail to receive treatment and just turn into full-blown AIDS after the period $\tau_{1}$; the term $(1-\varepsilon_{0})R_{1}E(e^{-\mu_{I}\tau_{2}})-\alpha_{TI}(1-\varepsilon_{0})R_{1}E(e^{-\mu_{I}\tau_{2}}) \frac{1}{(\alpha_{TI}+\mu_{T}+d_{T}+\alpha_{TA})}$ represents the net number of newly infected individuals who remain in treatment (note that $\alpha_{TI}(1-\varepsilon_{0})R_{1}E(e^{-\mu_{I}\tau_{2}})\frac{1}{(\alpha_{TI}+\mu_{T}+d_{T}+\alpha_{TA})}$ represents all newly infected individuals who have either stopped treatment and currently returning to the infectious state, or those in whom treatment fails, and they are currently changing into full-blown AIDS); the term $\frac{1}{G^{'}(0)}$ is the effect of nonlinearity in the incidence rate of the disease. Therefore, $P(1)$ is approximately the probability of finding an infectious individual who would either become infectious, receive treatment and remain treated, or progress without treatment into full-blown AIDS.

Thus, the basic reproduction number $R_{0}$ in (6.4) is the average number of secondary infections in the complete disease free population $E_{0}$, that would proceed to the treatment class or to the full-blown AIDS class, over the average lifespan of $\frac{1}{\mu_{I}+d_{I}}$ of an infectious person in the population. Note that the BRN $\mathfrak{R}_{0}$ in (6.5) is interpreted similarly for the system (3.13) without ART treatment.

(2.) Observe from (6.4) and (6.5) that $R_{0} < R^{noTr}_{0}$. Which implies that $R_{0}\rightarrow 0$ faster than $R^{noTr}_{0}\rightarrow 0$, and $R_{0}\rightarrow 1$ at a slower rate than $R^{noTr}_{0}\rightarrow 1$. This observation suggests that the disease is more easily eradicated when there is ART treatment and behavioral modification via information intervention, than when there is only behavioral change without ART treatment.

The following lemmas will be used to show the stochastic stability results for the DFE $E_{0} = X_{0}^{*} = (S^{*}_{0}, \ldots, S^{*}_{n}, I^{*}, T^{*}, Z^{*}) = (\frac{B}{\mu_{S_{0}}}, 0\ldots, 0, 0, 0, 0)$ of the system (4.3)-(4.9). The decoupled system with positive solution $X(t), t\geq t_{0}$ in (4.12) is utilized. Observe that the HIV positive states in $X(t)$ are $I$ and $T$. Also observe from Assumptions 2.1 & 2.2 that for each $j\in I(0, n)$, $H_{j}$ and $G_{j}$ are continuous and bounded over their domains. Moreover, denote by

Clearly, from Assumption 2.1 (A5) and Assumption 2.2(A5), it is easy to see that for any $t\geq t_{0}$,

Lemma 6.1. Let Theorem 5.1 hold, and define the $C^{2, 1}$-function $V:\mathbb{R}_{+}^{n+4}\times \mathbb{R}_{+}\rightarrow \mathbb{R}_{+}$, where

and

Also, let $\phi_{j}, j\in I(0, n)$ and $\varphi_{k}, k\in \{I, T, Z\}$ be defined as follows:

and

The differential operator $dV$ applied to $V(t)$ with respect to the stochastic system (4.3)-(4.9) satisfies the following:

where $\overrightarrow{w(t)} = (w_{\varepsilon_{0}}, w_{\varepsilon_{1}}, \ldots, w_{\varepsilon_{n}}, w_{S_{0}}, \ldots, w_{S_{n}}, w_{I}, w_{T})^{T}$,

In addition,

where $\mathbb{E}_{\tau_{1}}$ is the expectation operator with respect to the random variable $\tau_{1}$.

Proof. The complete proof of Lemma 6.1 is given in 13. Furthermore, an estimate for the term $\sum^{n}_{j = 0}\left(2\varepsilon_{j}\beta_{j}\frac{1}{\lambda(\mu)} +\sigma^{2}_{\varepsilon_{j}}\beta^{2}_{j}\right)\mathbb{E}_{\tau_{1}}\left[e^{-2\mu_{I}\tau_{1}}\right] S^{2}_{j}(t)G^{2}_{j}(i(t))$ in (6.18) used in the proof of Lemma 6.2 is given in Note 13.1.

In the next subsection, we examine stochastic stability of the DFE $E_{0}$ in $D(\infty)$.

6.1.   Stochastic Stability in probability of the DFE in $D(\infty)$

The next result presents conditions for stochastic stability results in probability for the DFE $E_{0}$, when Theorem 5.1(b) and Theorem 6.1(b) hold and the basic reproduction number $R_{0} < 1$.

Lemma 6.2. Let the assumptions for Lemma 6.1 be satisfied and let $R_{0}$ denote the BRN in (6.4). Also, let Theorem 6.1(b.) hold, i.e. the intensities of (4.3)-(4.9) satisfy $\sigma_{\bar{\varepsilon}_{j}}\equiv\sigma_{\varepsilon_{j}} > 0$, $\sigma_{S_{j}} = 0, \forall j\in I(0, n)$ and $\sigma_{k} = 0, \forall k\in \{I, T, A, R\}$, then the DFE, given by (6.1), for the system (4.3)-(4.9) exists in the positive-self invariant space $D(\infty)$. Moreover, in $D(\infty)$ let $V(t)$ be as defined in (6.8) and define the functional

Also define the following

and

Furthermore, let

and for each $j\in I(1, n)$,

Define

and

It follows that in $D(\infty)$, there exists positive constants $\phi_{j} > 0, j\in I(0, n)$ and $\varphi_{k} > 0, k\in \{I, T, Z\}$, such that the drift part of $d(V(t)+V_{5}(t))$ satisfies

whenever $R_{0, \sigma_{\varepsilon}} < 1$, $R_{1}(\mu) < 1+(1-R_{0, \sigma_{\varepsilon}})$, $\max{(U_{j})}_{j\in I(0, n)} < 1$, and $\max{(V_{0}, W_{0})} < 1$.

Proof. The complete proof of Lemma 6.2 is given in 14.

Theorem 6.3. Suppose the conditions of Lemma 6.2 are satisfied, then it follows that the DFE $E_{0}$ of the stochastic systems (4.3)–(4.9) is stochastically asymptotically stable in the large in $D(\infty)$. That is, the disease is eliminated from the population asymptotically, whenever the conditions $R_{0, \sigma_{\varepsilon}} < 1$, $R_{1}(\mu) < 1+(1-R_{0, \sigma_{\varepsilon}})$, $\max{(U_{j})}_{j\in I(0, n)} < 1$, and $\max{(V_{0}, W_{0})} < 1$ in Lemma 2 are satisfied.

Proof. When the conditions of Lemma 6.2 are satisfied, then from Eq (6.26), the drift part of $d(V(t)+V_{5}(t))$ is negative definite. Hence, applying Theorem 5 in ^[28], the result follows immediately ^[28,33,34].

6.2.   Persistence of the susceptible state $S(t)=\sum^{n}_{j=0}S_{j}(t)$ in $D(\infty)$

Since from Eq (6.1) the DFE for the susceptible states $S_{j}, j\in I(1, n)$ are $0$, and only the DFE of $S_{0}$ is positive, it is necessary to emphasize that the paths of the combined susceptible population $S(t) = \sum^{n}_{j = 0}S_{j}(t)$ will persist for all time, and never completely die out over time.

Indeed, it is shown that $S(t)$ will persist and approach the DFE $S^{*}_{0} = \frac{B}{\mu_{S_{0}}}$, whenever the stochastic stability conditions for the DFE in Lemma 6.2 and Theorem 6.3 hold.

Recall ^[27], the following definition of the persistence of a species denoted by the process $y(t), t\geq t_{0}$ in a stochastic dynamic system:

Definition 6.1.

(1.) $y(t)$ is said to be stable in the mean if $\lim_{t\rightarrow \infty}{\frac{1}{t}\int^{t}_{t_{0}}y(s)}ds = c > 0$, a.s.

(2.) $y(t)$ is said to be weakly persistent in the mean, if $\limsup_{t\rightarrow \infty}{\frac{1}{t}\int^{t}_{t_{0}}y(s)}ds > 0$, a.s., and strongly persistent in the mean, if $\liminf_{t\rightarrow \infty}{\frac{1}{t}\int^{t}_{t_{0}}y(s)}ds > 0$, a.s.

Theorem 6.4. Suppose the conditions of Lemma 6.2 and Theorem 6.3 are satisfied, then for $S(t) = \sum^{n}_{j = 0}S_{j}(t)$

That is, all paths of the combined susceptible state $S(t)$ are persistent and stable in the mean, almost surely, and converge to the DFE $S^{*}_{0} = \frac{B}{\mu_S{0}}$.

Proof. The complete proof of Theorem 6.4 is given in 15.

6.3.   Discussions of the stability results and sensitivity analysis of the BRN $\mathfrak{R}_{0}$

Remark 6.2. Interpreting the stochastic stability results in $D(\infty)$

(a.) Probabilistically, Theorem 6.3 and Lemma 6.2 signify that when the conditions of Lemma 6.2 hold, then all paths of the SSP $\{Y(t), t\geq t_{0}\}$ of Eqs (4.3)–(4.9) originate, oscillate and remain bounded in the closed ball $D(\infty)$, a.s. Moreover, there exists a DFE $E_{0}$ for the system, and any path that starts near $E_{0}$ has a high chance to continue oscillating near $E_{0}$, and over time, the path certainly converges to $E_{0}$.

Biologically, these results suggest that the maximum spread of the disease can not exceed the bounds of $D(\infty)$, regardless of the intensity of the fluctuations in the supply of ODAs. Furthermore, the disease can be eliminated, whenever the conditions in Theorem 6.3 and Lemma 6.2 hold.

Also, note from Theorem 6.3 and Lemma 6.2 that the new BRN $\mathfrak{R}_{0}$ that is modified by the noise in the system is $R_{0, \sigma_{\varepsilon}}$ in Eq (6.20). And when $R_{0, \sigma_{\varepsilon}} < 1$, and the other threshold conditions $R_{1}(\mu) < 1+(1-R_{0, \sigma_{\varepsilon}})$, $\max{(U_{j})}_{j\in I(0, n)} < 1$, and $\max{(V_{0}, W_{0})} < 1$ are satisfied, the disease is eliminated.

Remark 6.3. Sensitivity of the BRN $\mathfrak{R}_{0}$ to the supply of ODAs, delayed ART treatment and relapse from treatment

(b.) Recall assumptions (A)-(G) in Model-Assumptions 2.1–2.8, the effects of the random supply of ODAs, poverty levels, IECs and delayed ART treatment are assessed via the parameters $\varepsilon_{j}, \bar{\varepsilon}_{j}, j\in I(0, n)$, $\mu_{k}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$, $\gamma_{j}\equiv \gamma_{S_{0}S_{j}}, j\in I(1, n)$, and $\tau_{2}$, respectively. The effects of continuous changes in these parameters on the new BRN in Eq (6.20) are examined below.

First, observe from Eqs (6.3) and (6.4) that the magnitude of the probability value $P(1)$ determines the magnitude of the BRN $R_{0}$ in absence of noise in the system. Furthermore,

where

It is easy to see that $\frac{\partial P(1)}{\partial \bar{\varepsilon}_{0}} < 0$, if

and $\frac{\partial P(1)}{\partial \bar{\varepsilon}_{0}} > 0$, otherwise.

In other words, when Eq (6.30) holds, then $P(1)$ continuously decreases to zero (i.e. $0 < P(1) < < 1$) as $\bar{\varepsilon}_{0}\rightarrow 1$, and consequently the BRN $R_{0}$ continuously decreases to zero (i.e. $0 < R_{0} < < 1$), as $\bar{\varepsilon}_{0}\rightarrow 1$, and vice versa. This fact is exhibited in an Figure 7.

Note from Eq (6.30), the left-hand side (LHS) given by $E(e^{-\mu_{I}\tau_{2}})\left(1-\alpha_{TI} \frac{1}{(\alpha_{TI}+\mu_{T}+d_{T}+\alpha_{TA})}\right)$, is the expected conditional probability that a newly exposed HIV person, survives natural death (with probability $E(e^{-\mu_{I}\tau_{2}})$) and remains in ART treatment without relapsing to infectiousness (with probability $\left(1-\alpha_{TI} \frac{1}{(\alpha_{TI}+\mu_{T}+d_{T}+\alpha_{TA})}\right)$), given that the infected person is tested and begins ART treatment after $\tau_{2}$ time units. In other words, the LHS is the expected probability that an individual infected with HIV will remain healthy, and alive during the ART treatment.

The term, $E(e^{-\mu_{I}\tau_{1}})$, on the right-hand side (RHS) of Eq (6.30) is clearly the expected conditional probability that an exposed person will progress into full-blown AIDS after $\tau_{1}$ time units, given that the person fails to be tested and get ART treatment.

Thus, Eqs (6.28)–(6.30) suggest that when ODAs are readily available, and more people tend to get ART treatment (i.e. $\bar{\varepsilon}_{0}\rightarrow 1$), then the BRN $R_{0} < < 1$ and the disease is more easily eliminated.

Also, observe from Eq (6.30) that when either $\tau_{2}$ is small (i.e. $\tau_{2}\rightarrow 0$), or $\alpha_{TI}$ is small (i.e. $0 < \alpha_{TI} < < 1$), then $LHS > > RHS$, assuming all other constants are fixed. This implies that when either $\tau_{2}\rightarrow 0$, or $0 < \alpha_{TI} < < 1$, then $P(1)\rightarrow 0$, and consequently $R_{0} < < 1$. In other words, when newly infected people get tested soon after infection, and begin early ART treatment (i.e. $\tau_{2}\rightarrow 0$), then increasing the supply of ODAs (i.e. $\bar{\varepsilon}_{0}\rightarrow 1$) will result in more people getting tested and paying for ART treatment, which makes the disease eradication process easier. This fact is illustrates in the example in Figure 8.

The alternative when $\alpha_{TI}$ is small and $R_{0} < < 1$ signifies that when fewer number of people who are receiving ART treatment relapse from treatment into the infectious state, then increasing the supply of ODAs (i.e. $\bar{\varepsilon}_{0}\rightarrow 1$), so that more people get tested and pay for ART treatment, will make the disease eradication process easier.

Remark 6.4. Sensitivity of $R_{0}$ to the distribution of the delays:

(c.) Some information is known about the delay, $\tau_{1}$, after infection until full-blown AIDS occurs ^[1,2]. In fact, over 2 to 15 years, untreated HIV individuals develop full-blown AIDS ^[1,2]. Using Eq (6.30), more specific sufficient conditions for $\tau_{2}$ leading to $R_{0} < < 1$ are derived, whenever the distributions of the random delays $\tau_{1}$ and $\tau_{2}$ are specified. Recall, the delays $\tau_{1}$ and $\tau_{2}$ are distributed with densities $f_{\tau_{1}}, f_{\tau_{2}}$. One example is considered below, whenere the delays $\tau_{1}$ and $\tau_{2}$ have exponential distribution.

Delays $\tau_{1}$ and $\tau_{2}$ have exponential distributions:

Assuming that those who are newly infected proceed into the ART class or AIDS class independently at constant rates $\theta_{2}$ and $\theta_{1}$, respectively, then the time until treatment begins $\tau_{2}$, and the time until full-blown AIDS develops $\tau_{1}$, follow exponential distributions with means $E[\tau_{2}] = \frac{1}{\theta_{2}}$ time units, and $E[\tau_{1}] = \frac{1}{\theta_{1}}$ time units, respectively. In this scenario, two independent Poisson processes describe the number of people getting ART treatment and developing full-blown AIDS over time.

It is easy to see that Eq (6.30) becomes

where given that a newly infected person will begin ART treatment, $\frac{\theta_{2} \frac{1}{\mu} }{\theta_{2} \frac{1}{\mu}+ \mu \frac{1}{\mu}}$ is the conditional fraction of the newly infected population that successfully gets into ART treatment. This implies that the fraction $1-\frac{\theta_{2} \frac{1}{\mu} }{\theta_{2} \frac{1}{\mu}+ \mu \frac{1}{\mu}}$ will die naturally.

Also, $\frac{\theta_{1} \frac{1}{\mu} }{\theta_{1} \frac{1}{\mu}+ \mu \frac{1}{\mu}}$ is the conditional fraction of those who will turn into full-blown AIDS, given that they do not receive ART, and do not die naturally. Therefore, the LHS of Eq (6.31) is the fraction of the newly infected population that survives natural death, gets ART treatment, and does not relapse from treatment. The RHS is the fraction of the newly infected population that fails to get ART treatment, and develops full-blown AIDS.

This special example confirms the observation that, when more newly infected people are likely to get tested, begin early ART treatment and do not relapse from treatment, then increasing the supply of ODAs (i.e. $\bar{\varepsilon}_{0}\rightarrow 1$) so that more people are able to afford ART treatment, makes it easier to control the disease.

Remark 6.5. Sensitivity of the BRN $\mathfrak{R}_{0}$ to the delays when the supply of ODAs is fixed

(d) Now assuming the supply of ODAs is fixed (i.e. $\bar{\varepsilon}_{0}$ and ${\varepsilon}_{0}$ are fixed), to examine the direct effects of $\tau_{1}$ or $\tau_{2}$ on the BRN $R_{0}$, it is easy to see from (6.4)-(6.3) that $\frac{\partial R_{0} }{\partial E[e^{-\mu_{I}\tau_{1}}]} < 0$ and $\frac{\partial R_{0} }{\partial E[e^{-\mu_{I}\tau_{2}}]} < 0$. That is, the BRN $\mathfrak{R}_{0}$ decreases as $\tau_{1}$ and $\tau_{2}$ decrease.

One plausible reason why the BRN $R_{0}$ decreases with decrease in $\tau_{1}$, the time until full-blown AIDS occurs for those who are not receiving ART, is because the AIDS population in this model is not an active spreader of the disease. Thus, when the incubation period of HIV without treatment $\tau_{1}$ is shorter, then more people not receiving treatment quickly proceed into the non-spreading infected state $A(t)$.

Similarly, the decrease in the BRN $R_{0}$ with $\tau_{2}$ is plausible because when all newly infected persons who are able to afford ART treatment begin early treatment, then there will be reduced spread of the disease, since the model does not allow spread between the treated state $T(t)$, and other states in the population. An example of this scenario is exhibited in Figure 8.

Remark 6.6. Sensitivity of the BRN $\mathfrak{R}_{0}$ to the noises in the supply of ODAs

(e) Introducing noise into the disease dynamics by making the supply of ODAs random over time, it follows that the random effect of the supply of the ODAs inflates the BRN $R_{0}$, leading to the new BRN $R_{0, \sigma_{\varepsilon}}$ in (6.20), where the term $\frac{r^{2} \mathbb{E}_{\tau_{1}}\left[e^{-2 \mu_{I} \tau_{1}}\right]}{\left(1-\frac{1}{2} \lambda_{2}(\mu)\right)\left(\mu_{I}+d_{I}\right)} \sum_{j = 0}^{n} \beta_{j}^{2} \sigma_{\varepsilon_{j}}^{2}$ is the increment that depends on the intensities of the noises, whenever the ODA parameters $\varepsilon_{j}, \bar{\varepsilon}_{j}, j\in I(0, n)$ become random variables over time.

Clearly, the condition $R_{0, \sigma_{\varepsilon}} < 1$ is constrained by the intensity of the noise in the system. Since the modified BRN $R_{0, \sigma_{\varepsilon}}$ decreases, as $\sigma_{\varepsilon_{j}}\rightarrow 0$, and $R_{0, \sigma_{\varepsilon}}$ increases, otherwise, this suggests that a condition for the disease dynamics, where there is a constant significant supply of ODAs (i.e. the expected values $E[\bar{\varepsilon}_{j}]$ are large, and $\sigma_{\varepsilon_{j}}\rightarrow 0$) is most favorable to the elimination of the disease from the system, and vice versa. An example for this scenario is exhibited in Figure 8.

7.   Numerical simulation results for HIV/AIDS epidemic in Uganda

(a.) HIV/AIDS data:

In this section (1) ideas in Joshi et al. ^[4], (2) data reported for Uganda in the sources ^{[1,4,35,36,37]}, over time intervals between 1992 through 2018, and (3) data from clinical trials ^[38,39,40] on ART treatment, are used to design parameter estimates for the HIV/AIDS models (4.3)–(4.9), for a hypothetical population with comparable characteristics as Uganda over 1992-2018.

Note that access to ART treatment, and the use of PrEP in Uganda, are relatively recent, in fact, the Ugandan government launched universal access to ART treatment in 2004 ^[41]. Thus, given the apparent time and missing information gap on ART treatment and PrEP in Uganda between 1992 to 2004, we resort to theoretical arguments in many instances, to select plausible estimates of some parameters of the models (4.3)–(4.9).

Therefore, the primary objective of this section is mainly to numerically explore the effects of ART treatment, the supply of ODAs, early HIV testing-and-onset of ART treatment and noise in the supply of ODAs, on the prevalence of HIV/AIDS in the population.

(b.) Initial conditions:

Similarly to the Uganda UN data in 1991 ^[4], we assume the initial population consists of 18.38 million people. Furthermore, according to the UN population report for 1991 on Uganda ^[4], the estimated percentage of adults between the ages of 15 to 59 in the population was about 32%. Using this estimate, it is assumed in this example that the total sexually active adult population is 32%, consisting of individuals who are either vulnerable, removed or infected with HIV. Thus, the 5, 899, 980 million adults consist of the states $S_{0}$, $S_{1}$, $S_{2}$, $I$, $T$, $A$ and $R$. In ^[4], the prevalence rate of HIV assumed at that time was 15% in Uganda. This assumption led to 884, 997 initially infected ($I(0)$) and 5, 014, 980 initially susceptible $S(0)$. Since no treatment ($T$) and full-blown AIDS ($A$) states were considered in ^[4], we make modifications. In the 5, 899, 980 initial adult population, 15% is assumed initially infected or removed, i.e. about 884, 997 adults in $I, T, A$ or $R$ states. Moreover, there is a 12% prevalence rate of HIV (i.e. $I(0)$ is about 707, 997 adults) and 2% treatment rate (i.e. $T(0)$ is about 117, 999 adults); 1% of the initial adult population is assumed to be fully trained (by organizations in the community initially involved in IEC) in the proper application of strong HIV preventive measures such as the daily use of PrEP ^[24,25,26], and as a result of their training, they effectively adhere to this practice, and are completely removed from HIV transmission (i.e. $R(0)$ is about 58, 999 adults); 0% of the initial adult population is assumed to be in the AIDS state (i.e. $A(0) = 0$). Thus, in the remaining 85% of the adult population, all susceptibles (about 5014983 adults), we assume 10% are $S_1(0)$ (about 501, 498 adults), susceptibles adopting safe practices such as abstinence and being in a mutually monogamous relationship; 10% are $S_2(0)$ (about 501, 498 adults), susceptibles adopting all other medically advised practices such as condom use; the remaining 80% are susceptibles $S_{0}$ (about 4, 011, 986 adults), not practicing any safe sex measures.

Also, similarly to ^[4], it is assumed that there are initially about 240 organizations (20%) involved in the IEC. That is, $Z(0) = 0.2$. Applying the percentages above and scaling the initial values (by 1 million), lead to the initial conditions (7.1).

(c.) Information rates:

Similarly to ^[4], the size of information in the population $Z(t)$ is based on the number of organizations giving out information on HIV. Furthermore, the breakdown into the two states $S_{1}(t)$ and $S_{2}(t)$ at any time depends on the proportion of organizations giving information about the respective behaviors.

Therefore, as in ^[4], it is assumed 10% per unit time of $S_0(t)$ interact with $Z(t)$ and become $S_1(t)$, and 80% per unit time interact with $Z(t)$ and become $S_2(t)$. Also, another 10% per unit time of $S_0(t)$ interact with $Z(t)$ and become $R(t)$. That is, $\gamma_0 = 0.1$ per year, $\gamma_1 = 0.1$ per year and $\gamma_2 = 0.8$ per year.

(d.) Infection rates:

Some studies ^[6] argue that the effects of behavioral changes such as abstinence, delayed sexual initiation, mutual monogamy and correct and consistent condom use on reducing HIV prevalence rates are complex to understand and debatable. However, UNAIDS ^[37], asserts that condom use effectively reduces the risk of transmission of HIV by the range of 69-94%. Siding with UNAIDS, we utilize some information from ^[4] to estimate $\beta_{j}, j = 0, 1, 2$. That is, it is assumed that $\beta_{0}\in [0.01, 0.1]$, and that $\beta_{1}$ and $\beta_{2}$ are proportional to $\beta_{0}$. Moreover, $\beta_{1} < < \beta_{0}$ and $\beta_{1} < < \beta_{2} < \beta_{0}$, where $\beta_{0}$, $\beta_{1}$ and $\beta_{2}$ are respective disease transmission rates in the susceptible class $S_{0}$ (those not practicing any safe sex measures), in $S_{1}$ (those practicing abstinence and be faithful), and in $S_{2}$ (those practicing condom use). Based on the above relationships, it is assumed that $\beta_1 = 0.05 \beta_0$, and $\beta_2 = 0.4\beta_0$.

The Uganda UNAIDS ^[35] HIV incidence per 1000 population time series data for adults (15-49 years) in Figure 3 is used to find an estimate for $\beta_{0}$ over time. Note that the values for the HIV incidence per person obtained from UNAIDS ^[35] approximate the true values of $\beta_{0}$. Also, in real life, the value of $\beta_{0}$ typically changes over time. Thus, the estimate for $\beta_{0}$ is obtained over time $t\geq 0$ by fitting the multiple linear regression model $\beta_{0}(t) = b_{0}+b_{1}t+b_{2}t^{2}+b_{3}t^{3}+b_{4}t^{4}+Error$ to the data in Figure 3. The least squares fit is given by

with p-values for the significance of regression coefficients: $b_{1}(p-value = 1.73e-14 < 5\%)$, $b_{2}(p-value = 1.01e-09 < 5\%)$, $b_{3}(p-value = 6.98e-08 < 5\%)$, and $b_{4}(p-value = 2.96e-06 < 5\%)$, and multiple R-squared $R^{2} = 0.9964$. Moreover, the residual plots for the model are given in Figure 4, and the Prediction Multiple R-squared based on the PRESS statistic is $R^{2}_{Prediction} = 0.9937$. Clearly, the regression model has a good fit from the residual plots in Figure 4, and also has a high predictive capacity from the Prediction Multiple R-squared $R^{2}_{Prediction} = 0.9937$.

It is also easy to see that $\hat{\beta}_1(t) = 0.05 \hat{\beta}_0 (t)$ and $\hat{\beta}_2(t) = 0.4\hat{\beta}_0 (t)$.

(e.)Influx of new adults, natural and disease related death rates:

Using UN data for Uganda ^[36], the natural death rate per year is estimated using the mean lifespan of individuals in Uganda over the years $1997$-$2005$. It is easy to see that the natural deathrate for a mean lifespan of $47.51667$ years is $\mu_{k} = \mu = 1/47.51667 = 0.02104525$ per year, $k\in \{S_{j}, I, T, A, R\}, j\in I(0, 2)$.

Similarly, the UNAIDS Uganda report ^[37], indicates that HIV related death rates per year were $17.2\%$ in $1997$, $13.36\%$ in $1999$, $14\%$ in $2001$, $14.42\%$ in $2003$. This leads to the mean death rate $d_I = 0.1474$ per year.

According to Kasamba et. al. ^[42], ART led to significant reduction of the mortality rate of HIV positive people in Uganda. It was exhibited (cf. ^[42]) that in the same HIV positive group with death rate of 116.4 deaths per 1, 000 persons per year prior to ART, a reduction by 25% mortality rate is obtained after ART, leading to 87.4 deaths per 1, 000 person per year. Thus, based on this discussion, we assume $d_T = 0.25 d_I = 0.03685$ per year.

According to ^[1], a white blood cell (CD4) count, lower than 200 cells/$mm^3$, is classified as full blown AIDS. Furthermore, if not treated immediately, the infected person cannot survive beyond 3 years. This suggests that the average lifetimes $\frac{1}{d_{A}} < < \frac{1}{d_{I}}$. Given the apparent challenge of tracking and recording the progression of HIV infection until full-blown AIDS develops, data for HIV related deaths are recorded in many cases without a clear distinction between deaths that occur when the white blood cell (CD4) count is either lower than or above 200 cells/$mm^3$ ^[37]. Thus, using the estimated mean death rate $d_I = 0.1474$ per year above, and the relationship $\frac{1}{d_{A}} < < \frac{1}{d_{I}}$, we assume that $d_A$ is twice larger than $d_I$. That is, $d_A = 2d_{I} = 0.2948$ per year.

In ^[4], data for new susceptible adults entering the population over years 1990 to 2005 is utilized to calculate the mean influx rate of $0.55$ per year of new susceptibles into the population. Using this estimate, we assume that $B = 0.55$ per year new "naive" adult susceptibles of type $S_{0}$, enter the population at any time.

(f.) Delays until treatment and AIDS:

As mentioned in the introduction, without proper ART, the infected individual progresses to AIDS in 2 to 15 years. That is, $\tau_{1}\in [2, 15]$ years. Thus, for the purpose of simulations, assume that $\tau_1 = 8$ years. Also, according to ^[1], it is suggested that within 2 to 12 weeks a newly infected person develops HIV antibodies, and over 6 months nearly every infected person can be diagnosed for HIV. Thus, early ART is assumed between $\tau_{2}\in [0.038, 0.5]$ years. Also, for various simulations presented below, $\tau_{2}$ will be selected in a time range until an individual develops full-blown AIDS, i.e. $\tau_{2}\in [0.038, 15]$ years.

(g.) Nonlinear incidence rates and nonlinear information rates:

Figure 5(a, b) shows the time series of Uganda yearly cumulative HIV population and new HIV infections, respectively, in adults (15–49 years) over the years 1990-2018 ^[35]. Clearly, Figure 5(a, b) shows that the yearly cumulative HIV population is increasing at a decreasing rate, and the new HIV cases in Figure 5(b) decrease with time. That is, the rate of new infections slowly decreases as the yearly cumulative HIV population increases. This observation suggests that a bilinear incidence rate $\beta S(t)G(I(t))$, where the bilinear incidence function $G(I(t)) = I(t)$ depicted in Figure 5(c) is not appropriate. Also, the standard incidence rate $\beta S(t)G(i(t))$, where the standard incidence function $G(i(t)) = \frac{I(t)}{N(t)}$ depicted in Figure 5(d) is not appropriate. The most appropriate incidence rate for the Uganda HIV epidemic (considering all three options) is the modified standard incidence rate $\beta S(t)G(i(t))$ depicted in Figure 5(d), where the modified standard incidence function saturates, for example, in Figure 5(d), $G(i(t)) = \frac{b i(t)}{1+c i(t)}, b = 0.05, c = 5$.

Thus, from the argument in Figure 5, the nonlinear functions in Assumption 2.1 and Assumption 2.2 are taken to be

where $0 \leq a_j, b_j \leq1, c_{j}, d_{j}\geq 1$ for $j = 0, 1, 2$ are considered.

For simulating system (4.3)-(4.9), it is assumed that $a_j = b_j = 0.1 \quad \forall j = 0, 1, 2$. Similarly, for the growth rate function given in (3.10), it is assumed that $\phi_I = 0.05, \phi_T = 0.1, \phi_A = 0.15, \hat{\phi_I} = 0.03, \hat{\phi_T} = 0.06, \hat{\phi_A} = 0.09$.

(h.) Withdrawal from ART treatment rates:

Clinical trials have been conducted to determine factors influencing withdrawal from ART treatment, and also to estimate the time until withdrawal, particularly for the first highly active antiretroviral(HAART) regimen ^{[38,39,43,44]}. Some of the factors whose effects on withdrawal have been quantified in these studies include: drug toxicity, failure of treatment, gender, non-compliance and specific regimens of HAART containing combinations of inhibitors namely- protease inhibitors(PI), non-nucleoside reverse transcriptase inhibitors(NNRTI), and nucleoside reverse inhibitors (NRTI) ^[38,39]. It is estimated that when no toxicity occurs due to a HAART regimen, then less than 10% naive patients withdraw from treatment within 1 year because of failure of treatment ^[39], while taking into account all factors affecting withdrawal above, the cumulative probability of withdrawing from HAART at a one year period after starting ART is 37.6%, with a 95% CI of (25.3%, 39.9%) ^[38]. Toxicity is the commonest reason for withdrawal or modification of ART. Indeed, more recently, Torres et. al. ^[43] observed in a clinical trial involving 1558 patients of various age groups, that toxicity was the culprit in 95% (228 out of 239 total withdrawal) of all cases that withdrew from ART. Azevedo et. al. ^[44] observed that the main adverse events leading to modification or withdrawal from ART in the first year of treatment, are dermatological, neuropsychiatric and gastrointestinal events. Moreover, the estimated median time between initiating ART and the first modification or withdrawal due to adverse events is 70.5 days (95% CI: 26-161 days).

In this example, we assume that withdrawal from ART is influenced by the combination of factors above, and not only toxicity. Hence, as in ^[38], we estimate the rate of withdrawal from ART to be $\alpha_{TI} = 0.376$.

(i.) Failure rate of ART treatment:

Clinical trials have also been conducted to determine predictors of first line failure of HAART ^[45], and the failure rate of HAART for individuals who began treatment in the advanced stages of their HIV infectiousness ^[40]. Indeed, in the retrospective study ^[45] involving 195 patients initiated on first line ART, 7.69% failure rate is observed, and the significant predictors of failure were BMI, CD4 count at ART onset, and opportunistic infection presence. Also, in the clinical trial ^[40], over a period of 8 months of treating 250 HIV-infected patients with HAART, most of whom had a CD4 count less than $200\times 10^{6}/l$, it was observed that 30% of the patients developed AIDS-defining events due to failure of HAART ^[40,45]. Thus, since from Model-Assumption 6 failure of treatment occurs for individuals who begin ART in the advanced stages of the HIV infection, in this example, we use the failure rate in ^[40]. That is, we assume that the failure rate of ART at which treated individuals develop full-blown AIDS is given by $\alpha_{TA} = 0.3$.

Note that all parameter estimates given in (a.)–(i.) above are summarized in Table 1 in years, and the trajectories of the different states $S_0(t), S_1(t), S_2(t), T(t), T(t)$ and $A(t)$ are generated by applying the Euler-Maruyama stochastic approximation scheme over 40 years from 1991 to 2031. For the parameter estimates in Table 1, the DFE $E_0 = (26.13417, 0, 0, 0, \dots, 0)$.

Note that all practical simulation interpretations will be extended by extrapolation only over 40 years from 1991 until 2031, with initial conditions in Eq (7.1).

7.1.   Prediction and sensitivity of the BRN $\mathfrak{R}_{0}$

Figure 6 shows the predicted BRN $\mathfrak{R}_{0}$ in Eq (6.4) over the years 1990–2022 obtained by applying the parameter estimates in Table 1 and the estimated disease transmission rate $\hat{\beta}_{0}(t)$ in Eq (7.2). Note that the prediction values for $R_{0}$ over 2019–2022 are extrapolations, while all other predictions are interpolations. Assuming that the trend of the HIV incidence rate over time does not change from Figure 3, then extrapolation poses no major cause for concern, since the fit for Eq (7.2) has a high predictive capacity for new HIV incidence rates ($R^{2}_{Prediction} = 0.9937$). Clearly, Figure 6 depicts the BRN $\mathfrak{R}_{0}$ continuously decreasing over time, assuming that all the other parameters in Table 1 remain fixed. Therefore, there is more control over the HIV/AIDS epidemic in the population over time in this scenario.

Figure 7 depicts the sensitivity of the BRN in Eqs (6.4) and (6.20) to the continuous changes in (1) the delay after infection, until the onset of ART $\tau_{2}$, (2) the proportion getting treatment $1-\varepsilon_{0}$, and (3) the intensity of the noise in the supply of ODAs $\sigma_{\varepsilon_{j}}, j\in I(0, 2)$. Note that for simplicity, the estimated value for $\hat{\beta}_{0} = 0.0211$ obtained from (7.2) in 1990 is used for this purpose. It is easy to see that $\beta_1 = 0.05 \beta_0 = 0.00106$ and $\beta_2 = 0.4\beta_0 = 0.00844$.

As remarked in Remark 2, observe in Figure 7(i) that early treatment (small $\tau_{2}$ values) leads to $R_{0} < < 1$. In Figure 7(ii), when ODAs are readily available and more people get ART (i.e. $1-\varepsilon_{0}\rightarrow 1$), then $R_{0} < < 1$. In the Figure 7(iii), when ODA supply varies less (i.e. the intensities $\sigma < < 1$), then $R_{0} < < 1$.

These observations confirm the conclusions in Remark 2 that the disease is more easily eradicated, whenever (1) infected individuals begin early ART treatment, (2) more people get ART, and (3) the supply of ODAs does not fluctuate too much over time.

7.2.   Numerical results for the stochastic stability of the DFE and disease eradication

Applying the parameters in Table 1 to (14.6)-(14.7), and selecting $0 < K < 2.6836$, leads to $\lambda(\mu) = 0.25286$. Furthermore, the modified BRN $\mathfrak{R}_{0}$ and other disease control parameters in (6.20)-(6.25) are respectively, $R_{0, \varepsilon} = 0.0239+2.18285\times 10^{-15} < 1$, $U_{0} = 0.0857 < 1$, $U_{1} = 0.0010629 < 1$, $U_{2} = 0.0010637 < 1$, $V_{0} = 0.08175 < 1$ and $W_{0} = 0.36623 < 1$.

Thus, the conditions of Theorem 6.3 and Theorem 6.4 are satisfied and as a result, the DFE $E_0$ is stochastically stable in probability in the large, and the disease dies out asymptotically. Moreover, from Theorem 6.4, $S(t) = S_{0}(t)+S_{1}(t)+S_{2}(t)$ is persistent and stable in the mean. These results are exhibited in Figure 8. Clearly in Figure 8(e-1, f-1, g-1), the paths of $I, T, A$ die out asymptotically. In Figure 8(d-1), the path of the total susceptible state $S(t) = S_{0}(t)+S_{1}(t)+S_{2}(t)$ is persistent in the mean, and approaches the DFE $S^{*}_{0} = 26.13417$.

7.3.   Effects of the information interventions with and without ART treatment on the HIV/AIDS epidemic

The comparative effects of information intervention and ART treatment to control the HIV/AIDS epidemic in the population are exhibited in Figure 9, where the predicted BRNs $R_{0}$ and $R^{noTr}_{0}$ in Eqs (6.4) and (6.5) for the HIV/AIDS dynamics without ART treatment in Eq (3.13) and HIV/AIDS dynamics with ART treatment in Eqs (4.3)–(4.9). Clearly, the predicted BRN $\mathfrak{R}_{0}$ when ART treatment is available is continuously lower than the BRN $\mathfrak{R}_{0}$ when there is no ART treatment available as remarked in Remark 1. Thus, the disease is more easily eradicated, whenever IECs and also ART treatment are available in the population.

7.4.   A statistical analysis of the power of the model to predict the BRN $\mathfrak{R}_{0}$

To give more credibility to the model to predict the BRN $\mathfrak{R}_{0}$ over time using the fit in Eq (7.2), we apply the statistical technique of data imputation to assess the predictive power of Eq (7.2) for the BRN $\mathfrak{R}_{0}$ over time. The method consists of the following.

(1) The data for HIV incidence per 1000 people depicted in Figure 3 is subdivided into two parts namely: Part A (Figure 10) called the "training" data, consists of measurements for the HIV incidence per 1000 over 1990–2013. And Part B (Figure 10) called the "validation" or "testing" data, consists of measurements for the HIV incidence per 1000 over 2014–2018. The data in Part A will be treated as an incomplete data set with missing values for the HIV incidence per 1000 over 2014–2018.

(2) The data in Part A is used to fit the model $\beta_{0}(t) = b_{0}+b_{1}t+b_{2}t^{2}+b_{3}t^{3}+b_{4}t^{4}+Error$, similarly to Eq (7.2). The new fit is given by

(3) The new fit in Eq (7.4) is used to predict values for the HIV incidence/1000 people by extrapolation over the years 2014–2018, and the predicted values are used to impute the missing data in Part A over 2014–2018. Figure 10 (bottom figure) depicts the actual observed data that is also given in Figure 3, and the imputed/complete Part A data over over 1990–2018. Indeed, in the bottom figure in Figure 10, the purple circular dots are the actual observed data for the HIV incidence/1000 people shown in Figure 6. The cross black dots represent the imputed or estimated values obtained from (7.4) over the years 2014–2018. Observe that the actual observed and predicted values are close to each other, for the given sample data used. However, point estimation is limited; we use confidence intervals to affirm that predicted and actual observed values are close.

(4) The imputed/complete Part A data set (now complete over 1990–2018) is used to generate predicted 95% confidence intervals (CI) for new responses $\beta_{0}(t)$ over the years 2014–2018. The 95% CI's for new responses $\beta_{0}(t)$ over the years 1990–2018 are depicted in Figure 11. The outer thin blue dotted lines are the lower and upper 95% CI limits for the new responses $\beta_{0}(t)$ over the years 1990–2018. The inner thick red dashed lines are the lower and upper 95% CI limits for the conditional mean response $\mathbb{E}[\beta_{0}(t)|t]$ over the years 1990–2018. The thick solid black line is the fitted values of the regression model, computed on the imputed/complete Part A data over 1990–2018 in Figure 10. The purple dots are the actual observed values for the HIV incidence/1000 people, which are also depicted in Figure 3.

(5) The predicted 95% CI's for $\beta_{0}(t)$ over the years 2014–2018 in Figure 11 are compared with the actual observed values of Part B in Figure 10. Observe that all actual observed values of Part B fall within the 95% CI limits for the new responses $\beta_{0}(t)$ over the years 2014–2018. This implies that at the 5% significance level, there is sufficient evidence that predictions from the fitted model in Eq (7.4) lead to values that are close to the actual observed values of $\beta_{0}(t)$. Hence, adding the fitted model (7.2) into Eq (6.4), and predicting new values for the BRN $\mathfrak{R}_{0}$ over time, such as in Figure 6, leads to trustworthy results, provided that all other parameters in Eq (6.4) are correctly estimated.

However, observe from Figure 11 that as the years increase, the observed values tend to move further away from the 95% CI limits for the new responses $\beta_{0}(t)$ over the years 2014–2018. This suggests that extrapolation beyond four years with the model (7.2) may lead to errors in the predicted BRN $\mathfrak{R}_{0}$.

8.   Conclusions

In this study, a multipopulation behaviorial HIV/AIDS epidemic model is derived and studied. Human behavior is influenced by IECs in the population. White noise perturbations of the random supply of ODAs and also in the living standards of people in the population are considered. Moreover, the delay time to begin ART treatment after HIV infection is also studied. Stochastic stability in probability of the DFE is investigated, and conditions for disease eradication are found. The model is applied to Uganda HIV/AIDS epidemic, and simulation results are given. {The results of this study show that HIV/AIDS is more easily eradicated, whenever the two HIV/AIDS control strategies: IECs and ART treatment are combined together in the population. Moreover, when the supply of ODAs is optimal and less random, more people are encouraged to test for HIV early, and to afford ART treatment, which also controls the disease more easily.

Conflict of interests

The author declares that there is no conflict of interest in this paper.

Appendix.   9.Derivation and interpretation of the deterministic distributed delay epidemic models (3.1)–(3.7)

Indeed, the delay terms in Eqs (3.1)–(3.7) are derived as follows. The term $\beta_{j} S_{j} G_{j}(i(t)), j\in I(0, n)$ is the newly infected individuals from the susceptible state $S_{j}$ at any time $t\geq t_{0}$. The fraction $\varepsilon_{j}, j\in I(0, n)$ of the newly infected population fail to get ART and progress into full-blown AIDS after $\tau_{1}$ time units. Also, the other proportion $\bar{\varepsilon}_{j} = 1-\varepsilon_{j}, j\in I(0, n)$ of the newly infected population will get ART after $\tau_{2}$ time units, and enter into the treatment state $T(t)$.

Let $h = \max{(\tau_{1}, \tau_{2})}$. Thus, assuming exponential survival lifetime during HIV infectiousness, the total number of infectious persons at any time $t\geq t_{0}$ is

where $I_{0}$ is the initial infectious population, and $p_{\varepsilon_{j}}(t), j\in I(0, n)$ is the probability that a newly infected person who will not receive ART, remains infectious for $t$ time units. Also, $p_{\bar{\varepsilon_{j}}}(t), j\in I(0, n)$ is the probability that a newly infected person who will receive ART, remains infectious for $t$ time units. Clearly,

and

and

It is easy to see from Eqs (9.1)–(9.4) that for any time $t\geq h$, where all initial perturbations have already been converted, then $I(t)$ becomes

Recall Model-Assumption 2.5, the delays $\tau_{1}$ and $\tau_{2}$ are random variables with probability densities $f_{\tau_1}, t_0 \leq \tau_1 \leq h_1$ and $f_{\tau_2}, t_0 \leq \tau_2 \leq h_2$. Therefore, taking the average of Eq (9.5) with respect to the delays $\tau_{1}$ and $\tau_{2}$, and differentiating the result leads to Eq (3.3).

Hence, in the models (3.1)–(3.7), the delay term $\sum_{j = 0}^{n} \int_{t_{0}}^{h_{1}} \varepsilon_{j} \beta_{j} S_{j}(t-r) G_{j}(i(t-r))e^{-\mu_{I} r} f_{\tau_{1}}(r) d r$ represents the average number per unit time of individuals newly infected at earlier times $t-r, r\in [t_{0}, h_{1}]$, who survived natural death rate over the incubation period, $\tau_{1}$, of HIV, with exponential survival distribution $e^{-\mu_{I} r}, \forall r\in [t_{0}, h_{1}]$, and are currently converted into full-blown AIDS class.

The other delay term $\sum_{j = 0}^{n} \int_{t_{0}}^{h_{2}}\left(1-\varepsilon_{j}\right) \beta_{j} S_{j}(t-u) G_{j}(i(t-u))e^{-\mu_{I}u}f_{\tau_{2}}(u)du$ represents the average number per unit time of individuals newly infected at earlier times $t-u, u\in [t_{0}, h_{2}]$, who survived natural death rate over the time lapse between initial infection and the on set of treatment, with exponential survival distribution $e^{-\mu_{I} u}, \forall u\in [t_{0}, h_{1}]$, and are now newly converted into the treatment class. All other parameters of the model are non-negative and defined in Assumptions (A)–(G) in Model-Assumptions 2.1–2.8.

10.   The white noise processes in the HIV/AIDS epidemic dynamics

From Assumptions (D) and (G) in Model-Assumptions 2.1–2.8, it is assumed there are noises in the HIV/AIDS dynamics from the random supply of ODAs and random poverty levels/living standards in the community reflected by life expectancy. That is, the proportions per unit time $\bar{\varepsilon}_{j} = 1-\varepsilon_{j}$ and $\varepsilon_{j}, j \in I(0, n)$ are random variables, and likewise the natural death rates $\mu_{k}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$ are random variables per unit time. Denote these random variables, respectively, by $\tilde{\bar{\varepsilon}}_{j}$, $\tilde{\varepsilon_{j}}, j \in I(0, n)$ and $\tilde{\mu_{k}}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$. We represent the environmental variabilities as independent white noise processes applying similar techniques in the earlier studies ^[27,28].

For $t\geq t_{0}$, let $(\Omega, \mathfrak{F}, P)$ be a complete probability space, and $\mathfrak{F}_{t}$ be a filtration (that is, sub $\sigma$- algebra $\mathfrak{F}_{t}$ that satisfies the following: given $t_{1}\leq t_{2} \Rightarrow \mathfrak{F}_{t_{1}}\subset \mathfrak{F}_{t_{2}}; E\in \mathfrak{F}_{t}$ and $P(E) = 0 \Rightarrow E\in \mathfrak{F}_{0}$). The variabilities in $\tilde{\bar{\varepsilon}}_{j}$, $\tilde{\varepsilon_{j}}, j \in I(0, n)$ and $\tilde{\mu_{k}}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$ in any small time interval of length $dt$ are represented by independent white noise processes as follows:

and the $w_k(t)$'s are the normalized Wiener processes for the $k^{th}$ state at time $t$ ($k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$), with the following properties: $w_{k}(0) = 0, E(w_{k}(t)) = 0, Var(w_{k}(t)) = t$.

Note from Eqs (10.1) and (10.2) that the random variables $\tilde{\bar{\varepsilon}}_{j}\in \mathbb{R}$, $\tilde{\varepsilon}_{j}\in \mathbb{R}, j \in I(0, n)$ and their means $\mathbb{E}[\tilde{\bar{\varepsilon}}_{j}] = \bar{\varepsilon}_{j} = 1-\varepsilon_{j}\in (0, 1)$ and $\mathbb{E}[\tilde{\varepsilon}_{j}] = \varepsilon_{j}\in (0, 1)$.

Furthermore, to emphasize that the intensities $\sigma_{\varepsilon_j}$ and $\sigma_{\bar{\varepsilon}_j}$ are identical, it is easy to see from Eqs (10.1) and (10.2) that the variances of the random fluctuations in the random variables $\tilde{\bar{\varepsilon}}_{j}$, $\tilde{\varepsilon_{j}}, j \in I(0, n)$ and $\tilde{\mu_{k}}, k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$ in any small time interval of length $dt$, are given by $\mathbb{V}ar[\tilde{\mu_{k}}dt] = \sigma^{2}_{k}dt$, where $\sigma_k^{2}$ is the intensity of the environmental white noise in the natural death rate of the $k^{th}$ state; $\mathbb{V}ar[(1-\tilde{\varepsilon_{j}})dt] = \mathbb{V}ar[\tilde{\varepsilon_{j}}dt] = \sigma^2_{\varepsilon_j}dt\equiv \sigma^2_{\bar{\varepsilon}_j}dt$, where $\sigma^2_{\varepsilon_j}\equiv \sigma^2_{\bar{\varepsilon}_j}$ is the intensity of the white noise in the random variable $\tilde{\varepsilon_j}, j\in I(0, n)$. Note, $\sigma_{\varepsilon_j}$ and $\sigma_{\bar{\varepsilon}_j}$ are identical, however, the distinct notations are utilized to emphasize the origins of the noises.}

Note, as remarked in ^[28], the traditional use of the white noise process in the form as additive noise, given above in Eqs (10.1) and (10.2), to represent environmental variabilities in the biological system, is only for convenience. Some authors ^[46,47] have raised important concerns about the suitability of this form to represent biological processes. More advanced extensions to this technique of modeling environmental variabilities with the white noise process have been proposed, e.g. the mean-reverting approach ^[46,47].

Remark 10.1. Note that from Model-Assumption 2.7, it makes sense that the expected value $\mathbb{E}\left(\mu_{k}dt + \sigma_{k}dw_{k}(t)\right) = \mu dt, \forall k \in\left\{S_{j}, I, T, A, R\right\}, j\in I(0, n)$, where $\mu$ is the mean of the random variable natural death rates per unit time $\tilde{\mu_{k}}$ of every $k^{th}$ state, $k \in\left\{S_{j}, I, T, A, R\right\}, j\in I(0, n)$. That is, Eq (10.1) can be rewritten as

To emphasize the origins of the deathrates, the expressions in Eq (10.1) will be used, and wherever necessary, remarks for Eq (10.3) will be given.

To avoid repetition, the reader is referred to Wanduku ^[28], where it is shown that for any state $k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$, even when the natural death rates are influenced by random fluctuations of white noise types exhibited in Eqs (10.1) and (10.3), it follows that under the assumption of independent natural deaths in any interval $[t_{0}, t_{0}+t]$ of length $t$, the number of deaths that occur in state $k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$ are best described by independent Poisson processes with means $E(\tilde{\mu}_{k}) = \mu_{k}$, where $\tilde{\mu}_{k}$ and $\mu_{k}$ are defined in Eq (4.1). Thus, the survival lifetimes until natural death occurs in the $k^{th}$ state at time $t$ ($k \in\left\{S_{j}, I, T, A, R\right\}, j \in I(0, n)$) are given by

11.   Complete proof and interpretation of Theorem 5.1

Proof.

Observe from Eqs (4.3)–(4.9) and (4.12) that for the intensities $\sigma_{\varepsilon} = \sigma_{\bar{\varepsilon}}\geq 0$, $\sigma_{S_{j}}\geq 0$ and $\sigma_{k}\geq 0, k\in \{I, T, A, R\}$, then

When the solution of Eqs (4.3)-(4.9) is positive, i.e. $Y(t)\in \mathbb{R}^{n+6}_{+}$, then it is easy to see from Eq (4.12) that when Theorem 5.1 (a.) and (b.) holds (i.e. $\sigma_{k} = 0$, $\forall k\in \{S_{j}, I, T, A, R\}$, and $\forall j\in I(0, n)$, and also $\sigma_{\varepsilon_{j}}\equiv\sigma_{\bar{\varepsilon}_{j}}\geq 0$), then given that $N(t_{0})\in [0, \frac{B}{\mu_{min}}$], it follows that

where

Furthermore, it is easy to see from Eqs (11.2), (3.10) and (4.9) that for $Z(t_{0})\in [0, \frac{1}{\mu_{Z}}\frac{\phi_{max}}{\min{(1, \hat{\phi}_{min})}}\frac{B}{\mu_{min}}]$,

where $\phi_{max} = \max(\phi_{I}, \phi_{T}, \phi_{A})$, and $\hat{\phi}_{min} = \min(\hat{\phi}_{I}, \hat{\phi}_{T}, \hat{\phi}_{A})$.

Thus, it follows from Eqs (11.2) and (11.4) that the closed ball centered at the origin with radius

defined as follows

is positive self-invariant with respect to the stochastic systems (4.3)–(4.9), whenever the conditions of Theorem 5.1 (a.) and (b.) holds holds.

On the other hand, when Theorem 5.1 (c.) holds (i.e. at least one of $\sigma_{k} > 0$, $k\in \{S_{j}, I, T, A, R\}$, and $j\in I(0, n)$, and also $\sigma_{\bar{\varepsilon}_{j}}\equiv\sigma_{\varepsilon_{j}}\geq 0$), then the closed ball $\mathbb{\bar{B}}_{\mathbb{R}^{n+6}}(0, r)$ is no longer positively self-invariant Eqs (4.3)–(4.9). However, the positive paths of the SSP $Y(t)\in \mathbb{R}^{n+6}_{+}, t\geq t_{0}$ continue to oscillate in $\mathbb{R}^{n+6}_{+}$ for all time $t\geq t_{0}$.

Remark 11.1.

(1.) Theorem 5.1 signifies that the stochastic systems (4.3)–(4.9) always has a positive solution $\{Y(t, w), t\geq t_{0}, w\in \Omega\}$. Moreover, Theorem 5.1(a)-(b) signify that all paths of $\{Y(t, w), t\geq t_{0}, w\in \Omega\}$ that start in the closed ball $\mathbb{\bar{B}}_{\mathbb{R}^{n+6}}(0, r)$, continue to oscillate in the closed ball, almost surely, provided that the intensities of the noises in the system are infinitesimally small (i.e. $\sigma_{\varepsilon_{j}} = 0$, $\sigma_{S_{j}} = 0, \forall j\in I(0, n)$ and $\sigma_{k} = 0, \forall k\in \{I, T, A, R\}$), or provided that the only source of variability in the system is from the random supplies of ODAs in the system with intensities $\sigma_{\varepsilon_{j}} > 0$ (i.e. $\sigma_{\varepsilon_{j}} > 0$, $\sigma_{S_{j}} = 0, \forall j\in I(0, n)$ and $\sigma_{k} = 0, \forall k\in \{I, T, A, R\}$).

(2.) When additional sources of noise in the system occur e.g. the noises in the living standards of people in system, i.e. at least one of $\sigma_{k} > 0$, $\forall k\in \{S_{j}, I, T, A, R\}, j\in I(0, n)$, whenever $\sigma_{\bar{\varepsilon}_{j}}\equiv\sigma_{\varepsilon_{j}}\geq 0, , j\in I(0, n)$, then the paths of $\{Y(t, w), t\geq t_{0}, w\in \Omega\}$ are inflated out of bounds by the additional noises, to continue oscillating outside of the closed ball $\mathbb{\bar{B}}_{\mathbb{R}^{n+6}}(0, r)$, but remain in the positive phase space $\mathbb{R}^{n+6}_{+}$.

(3.) The observations in (1)-(2) above suggest that the stochastic HIV/AIDS epidemic model in Eqs (4.3)–(4.9) is more profound than the corresponding deterministic systems (3.1)–(3.7). For instance, (i.) randomness exclusively in the supply of ODAs does not perturb the system much different from the corresponding deterministic scenario (this fact is exhibited in the example in Section 7). In fact, in this case, the disease outbreak cannot grow beyond bounds defined by $\mathbb{\bar{B}}_{\mathbb{R}^{n+6}}(0, r)$ in Eq (5.1). (ii.) The occurrence of additional variability e.g. from the living conditions of people in the system, tends to inflate paths of the disease dynamics out of bounds into the unbounded space $\mathbb{R}^{n+6}_{+}$. This case suggests the possibility of extinction of the population occurring, whenever for instance, the living standards fluctuate with very high intensities in the HIV/AIDS epidemic.

12.   Complete proof and interpretation of Theorem 6.1

Proof. When the stochastic systems (4.3)–(4.9) is in equilibrium, then (1) the drift and diffusion components corresponding to the same equation in the system (4.3)-(4.9) must be equal. Moreover, the disease-free steady state requires that (2) the corresponding disease related states in the equilibrium point $E_{0} = X_{0}^{*} = (S^{*}_{0}, \ldots, S^{*}_{n}, I^{*}, T^{*}, Z^{*})$ satisfy $I^{*} = T^{*} = Z^{*} = 0$. It is easy to see that solving the resulting system of equations obtained by applying (1)-(2) above leads to $E_{0}$ in Eq (6.1), whenever the conditions in Theorem 6.1(a.)–(c.) hold, and there is no solution, whenever Theorem 6.1(d) holds ^[28].

Remark 12.1. In essence, Theorem 6.1 signifies that the stochastic and the corresponding deterministic systems (4.3)–(4.9) and (3.1)–(3.7), respectively, have the same DFE $E_{0}$, whenever $\sigma_{S_{0}} = 0$. Furthermore, when $\sigma_{S_{0}} > 0$, the DFE no longer exists. This result suggests that the magnitude of the intensity ($\sigma_{S_{0}}$) of the fluctuations in living standards of the susceptible people of type $S_{0}$, who have not yet modified their sexual behaviors, determines the existence of a disease-free steady state in the population, wherein the disease can be eradicated. If during the HIV/AIDS epidemic outbreak, the majority of the susceptible population ($S_{0}$) live in unhealthy conditions, so that more people tend to die naturally, then the intensity, $\sigma_{S_{0}} > 0$, of the noise in $S_{0}$ becomes significantly high, and the DFE ceases to exist. This implies that HIV/AIDS will persist in the population.

13.   Complete proof of Lemma 6.1

Proof. It is easy to see that the differential operator applied to $V_{1}$ with respect to Eqs (4.3)-(4.9), leads to

where

Similarly,

where

Also,

where

And

In the following, Cauchy-Schwarz inequality, Holder inequality, $(a+b)^{2}\leq 2a^{2} +2b^{2}$, and the following algebraic inequality Eq (13.8) will be used to obtained the results.

Indeed, it follows from Eqs (13.2) and (13.4) that

Also, it follows from Eq (13.6) that

where $\lambda(\mu)$ and $\lambda_{2}(\mu)$ are positive real valued functions of $\mu$.

Also, from Eq (13.7), it is easy to see that

From Eqs (13.9) and (6.7), observe that for each $j\in I(0, n)$, the terms $-2 \gamma_{j} H_{j}(Z(t)) \leq \gamma_{j}\left(1+\left(H_{j}^{*}\right)^{2}\right)$, $\gamma_{0} S_{0}^{*} H_{j}^{2}(Z(t))\leq \gamma_{0} S_{0}^{*}Z^{2}(t)$, and $\beta_{0} S_{0}^{*}G_{0}^{2}(i(t))\leq \beta_{0} S_{0}^{*}(i(t))^{2}\leq \beta_{0} S_{0}^{*}I^{2}(t)$. Now, combining Eqs (13.2)–(13.11) and grouping like-terms, and also applying Eq (6.7) and some further algebraic manipulations and simplifications, it is easy to see that

where $\phi_{j}, j\in I(0, n)$ and $\varphi_{k}, k\in \{I, T, Z\}$ are defined in Eqs (6.11)–(6.15).

Note 13.1.

(1.) Recall Theorems 6.1 & 5.1 assert that the DFE exits in $D(\infty)$ and $\mathbb{R}^{n+6}_{+}$ under a given set of restrictions for the intensities of the noises in the system. Furthermore, the modes of estimating the delay-term in Eq (6.18), $\sum^{n}_{j = 0}\left(2\varepsilon_{j}\beta_{j}\frac{1}{\lambda(\mu)} +\sigma^{2}_{\varepsilon_{j}}\beta^{2}_{j}\right)\mathbb{E}_{\tau_{1}}\left[S^{2}_{j}(t-\tau_{1})G^{2}_{j}(i(t-\tau_{1}))e^{-2\mu_{I}\tau_{1}}\right]$ give rise to different pictures of the disease dynamics inside the different spaces $D(\infty)$ and $\mathbb{R}^{n+6}_{+}$.

Indeed, in $D(\infty)$, using Eqs (5.2) and (6.7)

where $r$ is defined in Eq (5.2).

Due to space limitations, the second estimation of (13.13) and corresponding interpretation of the disease dynamics in $\mathbb{R}^{n+6}_{+}$, whenever Theorem 6.1(c.)-(d.) hold appears in a sequel to this paper. Only the stochastic stability in probability of the DFE $E_{0}$ in the space $D(\infty)$ will be discussed in this paper. That is, stability in probability, whenever Theorem 6.1(a.)-(b.) hold.

(2.) In subsequent analysis, it is assumed that Theorem 5.1(a.)-(b.), Theorem 6.1(a.)-(b.) and Theorem 6.2 hold. That is, the stochastic systems (4.3)-(4.9) has a positive solution and a DFE $E_{0}$ inside of $D(\infty)$.

14.   The complete proof of Lemma 6.2

Proof. Let $R_{0, \sigma_{\varepsilon}} < 1$, $R_{1}(\mu) < 1+(1-R_{0, \sigma_{\varepsilon}})$, $\max{(U_{j})}_{j\in I(0, n)} < 1$, and $\max{(V_{0}, W_{0})} < 1$. When the intensities of Eqs (4.3)-(4.9) satisfy $\sigma_{\bar{\varepsilon}_{j}}\equiv\sigma_{\varepsilon_{j}}> 0$, $\sigma_{S_{j}} = 0, \forall j\in I(0, n)$ and $\sigma_{k} = 0, \forall k\in \{I, T, A, R\}$, observe from Eqs (6.16)-(6.19) and (13.13) that

where $\phi_{j}, j\in I(0, 1)$ and $\varphi_{k}, k\in \{I, T, Z\}$ are obtained from Eqs (6.11)-(6.15) by setting the intensities $\sigma_{\bar{\varepsilon}_{j}}\equiv\sigma_{\varepsilon_{j}}> 0$, $\sigma_{S_{j}} = 0, \forall j\in I(0, n)$ and $\sigma_{k} = 0, \forall k\in \{I, T, A, R\}$. Moreover, it is easy to see from Eqs (6.11)-(6.15) and (14.1) that

and

Select the functions $\lambda(\mu)$ and $\lambda_{2}(\mu)$ such that in Eqs (14.4), (6.4) and (6.3), $P(1)$ is proportional to $\frac{\left[(1+2\frac{1}{\lambda(\mu)}+\lambda(\mu))\right]}{(1-\frac{1}{2}\lambda_{2}(\mu))}$. That is,

where

and $\lambda(\mu)$ satisfies

Note that when $\lambda_{2}(\mu)\in (0, 2)$ is given, then the solutions of Eq (14.7) are real and positive. Also, clearly from Eqs (14.4) and (6.4)-(6.3),

It is easy to see from Eq (14.4) that

15.   Complete proof of Theorem 6.4

Proof. It is easy to see from Eq (6.26) that

where

Integrating both sides of Eq (15.1) over $[t_{0}, t]$, and diving both sides by $t$, it is easy to see that

Observe that for all $Y(t)\in D(\infty)$, and from Eq (6.17), the quadratic variation of the local martingale $m(t) = \int^{t}_{t_{0}}\overrightarrow{g}(X(\zeta))d\overrightarrow{w(\zeta)}$ denoted by $< m(t)>$ is given by

For for all $Y(t)\in D(\infty)$, it is easy to see that $\limsup_{t\rightarrow \infty}\frac{1}{t} < m(t)> < \infty$, a.s., and by the strong law of large numbers for local martingales ^[48], $\lim_{t\rightarrow \infty}\frac{1}{t} m(t) = 0$, a.s.

Thus, dividing both sides of Eq (15.3) by the constant $m_{0}\left(\frac{1}{(n+1)}\right)>0$, and taking $\limsup_{t\rightarrow \infty}$ of both sides, we obtain

Clearly, from Eq (15.5), $\limsup_{t\rightarrow \infty}{\frac{1}{t}\int^{t}_{t_{0}}{\left(S(\zeta)-S^{*}_{0}\right)^{2}d\zeta }} = 0$, a.s. Furthermore, since the integrand is nonnegative, i.e. $\left(S(t)-S^{*}_{0}\right)^{2}\geq 0, \forall t\geq t_{0}$, and the random variable $S(t, w), w\in \Omega, t\geq t_{0}$ is continuous with respect to $t\geq t_{0}$, but nowhere smooth, it implies that $S(t)$ must be equal to $S^{*}_{0}$ almost everywhere with respect to $t\geq t_{0}$, and (6.27) is obtained.