Persistence and extinction of infection in stochastic population model with horizontal and imperfect vertical disease transmissions

1.   Introduction

Epidemic models deal with the transmission of infection in a population. They are frequently used to gain insights into disease dynamics and control mechanisms. Susceptible-infected (SI) type deterministic compartmental epidemic models with constant rate parameters have been extensively used to explore infection-related issues in different populations ^[1,2,3]. However, disease transmission is a random process that varies randomly due to exogenous factors. For example, disease transmission might be significantly affected by environmental factors like temperature, humidity, rainfall, etc. ^[4,5]. We cannot fully explain all the unknown factors in disease transmission through deterministic rules alone. Therefore, it is natural to use epidemic models incorporating random variations to account for environmental unpredictability. Incorporating such random variation in mathematical models leads to stochastic differential equations. Stochastic models can be derived by perturbing the model parameters randomly. An alternative approach consists of perturbing the state variables' equations, but the interpretability of the effect of the noise is more complicated here ^[6,7,8].

Disease transmission from one infected individual to another susceptible individual of a population is a significant issue, and the persistence of parasites and their virulence largely depend on this transmission mechanism ^[9,10,11]. Horizontal and vertical transmissions are two utterly distinct disease transmission modes usually followed by different pathogens ^[12]. Here, we consider a homogeneous mixture of susceptible and infected populations, where a disease spreads following both horizontal and vertical transmissions. In horizontal transmission, the infection spreads through contact. On the contrary, vertical transmission occurs through birth only. Though many parasites spread disease through multiple pathways, horizontal transmission is the predominant infection-spreading mode ^[13]. Infected hosts can give birth to susceptible and infected individuals, so the vertical transmission may be imperfect ^[9]. Stochastic population models with only horizontal disease transmission have been proposed in ^[14,15]. These researchers proved the existence of a stationary distribution condition for eradicating the disease from the population. A simple stochastic susceptible-infected-recovered (SIR) model with horizontal and vertical transmissions was considered in ^[16]. This model, however, entailed the linear birth rate of all populations, contrary to the realistic density-dependent birth. It is shown that considerable noise is conducive to controlling the infection. Li et al. ^[17] analyzed an SIRS epidemic model with a nonlinear incidence rate and environmental stochasticity and only horizontal transmission of infection with density-independent growth of susceptible populations. A set of sufficient conditions for the disease extinction and persistence is provided, and simulation results are presented.

Our main objective is to study a susceptible-infected-susceptible type compartmental population model with both disease transmission modes and uncertainty/stochasticity in the horizontal disease transmission term. The deterministic version of this model has similarities to the model studied by ^[2,18,19]. These researchers obtained the local stability results of the deterministic models and studied the trade-offs between modes of transmission and virulence. We determine the local and global stabilities of the corresponding deterministic system to compare its results with the stochastic model. A set of conditions is provided under which the deterministic solution becomes a limiting case of the stochastic solution. We also study the law of the disease extinction time. It is experimentally demonstrated that the relative fecundity of parasites plays a significant role in the persistence of infection and survival of both the host and parasite populations ^[19,20,21]. One way to measure the parasites' fitness is to determine the extinction time of the infected host population. We show how the relative fecundity and the disease transmissibility jointly affect the parasite fitness.

The rest of the paper is organized as follows. The SIS epidemic model is presented in Section 2. A brief study of the deterministic model is presented in Section 3. The stochastic system analysis is given in Section 4. The analytical results are illustrated in Section 5 including local and global sensitivity analyses. A case study for the imperfect vertical and horizontal infection transmission is provided in Section 6 to fortify our observed theoretical findings, and the paper ends with a discussion in Section 7.

2.   A Stochastic SIS model with horizontal and vertical transmission

The host population is divided into susceptible and infected classes in the presence of some parasite infection. Let $S(t)$ and $I(t)$ denote the densities of the susceptible and infected host populations at any time $t$, respectively. Infected hosts give birth to susceptible and infected individuals at rates $b_I$ and $e$, respectively. Parameter $e$ thus represents the vertical transmission. The birth rate of the susceptible host is denoted by $b_S$. Parasites may affect the fecundity and morbidity rates of the host population ^[22,23]. It is therefore assumed that the death rate of infected hosts is never less than that of susceptible hosts, i.e., $u_I\geq u_S$, and the birth rate of susceptible hosts ($b_S$) is higher than the total birth rate of infected hosts, i.e., $b_S\geq b_I+e$. An infected individual may recover from the infection and join the susceptible class to be reinfected. Let $\beta$ be the horizontal disease transmission coefficient, $K$ be the environment's carrying capacity, and $\mu$ be the recovery rate of infected individuals. The following coupled nonlinear differential equations can represent the SIS epidemic model:

This deterministic model has similarities to the model studied by Lipsitch et al. ^[2], where they considered $\mu = 0$ so that the underlying epidemic system is of SI type. They obtained the local stability results of the SI epidemic model and observed the implications of different transmission modes for the evolution of virulence. Mangin et al. ^[18], and Ebert et al. ^[19] considered a similar deterministic SI epidemic model with horizontal transmission only (where $e = 0$ and $\mu = 0$) to understand the trade-offs between modes of transmission and virulence. A discrete version of the model (2.1) with $\mu = 0$ was also analyzed in ^[24]. Thus, the SIS model with both transmission modes is a modification of earlier models. However, our main interest is in the stochastic version of the above model. In the following, we proved the global stability of the equilibrium points of the system (2.1), which has not been studied.

Horizontal disease transmission may occur directly through the effective contact between a susceptible and an infective or indirectly through the environment and intermediate host populations ^[25]. Such transmissions are primarily unknown and random. Assuming that the disease transmission process is random and directly affects the horizontal disease transmission parameter $\beta$, we replace the parameter $\beta$ by $\beta+\sigma d\xi_t$, where $\sigma$ is a real constant that measures the intensity of the noise and $\xi(t)_{t\geq 0}$ is a standard Wiener process defined on a complete probability space $(\Omega, \mathit{ F, P)}$ with a filtration $\mathit{ \{F_t\}}_{t \in \Re^+}$ and satisfies $< d \xi(t) > = 0$, $< d\xi(t), d\xi(t') > = \delta (t-t')$, where $\delta$ is the Dirac delta function. It induces a multiplicative noise and thus avoids the negativity of the solution due to initial negative fluctuations. Zero becomes a lower bound in this case, even if the initial population size is small and there is a negative fluctuation of noise ^[26]. Under this assumption, the stochastic SIS epidemic model (2.1) with horizontal transmission and imperfect vertical transmission becomes

Here, we analyze stochastic system (2.2), which considers density-dependent growth of host populations and horizontal as well as imperfect vertical transmissions of infection. We provide the local and global dynamics of the corresponding deterministic system with respect to the deterministic basic reproduction number. The persistence and extinction conditions of disease for the stochastic system are derived. The stochastic process' ergodicity and the stationary distribution criteria are proved. We present different simulation results to illustrate the theoretical results and provide insights into parasites' fitness.

3.   Study of the deterministic model

Deterministic system (2.1) includes an additional recovery term of the infected host with respect to the model studied in ^[2]. Local stability of different equilibrium points of the system (2.1) with no recovery (i.e., $\mu = 0$) has been proved ^[2,24]. Here, we start by recalling (without proof) the local dynamics of system (2.1), where recovery is considered, and then provide a new global stability result.

3.1.   Basic reproduction number

The basic reproductive number is an essential measure in epidemics and is critical in eradicating infection. It is defined as the number of new cases arising from a single infected individual when introduced into a susceptible population group ^[27]. The disease cannot be established in the host population if the basic reproduction number is less than unity. Otherwise, the disease can invade the population. We measure the basic reproduction number of deterministic system (2.1) through the next-generation matrix ^[28].

Proposition 1. The basic reproduction number, $R^{D}_0$, of deterministic system (2.1) is

Proof. In the absence of infection, the susceptible population has equilibrium density $\hat{S} = K\left(1-\frac{u_s}{b_s}\right)$. The Jacobian matrix of the epidemic model (2.1) evaluated at $(\hat{S}, 0)$ is

The sub-matrix of $J_{11}$ associated with the infectious compartment is a one-by-one matrix

where $F = b_I\left(1-\frac{\hat{S}}{K}\right)+\beta\hat{S}$ and $V = (u_I+\mu)$. The next-generation matrix is then given by

The basic reproduction number of system (2.1) is the spectral radius of the scalar matrix $FV^{-1}$ ^[28] and is thus

□

As mentioned earlier, infection will die out from the population if the basic reproduction number is less than $1$. A straightforward way to reduce the basic reproduction number is to decrease the horizontal transmission coefficient $\beta$. Increasing parameter $\mu$ can also reduce the value of the basic reproduction number and may be a parameter of interest.

3.2.   Local and global stability of the equilibrium points

From a biological point of view and for the model's applicability, solutions should be non-negative at any time, and the deterministic model solutions should stay bounded. For this, one can easily prove the following result in the line of ^[29].

Proposition 2. Solutions of system (2.1) are positively invariant and uniformly bounded in the domain $G \subset \mathbb{R}^{2}_{+}$, where $G = \bigg\{(S, I)\in \mathbb{R}^{2}_{+}\; \rvert\; 0\leq S(t), I(t)\leq K\bigg\}.$

The local stability results (without proof) may be summarized as follows.

Theorem 1. Model (2.1) has three non-negative equilibrium points.

(i) The trivial equilibrium point $E_0^{(1)}$ = $(0, 0)$ is locally asymptotically stable if $b_S < u_S$ and $b_I < u_I+\mu$.

(ii) The disease-free equilibrium point $E_1^{(1)} = (\hat{S}, 0)$, where $\hat{S} = K(1-\frac{u_S}{b_S}),$ $u_S < b_S$, is locally asymptotically stable if $R^{D}_0 < 1$.

(iii) The infected (interior) equilibrium point $E_*^{(1)} = (S^*, I^*)$, having equilibrium densities

exists and becomes locally asymptotically stable under conditions

where $A = \frac{\beta K}{b_I ^2}\; [\beta K(b_I+e)+b_I(b_S-b_I-e)]$, $B = -K(b_S-u_S)+\frac{K}{b_I}[(b_I-u_I-\mu)(K\beta+b_S+e)-(e+\mu)(\beta K-b_I)]+\frac{2Ke(\beta K-b_I)}{{b_I}^2}(b_I-u_I-\mu)$, $C = - \frac{K^2}{{b_I}^2}[b_I-u_I-\mu][\mu (b_I+e)+eu_I].$

Biologically, all populations go extinct if the birth rate of susceptible hosts is less than its death rate and the ensemble clearance rate of infected hosts exceeds its birth rate. Noticeably, the stability of the trivial equilibrium ceases the existence of the other two equilibrium points.

The disease-free equilibrium point exists if the birth rate of susceptible hosts is larger than its death rate and becomes stable if the basic reproduction number is less than unity. It is worth mentioning that local stability criteria of $E_1^{(1)}$ makes $E_*^{(1)}$ and $E_0^{(1)}$ unstable.

The infected equilibrium $E_*^{(1)}$ exists, and both populations persist in a stable state if the birth rate of susceptible hosts is higher than the total birth rate of infected hosts, i.e., $b_S\geq b_I+e$. In addition, the vertical birth rate is higher than its ensemble clearance rate of infected hosts but lower than the maximum disease transmission rate through horizontal transmission, i.e., $\mu+u_I < b_I < \beta K$, and the basic reproduction number is greater than unity. Thus, the basic reproduction number $R^{(D)}_0$ must be less than unity to make the system disease-free. This is achievable through parameter $\mu$, which transforms an SI system into an SIS one.

Local asymptotic stability guarantees that the solutions will reach equilibrium if they start close to the equilibrium value. However, if the considered initial values are far from equilibrium, they may not arrive at the equilibrium point. To assure that the equilibrium point's stability does not depend on the initial values, it is necessary to show its global stability. We prove that the local stability results of different equilibrium points are sufficient for their global stability.

Theorem 2. Each equilibrium of the model (2.1) is locally and globally asymptotically stable throughout domain $G \subset \mathbb{R}^{2}_{+}$, unless the solution starts from the other two equilibrium points.

Proof. We first show that there is no periodic orbit in $D \subset G$, where $D$ = $\{ (S, I)$: $0 < S < K$, $0 < I < K \}$. For this, we define a continuously differentiable function $B(S, I)$ in $D$, where $B(S, I)$ = $1/(SI)$. Defining $F = (F_1, F_2)$, where $F_1 = b_S\{(1-\frac{S+I}{K})-u_S-\beta I \}S+eI[1-\frac{S+I}{K}]+\mu I$, $F_2 = b_I I[1-\frac{S+I}{K}]-u_I I+\beta SI-\mu I$, and noting that $I\leq K\; \mbox{by Proposition}\; 2,$ one then has

Hence, by Bendixson-Dulac criteria ^[30], system (2.1) has no periodic orbit in the interior of $D$. Since $E_0^{(1)}$ is the only stable equilibrium, there is no periodic orbit in the domain of definition. Therefore, $E_0^{(1)}$ is globally asymptotically stable whenever it is locally asymptotically stable. Similar arguments prove that local stability criteria of $E_1^{(1)}$ and $E_*^{(1)}$ also assure their global stability if the solutions are not started from the equilibrium point $E_*^{(1)}$ for the first case and $E_1^{(1)}$ for the second case. □

4.   Study of the stochastic model

In this section, we present our major results of the SIS epidemic model with both modes of disease transmission. We use the Lyapunov analysis method, Ito's formula, Chebyshev's inequality, the law of large numbers, and the other standard techniques to obtain new results for system (2.2).

A population explosion may occur in the case of a multiplicative noise ^[31]. It is, therefore, necessary to show that such an explosion does not happen here. Also, biological populations should always be nonnegative. For this, we prove the stochastic solutions' global existence, positivity, and boundedness.

Proposition 3. For any initial value $(S(0), I(0))$ $\in$ $\mathbb{R}^2_+$, there exists a unique solution $(S(t), I(t))$ of the system (2.2) for $t\geq 0$, and the solution remains in $\mathbb{R}^2_+$ with probability $1$.

Proof. Since the coefficients of system (2.2) satisfy local Lipschitz conditions, there is a unique positive local solution on $[0, \nu _e)$, where $\nu_e$ is the explosion time ^[32]. We show that this solution is global, meaning $\nu_e = \infty$. One gets from (2.2)

giving

Let $m_0 > 0$ be sufficiently large so that $S(0)$, $I(0)$ lie within the interval $[\frac{1}{m_0}, m_0]$. For any integer $m \geq m_0$, define a sequence of stopping times by

where we set inf $\Phi = \infty$ ($\Phi$ represents the empty set). Since $\nu_m$ is non-decreasing as $m \to \infty$, one gets $\nu_\infty = \lim_{n \to \infty}\nu_m.$ Then, $\nu_\infty \leq \nu_e$ almost surely (a.s.). Now, we prove that $\nu_\infty = \infty$ a.s. If this statement is violated, then there exists $M > 0$ and $\delta \in (0, 1)$ such that $P\{\nu_\infty \leq M\} > \delta.$ Thus, there is an integer $m_1 \geq m_0$ such that

Define a $C^2$ function $V: \mathbb{R}^2_+ \to \mathbb{R}_+$ by

Using Ito's formula, we have

Let $U = (u_S+u_I+\mu)+(b_S+b_I+e+\beta)K+\sigma^2 K^2$. Integrating both sides of (4.2) from 0 to $\nu_m \wedge M$, one gets

Expectation in both sides yields

Set $G_m = \{\nu_m \leq M\}$ for $m\geq m_1$ and from (4.1), we have $\mathbb{P}(G_m) \geq \delta.$ For every $\tau \in G_m$, $S(\nu_m, \tau)$, $I(\nu_m, \tau)$ are equal to $m$ or $\frac{1}{m}$, implying that $V(S(\nu_m, \tau), I(\nu_m, \tau))$ is no less than min $\{m-1-\; \mbox{ln} (m), 1/m-1- \; \mbox{ln} (1/m)\}$. Therefore, we have

where $1_{G_m(\tau)}$ is the indicator function of $G_m$. Then, $m \to \infty$ leads to the contradiction

Therefore, $\nu_\infty = \infty$ a.s. This completes the proof. □

Proposition 4. Solutions of system (2.2) are stochastically bounded on any time interval for any initial value $(S(0), I(0)) \in \mathbb{R}^2_+$.

Proof. From (2.2), we have

Therefore, $\lim_{t \to \infty} (S+I) \geq$ $\frac{K}{M}(M-N) = \bar{L}$ (say).

Define the function

for $(S, I) \in \mathbb{R}^2_+$, $\theta > 1$ and $U(S, I) = S^\theta+I^\theta$. By Ito's formula, we have

where $A = \theta K^{\theta-1}(b_S+e+\mu) K+\theta K^{\theta-1}(b_I+ \beta K) K+\theta (\theta-1)\sigma^2 K^{\theta+2}$.

Using Proposition 3 and from (4.3), one gets

Making $m \to \infty$, we have

implying

Defining $|Y(t)| = \left(S^2(t)+I^2(t)\right)^{\frac{1}{2}}$ and noting that $|Y(t)|^\theta = (S^2(t)+I^2(t))^\frac{\theta}{2} \leq 2^{\frac{\theta}{2}} \; \mbox{max}\; \left\{S^\theta (t), I^\theta(t)\right\} \leq 2^{\frac{\theta}{2}}\left(S^\theta+I^\theta\right),$ we obtain $\mathbb{E}\left(|Y(t)|^\theta\right) \leq 2^{\frac{\theta}{2}}\left(e^{-t}W(S(0), I(0))+A-Ae^{-t}\right).$ Thus, $\limsup_{t \to \infty} \mathbb{E}\left(|Y(t)|^\theta\right)\leq 2^\frac{\theta}{2} A < \infty.$ Therefore, there exists a positive constant $\eta$ such that $\limsup_{t \to \infty} \mathbb{E}\left(\sqrt{Y(t)}\right) < \eta.$ For any $\nu > 0$, $\omega = \frac{\eta^2}{\nu^2} > 0$, and using the Chebyshev's inequality, $\mathbb{P}(|Y(t)| > \omega)\leq \frac{\mathbb{E}\left(\sqrt{Y(t)}\right)}{\sqrt\omega},$ one gets

Hence, the result is proven. □

4.1.   Stochastic extinction and persistence of infection

One important concern in epidemiology is the eradication of infection from the population. The average infected population in the time interval $[0, t]$ is $\frac{1}{t} \int_{0}^{t}I(\kappa)\; d\kappa$. The infected population is said to be strongly non-persistent or extinct if the supremum of its limiting average population is zero. In this case, the infection is said to be removed from the system. On the other hand, if the infimum of its limiting average population is always positive, then the $I$ population is said to be strictly persistent ^[33]. In this case, the eventual average $I$ will always stay away from zero, and the infection consistently remains in the system. In the following two theorems, we provide sufficient conditions to eliminate the disease or persist almost surely in stochastic system (2.2).

Before proving these results, we state the following well-known lemma.

Lemma 1. ^[34] Suppose $w(t)\in C(\Omega\times[0, \infty), \mathbb{R}^0_+)$, where $\mathbb{R}^0_+ = \{p\; \rvert\; p > 0, p\in\mathbb{R}\}$.

(i) If there exist two positive constants $T$ and $\kappa_0$ such that

where $\alpha_i\; (1\leq i\leq n)$ are constants, then

(ⅱ) If there exist three positive constants $T, \; \kappa, \; \kappa_0$ such that

then

Using the above definition of extinction and the Lemma 1(ⅰ), we present a vital result that guarantees disease eradication in the stochastic system.

Theorem 3. If $R^S_0 = \frac{\frac{\beta^2}{2 \sigma^2}+b_I}{u_I+\mu} < 1$, then the infected population of model (2.2) goes to extinction almost surely.

Proof. Let $(S(t), I(t))$ be a solution of the system (2.2) with initial value $(S(0), I(0))$. Applying Ito's formula in the second equation of system (2.2), we have

Integration on the both sides of (4.3) from 0 to $t$, and division by $t$ yields

where, $M(t) = \int_0^t \sigma S(\tau)d\xi(\tau)$ is the local martingale with $M(0) = 0$. Moreover, $< M, M > _t = \left(\int_0^t \sigma S(\tau) d\xi(\tau)\right)^2 = \int_0^t \sigma^2 S^2(\tau) d\tau$ $\leq \sigma^2 K^2 t$. One then has $\limsup_{t \to \infty} \frac{ < M, M > _t}{t} \leq\sigma^2 K^2 < \infty$ a.s. Therefore, by the law of a large number ^[35], $\lim_{t \to \infty} \frac{M(t)}{t} = 0$ a.s. Taking limit superior on both sides of (4.3) and using the Lemma 1, we obtain

whenever $R^S_0 = \frac{\frac{\beta^2}{2 \sigma^2}+b_I}{u_I+\mu} < 1$. Thus, we have $\lim_{t \to \infty} I(t) = 0$ for $R^S_0 < 1.$ This completes the proof. □

The condition $R^S_0 < 1$ may be considered the equivalent basic reproduction number for the stochastic system. The stochastic basic reproduction number ($R^S_0$) prescribes some restrictions on the system parameters and noise. If the condition is satisfied, the infected population will go extinct, and the system will be disease-free. One can easily observe that $R^S_0$ is a decreasing function of the recovery rate, $\mu$, and the noise intensity, $\sigma$. The recovery rate may be increased by taking suitable external measures; consequently, $R^S_0$ can be made less than unity to make the eradication process possible. Similarly, additional noise may also be helpful in disease elimination. However, the disease transmission parameter $\beta$ and the birth rate $b_I$ of the infected host are positively correlated with $R^S_0$.

The following theorem is an ergodic behavior of $I$ population of stochastic system (2.2).

Theorem 4. If $\frac{b_I-u_I-\mu}{\frac{b_I}{K}+\frac{A(b_I+e)}{b_S}-\frac{A u_I}{b_S}} > 0$, the infected population $I$ persists in mean a.s., and

Proof. Integrating both the Equations of (2.2) from 0 to $t$ and dividing by $t$, one obtains

where $< x(t) > = \frac{1}{t}\int_0^t x(s)\; ds.$ Since $S$ and $I$ are stochastically ultimately bounded, $\lim_{t \to \infty}\frac{S(t)-S(0)}{t} = 0$ and $\lim_{t \to \infty}\frac{I(t)-I(0)}{t} = 0$. Eq (4.4) then becomes

Taking limit as $t$ goes to $\infty$, we have

Using Ito's formula in the second equation of system (2.2), one gets

where $A = \; \mbox{min}\; \left\{\beta-\frac{b_I}{K}, \frac{K \sigma^2}{2}\right\}$. Integrating both sides of (4.5) from $0$ to $t$, dividing by $t$, and using (4.4), we get

Since $I$ is bounded, $\lim_{t \to \infty}\frac{\mbox{ln}\; I(0)}{t} = 0$ and by the argument given previously, $\lim_{t \to \infty}\frac{\sigma}{t}\int_0^t S(\tau) d\xi(\tau) = 0$. Therefore, using Lemma 1, we have

Thus, if the condition stated in the theorem holds then the infected population persists almost surely for all future time. □

4.2.   Asymptotic stability

It is worth mentioning that stochastic system (2.2) has no explicit equilibrium density, unlike the deterministic system (2.1). Instead, stochastic solutions fluctuate around the deterministic equilibrium or a fixed value. We provide below the conditions under which the behavior of the stochastic system is similar to that of the deterministic equilibrium solution.

Theorem 5. Let $(S(t), I(t))$ be the solution of (2.2) with initial value $(S(0), I(0)) \in \mathbb{R}_+^2$ and the stability criteria of $E^*$ hold. If $\frac{b_S}{K}(S^*+1)+g_1 \bar{L}+\frac{b_I I^*}{K}+u_S-\left(\frac{3}{2}b_s+e+b_I\right) > 0$, $\frac{e}{K}(S^*+I^*)+\frac{(b_S+b_I) S^*}{K}+\frac{b_I}{K}I^*+g_2 \bar{L}+ \frac{c_1 b_I}{K}+u_I-\left(\frac{b_S}{2}+2(e+b_I)\right) > 0$ and $\frac{2S^* b_S }{K}+\frac{e(S^*+I^*)}{K}+\frac{b_I(S^*+I^*)}{K}-b_S-b_I-e > 0$, then $\limsup_{t\rightarrow \infty}\frac{1}{t} \int_0^t [(S(\tau)-S^*)^2+(I(\tau)-I^*)^2] d(\tau) \leq G_1 \sigma^2 \; \; \mathit{\mbox{a.s.}},$ where $E^* = (S^*, I^*)$ is the endemic equilibrium of the deterministic system and $G_1 = \frac{K^2\left(K^2+\frac{c_1}{2}I^*\right)}{H}$, $H = \mathit{\mbox{min}}\; \bigg\{\frac{b_S}{K}(S^*+1)+g_1 \bar{L}+\frac{b_I I^*}{K}+u_S-\left(\frac{3}{2}b_s+e+b_I\right),$ $\frac{e}{K}(S^*+I^*)+\frac{(b_S+b_I) S^*}{K}+\frac{b_I}{K}I^*+g_2 \bar{L}+ \frac{c_1 b_I}{K}+u_I-\left(\frac{b_S}{2}+2(e+b_I)\right) \bigg\}$, $g_1 = \mathit{\mbox{min}}\left\{\frac{e}{K}, \frac{b_S}{K}\right\}$, $g_2 = \mathit{\mbox{min}}\left\{\frac{e}{K}, \frac{b_I}{K}\right\}$ and $\bar{L} = \frac{K}{M}(M-N)$.

Proof. Define a positive function $G = \frac{1}{2}(S-S^*+I-I^*)^2+ c_1\left(I-I^*-I^*\mbox{log}\; \frac{I}{I^*}\right)$, where $c_1$ is a positive constant to be determined later. At $(S^*, I^*)$, we have

Using (4.7) and Ito's formula, one obtains

Choosing $c_1 = \frac{1}{\beta-\frac{b_I}{K}} \bigg[\frac{2S^* b_S }{K}+\frac{e(S^*+I^*)}{K}+\frac{b_I(S^*+I^*)}{K}-b_S-b_I-e\bigg] > 0$ and using $|(S-S^*)(I-I^*)| \leq \frac{1}{2}\left((S-S^*)^2+(I-I^*)^2\right)$, (4.7) becomes

where $\bar{L} = \frac{K}{M}(M-N)$.

Define $H = \mbox{min}\; \bigg\{\frac{b_S}{K}(S^*+1)+g_1 \bar{L}+\frac{b_I I^*}{K}+u_S-\left(\frac{3}{2}b_s+e+b_I\right)$, $\frac{e}{K}(S^*+I^*)+\frac{(b_S+b_I) S^*}{K}+\frac{b_I}{K}I^*+g_2 \bar{L}+ \frac{c_1 b_I}{K}+u_I-\left(\frac{b_S}{2}+2(e+b_I)\right) \bigg\}$. Integrating (4.7) from 0 to $t$, one gets

Let $N_1(t) = \int_0^t (I(\tau)-I^*)S(\tau) d\xi(\tau)$, which is a continuous martingale and $N_1(0) = 0$. Also, $< N_1, N_1 > _t = \left(\int_0^t (I(\tau)-I^*)S(\tau) d\xi(\tau)\right)^2 = \int_0^t (I(\tau)-I^*)^2S^2(\tau) d(\tau) \leq 4 K^4 t$ and $\limsup_{t \to \infty} \frac{ < N_1, N_1 > _t}{t}\leq 4K^4 < \infty \; \mbox{a.s}$. Therefore, by the law of large numbers ^[35], $\lim_{t \to \infty} \frac{N_1(t)}{t} = 0 \; \; \mbox{a.s}$. Combining these results and then dividing (4.7) by $t$ and taking limit superior, we have

Hence, the theorem is proven. □

Furthermore, if $\sigma \rightarrow 0$, then

Therefore, $\lim_{t\rightarrow \infty}S(t) \rightarrow S^*, \; \lim_{t\rightarrow \infty}I(t) \rightarrow I^*$, and the stochastic solution tends to the deterministic equilibrium solution. This implies that if the noise intensity is low, the stochastic system behaves similarly to the asymptotic solution of the deterministic system, provided the restrictions in the above theorem hold.

4.3.   Probability density function

In this part, we compute the exact expression of the density function of the stationary distribution. The $S$-compartment of deterministic model (2.1), under usual assumption of constant population $S(t) + I(t) = N$, reads

Under the stochastic perturbation of $\beta$ to $\beta + \sigma \xi(t)$, the above equation becomes

In the case of white noise perturbation with $\xi(t)$, the probability density function $p(S, t)$ of $S(t)$ is governed by the classic Fokker-Planck equation

where, $h(S) = b_S S \left(1-\frac{N}{K}\right) -u_S S - \beta S (N-S) + e(N-S)\left(1- \frac{N}{K}\right) + \mu(N - S)$ and $\sigma(S) = - \sigma S (N - S)$. The closed-form stationary solution $p_0(S) = \lim_{t \to \infty} p(S, t)$ of the Fokker-Planck equation is given by $p_0(S) = \frac{\mathcal{N}}{\sigma ^2(S)} \mbox{exp}\left(2 \int_{}^{S}\frac{h(\omega)}{\sigma^2(\omega)}d\omega \right),$ where $\mathcal{N}$ is the normalizing factor. Upon simplification, the closed-form expression is given by

For the $I$ compartment, we similarly assume that $q(I, t)$ is the pdf that is governed by the Fokker-Planck equation

where $h'(I) = b_I I \left(1 - \frac{N}{K}\right) - u_I I + \beta (N - I) I - \mu I$ and $\sigma'(I) = \sigma (N-I)I$. The closed form stationary solution is given by

where

and $\mathcal{N'}$ is the normalizing factor. Simplifying the above expression, we obtain

5.   Simulation results

We present simulation results of stochastic system (2.2). The following system parameters are considered fixed unless it is stated: $b_S = 0.2, \; b_I = 0.03$, $u_S = 0.1, \; u_I = 0.15, \; e = 0.02, \; {K = 1}, \; \beta = 0.5, \; \mbox{and}\; \mu = 0.1 \;$ ^[9]. Without loss of generality, we set $K = 1$, effectively selecting an appropriate spatial unit for measuring population densities, as outlined in Lipsitch et al. ^[9]. The initial value is assumed to be $(S(0), I(0)) = (0.6, 0.3).$

5.1.   Time evolution

The time evolution of a stochastic system is different for each run. However, that of the deterministic system is unique for a given parameter set. Thus, due to inherent randomness, a considerable difference may exist between the solutions of the stochastic system for different runs. To illustrate this, we have presented $10$ simulation results for two distinct noises in the first two rows of Figure 1. The first column shows that the disease persists for low values of the noise ($\sigma = 0.4$) when the parameter set satisfies the conditions of Theorem 4 with $\frac{b_I-u_I-\mu}{\frac{b_I}{K}+\frac{A(b_I+e)}{b_S}-\frac{A u_I}{b_S}} = 22 > 0$. The second column illustrates that the disease dies out for higher values of the noise ($\sigma = 0.8$) when the parameter set satisfies the conditions provided by Theorem 3 with $R_0^S = 0.9012 < 1$. The solution trajectories here differ in each simulation, though the initial condition and parameter values remain the same. It is, therefore, more justified to present the average behavior of the stochastic solutions. In the last row, we have given the mean value of $100$ solutions of stochastic system (2.2) for $\sigma = 0.4$ ((a), (b), (e)) and $\sigma = 0.8$ ((c), (d), (f)), which show similar behaviors.

5.2.   Steady-state

It is worth mentioning that the noise-induced stochastic system (2.2) has no steady-state value. In Theorem 5, we have proven that the stochastic and deterministic systems behave similarly if the noise intensity is low. However, if the noise intensity is high, the stochastic solution largely deviates from its deterministic counterpart. We illustrate such behavior of the system in Figure 2 for two noises: $\sigma = 0.05$ and $\sigma = 0.2$. The conditions of the stationary distribution (Theorem 5) are satisfied here.

We further examine the probabilistic smoke cloud around the deterministic steady state. We repeat the simulation $10,000$ times and plot the values of $S$ and $I$ populations at $t = 120$ in Figure 3. This figure shows that species densities are primarily concentrated in the neighbourhood of deterministic concentrated in the neighborhood of deterministic equilibrium value $(S^*, I^*) = (0.47, 0.014)$ (blue point) if the noise intensity is low $(\sigma = 0.05)$. Observe that the frequency distribution of the susceptible population ($S$) is distributed in the range $0.45-0.48$ around its equilibrium value $S^* = 0.47$ (Figure 3(b)), and the same for the infected population ($I$) is distributed in the range $0.005-0.025$ around its equilibrium value $I^* = 0.014$ (Figure 3(c)). This shows the strength of the stabilizing factors of the population interaction compared to the diffusive effects of the random environmental fluctuations ^[36].

Indeed, the probability cloud is compact if the interaction strength can overcome the noise. However, if the interaction strength is weak compared to the intensity of the environmental variance, the cloud covers a larger area (Figure 3(d)) if the noise intensity is high $(\sigma = 0.15)$. The corresponding frequency distributions (Figures 3(e), (f)) show that $S$ and $I$ are distributed over a more extensive range of $0.41-0.5$ and $0.002-0.045$, respectively, around their deterministic equilibrium values. This reveals that the stochastic system does not deviate too much from its deterministic system if the noise is low, but it does if it is high.

We further examined the evolution of histograms of $S$ and $I$ populations, as an approximation of their probability density functions. We repeat the simulation 50,000 times to plot the histograms. As time progresses, one can observe a shift in the mean and spread of both $S(t)$ and $I(t)$. At time $t = 20$, the distribution of $S(t)$ is centered around 0.31 and $I(t)$ around 0.07, indicating that a smaller proportion of the population is infected (Figure 4(a)). At $t = 40$, the distribution is shifted to around 0.40, and $I(t)$ has a peak around 0.02. This may indicate an increase in the susceptible individuals as the infected population temporarily decreases, possibly due to recovery or a reduction in transmission (Figure 4(b)). At $t = 80$, $S(t)$ has a mean near 0.46, and $I(t)$ near 0.01, with both distributions showing a greater variability (Figure 4(d). As time progresses, $S(t)$ tends to fluctuate near 0.48 and $I(t)$ near 0.01, which indicates that the system is approaching a stationary distribution, even though individual realizations of $S(t)$ and $I(t)$ still experience random fluctuations. This is useful for analyzing long-term behavior in stochastic models, as it suggests that the system has an average state around that it oscillates due to random environmental noise.

5.3.   Disease extinction time

It is observed that the noise intensity ($\sigma$) and the disease transmission rate ($\beta$) play a pivotal role in the extinction and persistence of the disease. We present, in Figure 5(a), the mean extinction time of the infected population of the stochastic system (2.2) with the variations in $\beta$ and $\sigma$. To find the extinction time numerically, we repeat the simulation $100$ times for each pair of ($\beta$, $\sigma$) and plot the mean of the extinction time of $100$ simulations. The infected population is considered extinct if it goes below $0.0001$. Moreover, $R^S_0 < 1$ is satisfied throughout the considered range of $\beta$ and $\sigma$. Figure 5(a) shows that the extinction time is longer if $\beta$ is high and $\sigma$ is low. However, the extinction time decreases if $\beta$ is low or $\sigma$ is high. We also investigate how disease extinction time is influenced by the parameter variations in $\beta$ and $b_I$. We select four values of $\sigma$, namely, $0.2, \; 0.4, \; 0.6, \; 0.8$. For each of these $\sigma$ values, we jointly vary $\beta$ and $b_I$ and record the mean extinction time for 100 simulations (see Figure 5(b)–(e)). It has been observed that disease extinction takes a shorter time if both $\beta$ and $b_I$ are low. However, disease extinction time is prolonged for higher values of $\beta$ and $b_I$. We further observe that when $\sigma$ is low, disease extinction time in the whole parametric space of $\beta$ and $b_I$ is high compared to the higher value of $\sigma$. This reveals that stronger environmental noise helps eradicate the disease from the system early, implying that it may act as a regulatory mechanism to control the disease.

To unveil the underlying law of extinction time, we plot the mean of extinction time of $100$ simulations with the variation of $\sigma$ for three particular values of $\beta$: $0.1, 0.2$, and $0.3$ (Figure 6(a)). Mean extinction time gradually decreases as $\sigma$ increases. An exponential curve fits it. A similar figure was drawn with three fixed values of $\sigma$, $0.2, 0.5$, and $0.8$, for the variation of $\beta$ (Figure 6(b)). In the latter case, the disease extinction time gradually increases with the increase of the disease transmission coefficient. The mean extinction time fits the positive exponential curve.

It is shown that $R^S_0$ is a decreasing function of the recovery rate, $\mu$. Here, we explore the disease extinction scenario for $\mu$ for different values of $\sigma$ and $\beta$. Figure 7 shows that the extinction time curves have two distinct natures for the variations in $\mu$. If $\sigma$ is much higher than $\beta$ (cyan curve, bottom), the extinction time follows the negative exponential law. As the difference between $\sigma$ and $\beta$ decreases while $\sigma$ remains higher (purple and green curves), the extinction time increases. On the other hand, if $\sigma$ and $\beta$ are both elevated and nearly equal in magnitude, the extinction time becomes substantially larger, exhibiting a slower rate of decline with respect to $\mu$. All these have a long tail, implying that a shorter extinction time of the disease is possible for an extended range of $\mu$.

5.4.   Parasite fitness

As mentioned earlier, parasite infectivity and the relative fecundity of infected hosts are critical for parasite fitness. The survival fitness of parasites depends on how long the infected hosts persist in the system. For this, we plot (Figure 8) the extinction time of infected hosts from the system for the simultaneous variations in $\beta$ and $\frac{b_I+e}{b_S}$. Recall the ensemble fecundity of the infected population as $b_{I}+e$, and that of the susceptible as $b_S$. The ratio $\frac{b_I+e}{b_S}$ denotes the relative fecundity of the infected host population. By our assumption, $b_S\geq b_I+e$ and therefore, the ratio $\frac{b_I+e}{b_S}$ is always less than or equal to $1$. If the force of infection is not too high ($\beta < 0.6$), then the infection eradication time increases with the increasing relative fecundity of the infected host for some fixed value of noise. However, for the higher force of infection ($\beta > 0.6$), the relative fecundity of infected hosts has a negligible effect on extinction time. The trend remains the same if the noise intensity is altered. Thus, the parasite fitness is low if the relative fecundity and disease transmissibility are low but gradually increase if any of them increase.

5.5.   Sensitivity analysis

5.5.1.   Local sensitivity analysis

To assess how variations in parameters control disease transmission and progression with pseudo-recovery, a sensitivity analysis of model (2.1) is conducted numerically, following the approach outlined in ^[37,38]. The normalized forward-sensitivity index of a variable, $v$, that depends differentiably on a parameter, $p$, is defined as

We numerically evaluate the local time-dependent normalized sensitivity index for the variables $S$ and $I$. The positive sensitivity of $S$ with respect to $\mu$ and negative sensitivity of $I$ is expected, as the higher recovery rate would naturally shift more individuals from the infected state back to the susceptible state. The oscillatory behavior here might be due to the model's feedback between $S$ and $I$, where changes in recovery (through $\mu$) indirectly affect $S$ and $I$ by influencing the rate at which individuals transition between susceptible and infected states.

5.5.2.   Global sensitivity analysis

In a sensitivity analysis based on variance, such as the Sobol method, the total variance in the model output is decomposed to attribute portions of this variance to different input parameters ^[39]. Let $Y$ represent the model output (either $S$ or $I$) and $\theta$ be a vector of model parameters ($\beta, b_i, b_s, e, K, \mu, u_i, u_s$). The variance of $Y$, $\mathrm{Var}(Y)$, can be decomposed as:

where $\mathrm{Var}_{\theta_i}(Y)$ is the partial variance due to parameter $\theta_i$ alone, while $\mathrm{Var}_{\theta_i, \theta_j}(Y)$ represents the variance due to the interaction between $\theta_i$ and $\theta_j$, and so on. The first-order Sobol index $S_i$ for each parameter $\theta_i$ is defined as

This index quantifies the direct effect of parameter $\theta_i$ on $Y$.

In the above pie chart, each segment represents the proportion of first-order variance in $S$ or $I$. A larger segment implies a higher influence of that parameter on the respective model output. In the first-order Sobol index of $S$, we observe that the parameters $\beta$ and $u_S$ strongly influence the variance of $S$, indicating that variation in these parameters significantly affects the susceptible population. In the right panel, $u_S$ and $u_I$ have the highest impact on the infected population.

6.   Case study

Here, we present a case study for the imperfect vertical and horizontal infection transmission that supports some of our observed theoretical results of the deterministic and stochastic systems. The parasite Holospora Undulata frequently infects the protozoa Paramecium caudatum and uses both the horizontal and vertical transmission pathways ^[40]. P. caudatum are filter feeders and ingest infectious forms of Holospora Undulata while taking food from the water. Lunn et al. ^[41] performed a microcosm study, where the 34 days time series data of P. caudatum and its specialist bacterial parasite Holospora Undulata under different nutrient concentrations (high and low) have been investigated. We fitted our SIS deterministic and stochastic models with these 34 days of experimental data corresponding to high nutrient concentrations following the optimization technique presented in Appendix Ⅰ. During the estimation of the parameter set, we maintain the model assumptions $u_I\geq u_S$ and $b_S\geq b_I+e$. The experimental values are plotted on a logarithmic scale for ease of fitting. The best-fit parameter set thus obtained by fitting the stochastic system (2.2) is as follows:

The deterministic basic reproduction number $R_0^{(1)}$ for this parameter set is $4.1801$, implying the disease is persistent. The average value (solid line) of 1000 solutions of the stochastic system (2.2) with the above parameter set has been plotted in Figure 11 with respect to time (day) with $95\%$ confidence interval (CI). This shows that the stochastic solution matches the experimental data except at some initial data points.

7.   Discussion

Mathematical models of epidemics have helped understand the disease spread in a population. After the benchmark work of Kermack and McKendrick ^[42], many mathematical models have been proposed and analyzed, considering the epidemiological demands. Most of these models are deterministic types represented by ordinary differential equations. Deterministic models have been criticized for viewing all the rate parameters as constant, even though these constants constantly fluctuate due to the various unknown environmental noises. Uncertainty is an integral part of modeling biological phenomena due to the complexity and lack of knowledge of microscopic events involved in the system. In particular, the disease transmission phenomena are random, and the transmission rate depends on various factors. Incorporating such randomness in a model leads to stochastic differential equation models.

In this paper, we consider an SIS-type stochastic epidemic model, where the disease is transmitted through horizontal and vertical transmission modes. To incorporate stochasticity in the model, we introduce a white noise in the horizontal disease transmission term, the predominant disease transmission mechanism. The fluctuation in the transmission rate usually occurs around some mean value. Therefore, the error term follows a normal distribution, enabling us to approximate it by a white noise ^[43].

Eradication or control of various emerging and reemerging diseases is a global challenge. Environmental factors like temperature, humidity, rainfall, and pollution significantly affect disease persistence. Knowing the different routes for theoretically eradicating an infection from the population is essential. Our mathematical analysis revealed that the population can be disease-free by modulating some parameters and the noise intensity. We have shown that if the stochastic basic reproduction number, $R^S_0$, can be made less than unity, the disease can be eradicated from the system. It is observed that noise intensity plays a pivotal role in eliminating and has an inverse relationship with $R^S_0$. The environmental noise can make a system disease-free if it significantly affects the transmission mechanism and makes $R^S_0 < 1$. On the contrary, $R^S_0$ is directly proportional to the disease transmission coefficient. Thus, $R^S_0$ may exceed the threshold value if the disease transmissibility increases. Therefore, the management strategy should be to reduce horizontal disease transmissibility. Our numerical computations revealed that the disease eradication time follows a negative exponential law with increasing noise intensity. In contrast, the disease transmission coefficient follows the positive exponential law.

Furthermore, disease eradication is also possible with respect to controllable parameter $\mu$. This parameter measures the recovery rate of the infected population and can be used in the disease eradication process. Moreover, an infected host may recover from the infection due to its immune mechanism. Recovery of infected hosts is also possible with the help of some external measures. Extinction time follows either a negative exponential or Gaussian law for some given noise and disease transmissibility values. In the former case, the noise must be high, or the transmission coefficient must be low. In the latter case, the effects of noise and disease transmissibility should be reversed. The eradication time will be significantly quicker if the recovery rate is high. On the other hand, disease management would be more challenging if the noise intensity is low and the disease transmissibility is high.

Several researchers have demonstrated that parasites have a detrimental effect on the host's fecundity. The survival of parasites depends on the infected host density and the extinction time of the infected host population. Thus, parasites are expected to evolve towards higher relative fecundity of the infected host to enhance fitness. Our simulation result for simultaneous variation of the disease transmission coefficient and relative fecundity rate shows that parasite fitness increases with increasing relative fecundity of parasites when the transmissibility is low. The relative fecundity of infected hosts, however, has a negligible effect on the extinction time of the infected host population and, consequently, on the parasite survival fitness if the disease transmissibility of parasites is high. A similar trend is also maintained if the noise intensity is altered.

We utilizie a susceptible-infected-susceptible framework, which is appropriate for diseases where recovered individuals do not gain permanent immunity and can reenter the susceptible population. However, many infectious diseases exhibit immunity following infection, aligning more closely with the SIR or SIRS model frameworks, where recovered individuals gain immunity before potentially returning to the susceptible class. Integrating stochastic effects and vertical transmission in these models could reveal important nuances in epidemic behavior under varying immunity dynamics with the increasing complexity of mathematical treatment ^[16,44,45]. One of the central limitations of the white noise perturbation is that increasing the white noise leads to quicker disease extinction, which is typically not observed in real-world epidemics. Therefore, the white noise model can underestimate the severity of disease transmissions ^[46,47]. This is also counterintuitive because the increase of noise in the transmission parameter can be a relaxation of the restrictions of contact among species. Thus, the increased noise that could lead to rapid disease eradication might be mathematically correct, but the real-world implications of such a conclusion must be carefully examined. However, temporally correlated noise (modeled via the Ornstein-Uhlenbeck process) introduces more realistic fluctuations in transmission rate and does not necessarily foster disease extinction due to increased stochastic noise ^[48,49].

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is founded by the Indo-French Centre For Applied Mathematics (IFCAM), Ref No. MA/IFCAM/18/49.

Conflict of interest

The authors declare there is no conflicts of interest.

Appendix

Deterministic model parameters were estimated from the experimental data set through the nonlinear least square technique. To begin, we considered an arbitrary parameter set ${\Lambda} = (\Lambda_1, \Lambda_2, ...., \Lambda_n)$ and then minimized the sum of squared differences between the experimental data $Y = (y_1, y_2, ..., y_n)$ and the corresponding model solution data. We used lsqcurvefit, an optimization toolbox from Matlab, to solve the nonlinear least square problem. Lsqcurvefit fits a function to a set of data points using the Levenberg-Marquardt (LM) algorithm, which is an iterative process used to find the minimum of the sum of squares of the nonlinear function ^[50]. It can be understood as an amalgamation of steepest descent and Gauss-Newton methods ^[51]. The assumed objective function is defined as

where $g(t_j, \Lambda)$ is the model solution at time step $t_j$ at which the solution is to be estimated, and $n$ is the total number of data points. The optimization process will provide the best fit solution $g(t_j, \tilde{\Lambda})$ to the experimental data corresponding to the optimal parameter set $\tilde{\Lambda} = (\tilde {\Lambda}_1, \tilde{\Lambda}_2, ...., \tilde {\Lambda}_n)$, which should satisfy the model's basic assumptions. After estimating parameters for the deterministic system, we estimated the noise intensity for the stochastic system to obtain a good agreement between the output and the experimental data. We computed the corresponding $r$-squared value for stochastic simulation data and experimental data to find the suitable noise strength. Here, we considered 10,000 different random values of the noise intensity $\sigma$ between 0 and 1 through Latin hypercube sampling. Then, for each of these $10,000$ values of $\sigma$, the stochastic system was simulated 100 times. Taking mean of the 100 simulations, the $r-$ squared was computed between the average stochastic simulation output and the experimental data. Noise intensity was selected for which the $r$-squared was maximum.