Determining the global threshold of an epidemic model with general interference function and high-order perturbation

1.   Introduction

Mathematical formulations play a central role in exploring and predicting the future of the communicable diseases ^[1]. For example, in the case of the coronavirus disease and of course its new mutations, many researchers have contributed to modeling the mechanisms of its spread and providing scientific recommendations to control its expansion ^[2,3,4,5]. It should be noted that most of the epidemiological models presented in these articles are an improved and adapted version of the SIR model ^[6]. This compartmentalized system was constructed by Kermack and McKendrick ^[7] and its philosophy is based on the idea of subdividing individuals according to their different characteristics. Explicitly, individuals are grouped into three main groups: susceptible class (${\bf C}_1$), infected class (${\bf C}_2$) and permanently recovered class (${\bf C}_3$). The transmit rates among these groups are determined by the following dynamical system:

where $\Theta > 0$ designates the inflow rate into ${\bf C}_1$, $\mathfrak{u} > 0$ indicates the normal death rate, $\mathfrak{b} > 0$ is the dissemination rate of the epidemic, $\mathfrak{a} > 0$ is the mortality ratio due to the infection, and $\mathfrak{c} > 0$ is the cure rate. We note that the above system is one of the most straightforward epidemiological models used to depict the first wave of COVID-19 ^[8]. Lately, Zhou et al. ^[9] proposed a general version of system (1.1) by including the pre-existing immunity presumption and the nonlinear incidence function $\mathfrak{b}{\bf C}_1(t)g\big({\bf C}_2(t)\big)$. By mandating certain conditions on $g$, they introduced the following system:

where $\mathfrak{z}$ is the pre-existing immunity rate, and $\mathfrak{u}_1$, $\mathfrak{u}_2$, $\mathfrak{u}_3$ are respectively the normal mortality rates of ${\bf C}_1$, ${\bf C}_2$, ${\bf C}_3$. The function $g$ covers some functional responses, for instance, $g\big({\bf C}_2\big) = {\bf C}_2$ ^[10], $g\big({\bf C}_2\big) = \frac{{\bf C}_2}{\mathfrak{m}+{\bf C}_2}$ ($\mathfrak{m} > 0$) ^[11] and ${ }\mathfrak{b}g\big({\bf C}_2(t)\big) = \mathfrak{b}-\frac{\mathfrak{b}_c}{\mathfrak{m}+{\bf C}_2(t)}$, where $m > 0$ is the media intrusion rate and $\mathfrak{b}_c > 0$ is the reduced active contact rate ^[12].

When dealing with epidemiological models, more characteristics can be considered such as cross-individual overlap ^[13]. It is worth to point out that the choice of functional response affects the prediction of illness behavior. Moreover, the previous setup of Zhou et al. ^[9] overlooks a large category of incidence rates that are often used in the literature. In this research, we exhibit an enhanced SIR illness system with an interference function that includes additional response examples (see Table 1). In line with this setting, the system (1.2) can be rewritten as follows:

We presume that the general interference response ${\mathcal H}\in \mathcal{C}^2({\mathbb R}_+\times{\mathbb R}_+, {\mathbb R}_+)$ verifies these two hypotheses:

$\bullet$ Assumption (a): There exists a constant $\varrho > 0$ such that for all ${\bf C}_1, {\bf C}_2 \geq 0$,

$\bullet$ Assumption (b): $\underset{{\bf C}_2\to 0}{\lim}\underset{{\bf C}_1 > 0}{\sup}\{|{\mathcal H}({\bf C}_1, {\bf C}_2)-{\mathcal H}({\bf C}_1, 0)|\} = 0$.

The properties $(a)$ and $(b)$ are readily satisfied by the typical examples listed in Table 1.

Substantially, the infection mechanism is random in nature at all scales ^{[16,17,18,19,20]}. From the attachment of the virus to the human cell to the encapsulation of repetitive genetic information, from the shedding of new virions to the transmission of a second individual, from individual behavior to global mobility, all the factors influence the diffusion of the infection and make it more uphill task to foretell its conduct ^[21,22,23]. The stochastic approach offers considerable advantages in providing insight into population dynamics under the said fluctuations ^[24,25,26]. To correctly describe this randomness, a series of perturbed compartmental models with different assumptions that simulate reality have been proposed ^{[27,28,29,30,31]}. Most of these models assume that random fluctuations can be modeled by integrating Brownian motions into the deterministic formulation ^[32]. By selecting this class of fluctuations in their linear shape, probabilistic systems are widely used in epidemiology to analyze disease prevalence ^[33,34,35]. In these works, the primary focus was the investigation of some biological long-run characteristics of the infection.

Since 2017, a new form of stochastic systems has emerged where the story begins with the work proposed by Liu and Jiang in ^[36]. By reason of the complexity of environmental changes, they claimed that the relative linear order of the disturbance can be raised to the second. Based on this assumption, scientific papers have proposed and analyzed various real systems with second-order fluctuations ^[37,38]. Recently, the autrhors of ^[9] suggested an enhanced type of perturbation in its general representation. This frame generalizes the previously mentioned studies and offers a new line of research. Motivated by their arguments, we consider the following polynomial perturbed system:

where ${\mathbb W}_1(t)$, ${\mathbb W}_2(t)$, ${\mathbb W}_3(t)$ are independent Brownian motions defined on a filtered probability space $\Omega^{\mathcal{E}, {\mathbb P}}\equiv(\Omega, \mathcal{E}, \{\mathcal{E}_t\}_{t\geq 0}, {\mathbb P})$ such that $\{\mathcal{E}_t\}_{t\geq 0}$ follows the usual assumptions, and $\mathfrak{q}_{k h} > 0$ $(k = 1, 2, 3, \; h = 0, 1, 2, \dots, N)$ are the high-order intensities of white noises.

In ^[9], the authors indicated that the model (1.4) in the case of incidence $\mathfrak{b}{\bf C}_1(t)g\big({\bf C}_2(t)\big)$ and polynomial perturbation is well-posed mathematically and biologically, then they treated its long-run behavior. The problem is that they obtained two separate critical conditions for extinction and stationarity which is not ideal when addressing epidemiological models. Moreover, {they mention that there was} a large gap between the defined criteria; and the corresponding threshold value is still an open question (for more details see the discussion part of ^[9]). Compared to the results presented in ^[9], in this article, we address the said problem from a global angle by considering an epidemic model with a general incidence function. Taking the latter into account makes the analysis very complex, which has prompted us to innovate alternative techniques. By adding the polynomial perturbation, we present the acute threshold value between stationarity and extinction of the infection, which offers an excellent insight into the possible scenarios of epidemic status in a given population.

Technically, we introduce an analytical method based on some long-term characteristics of an auxiliary boundary equation ^[39]. By using the ergodic theorem and the stochastic comparison lemma, we establish the well-defined threshold between stationarity and infection extinction. Our method differs from the one used in ^[9] by using the mutually exclusive possibilities lemma and other analytical tools. Specifically, we focus on the long-term characteristics of the Markov process ${\bf D}(t)$ that verifies

In accordance with Lemma 5 of ^[9], ${\bf D}$ admits the following single stable distribution:

where $C^{{\bf D}}$ is a specific constant that verifies ${ }\int_{{\mathbb R}_+}\pi^{{\bf D}}(y){\rm{d}} y = 1$. From the probabilistic comparison result ^[40], we can compare the processes ${\bf D}$ and ${\bf C}_1$ as follows: ${\bf D}(t)\geq {\bf C}_1(t)$ almost surely (a.s.). Furthermore, the time average of ${\bf D}(t)$ converges almost surely to ${ }\int_{{\mathbb R}_+}y\pi^{{\bf D}}({\rm{d}} y)$ as $t\to \infty$. In accordance with the above results, we clearly state that the present work aims to prove that the following quantity

is the sharp threshold between the disappearance of the illness and the ergodic characteristic of the system (1.4), and its sign provides a stellar overview of the potential scenarios of the epidemic situation. In other words, this article proposes a nice generalization of the article ^[9] and presents a new treatment applicable to other complex models. Importantly, we show numerically that the extreme amount of disturbance reduces the disease extinction time.

The rest of the paper is organized as follows: In Section 2, we show that the disappearance scenario occurs when $\mathcal{T}_{\circ}^{\Sigma} < 0$. In Section 3, we prove that the scenario of the ergodicity of system (1.4) occurs when $\mathcal{T}_{\circ}^{\Sigma} > 0$. In Section 4, we verify numerically the correctness of our outcomes. Finally, we conclude this paper in Section 5 with a discussion.

2.   Scenario 1: disease disappearance

This section aims to exhibit the criterion for the demise of the infection.

Theorem 2.1. The disappearance of the disease occurs if $\mathcal{T}_{\circ}^{\Sigma} < 0$.

Proof. By employing Itô's lemma for drift-diffusion processes ^[41], we obtain

In line with the probabilistic comparison lemma, we conclude that

After that, we make two operations on both sides of (2.1): integration from $0$ to $t$ and division by $t$, then the result is

The next step is based on the use of the exponential inequality for martingales ^[41], which leads to

for all $0 < \alpha_1 < 1$ and $T_1 > 0$. From the Borel-Cantelli result ^[41], we assure the existence of $T_{1, \omega} = T_1(\omega)$, $\forall \omega$ in $\Omega$, such that

holds for all $T_1\geq T_{1, \omega}$ and $T_1-1 < t\leq T_1$ a.s. Therefore,

By taking the limitsup on two sides of (2.2), we infer that

We let $\alpha_1$ tends to $0^{+}$, then the obtained result is

That implies the stochastic extinction $\underset{t\to\infty}{\lim} {\bf C}_2(t) = 0$ a.s. In other words, the class of individuals carrying the infection will disappear.

3.   Scenario 2: the ergodicity of the model (1.4)

This section introduces a new approach to establish the criterion of the ergodicity of our probabilistic system (1.4).

Theorem 3.1. The single ergodic stable distribution of the probabilistic model (1.4) exists if $\mathcal{T}_{\circ}^{\Sigma} > 0$.

Proof. In order to reduce notations and provide clear mathematical writing, we set

The Itô differential operator $\mathcal{L}$ associated with the stochastic equation of ${\bf C}_2(t)$ is given by

From Assumption (a), we have

Now, the application of $\mathcal{L}$ on $\big(\ln {\bf D}(t)-\ln {\bf C}_1(t)\big)$ gives

Since ${\bf D}(t)\geq {\bf C}_1(t)$ a.s., we obtain

We define the function $\Phi(t)$ as follows:

Based on (3.1) and (3.2), we get

We add and subtract at the same time the quantity ${ }\mathfrak{b} \int_{{\mathbb R}_+}{\mathcal H}(y, 0)\pi^{{\bf D}}({\rm{d}} y)$ in (3.3) as follows:

To eliminate the term associated with ${\bf C}_2(t)$, we set

where the positive constant $\Theta$ verifies

The application of $\mathcal{L}$ on $\Phi_{\Theta}(t)$ gives

In the same vein, we apply $\mathcal{L}$ on the function $\zeta^{-1}\big(1+{\bf C}_1(t)\big)^{\zeta}+\zeta^{-1}{\bf C}_2^{\zeta}(t)$, $\zeta\in(0, 1)$, then

Accordingly, we derive that

Now, we define a new function $\Phi_{\Theta, \zeta}$ as follows:

where $S > 0$ satisfies that $-S\mathcal{T}_{\circ}^{\Sigma}+\mathcal{Z}+2\leq 0$ and $\mathcal{Z}$ is given by

Clearly, the function $\Phi_{\Theta, \zeta}$ reaches its minimum value at a point $({\bf C}_1^{\ell}, {\bf C}_2^{\ell})$. For this reason, we will consider a new non-negative function defined as follows:

From the above calculation, we obtain

Now, we define five sets:

Here, ${\mathbb R}_+^{2, \star} = \{(x, y):\; x > 0, y > 0\}$, $a_{\star} = \min\{a_{\circ}, a\}$, where $a_{\circ} > 0$ verifies (3.9), and $a > 0$ is chosen carefully such that

where

Plainly, $\mathcal{J}_{a, a_{\star}}^c = {\mathbb R}^{2, \star}_{+}\setminus\mathcal{J}_{a, a_{\star}} = \mathcal{J}_{a, 1}\cup\mathcal{J}_{a_{\star}, 2}\cup\mathcal{J}_{a, 3}\cup\mathcal{J}_{a_{\star}, 4}$. In the following, we will verify that

for any $\big({\bf C}_1(t), {\bf C}_2(t)\big)\in\mathcal{J}_{a, a_{\star}}^c$ which is equivalent to showing it on $\mathcal{J}_{a, 1}$, $\mathcal{J}_{a_{\star}, 2}$, $\mathcal{J}_{a, 3}$ and $\mathcal{J}_{a_{\star}, 4}$, respectively. For this reason, we have the following situations:

$(1)$ Assume that $\big({\bf C}_1(t), {\bf C}_2(t)\big)\in \mathcal{J}_{a, 1}$. From (3.4), we obtain

$(2)$ Here, we use the uniform continuity at ${\bf C}_2 = 0$ of the function ${\mathcal H}\big({\bf C}_1(t), {\bf C}_2(t)\big)$. By Assumption (b), $\exists a_{\circ} > 0$ such that as $0 < {\bf C}_2\leq a_{\circ}$,

Consequently, if ${\bf C}_2 < a_{\star} = \min\{a_{\circ}, a\}$, we get from (3.5) that

$(3)$ Assume that $\big({\bf C}_1(t), {\bf C}_2(t)\big)\in \mathcal{J}_{a, 3}$. From (3.6), we have

$(4)$ Assume that $\big({\bf C}_1(t), {\bf C}_2(t)\big)\in \mathcal{J}_{a_{\star}, 4}$. From (3.7), we get

In summary, the assertion (3.8) is obtained. On the other hand, we can easily show that $\exists \mathcal{O} > 0$ such that $\Psi\big({\bf C}_1, {\bf C}_2\big)\leq \mathcal{O}$, for all $\big({\bf C}_1, {\bf C}_2\big)\in {\mathbb R}_+^{2, \star}$. Accordingly, we get

By using the ergodic property of ${\bf D}(t)$, we conclude that

Consequently,

Hence,

where ${\mathbb R}_+^{3, \star} = \{(x, y, z):\; x > 0, y > 0, z > 0\}$. From Lemma 3.2 of ^[42] and also the mutually exclusive possibilities lemma ^[43], we confirm the existence of a single invariant distribution $\pi^{\Sigma}$.

Remark 3.1. From Theorem 3.1 of this paper and also Theorem 2.6 of ^[44], we can deduce interesting indications on the stochastic permanence of the Markovian processes ${\bf C}_1$, ${\bf C}_2$ and ${\bf C}_3$. Explicitly, we obtain that

By way of illustration, this indicates the persistence of the infection over time.

4.   Numerical verification

In this section, we exhibit some simulations to shed some light on the exactitude of our global threshold. For that purpose, we present three situations of system (1.4), and in each case, we explore the complex long-run behavior of the illness. We will consider general saturated interference function introduced in Table 1. The model parameters are theoretically selected according to well audit the outcomes of this paper. By using the high-order discrete Milstein method, the associated discretization equations of our system are directly obtained as follows:

where the time increment $\Delta t > 0$. $\xi_k$, $\zeta_k$ and $\alpha_k$ are three independent random variables which follow the Gaussian distribution $\mathcal{N}(0, 1)$. $({\bf C}_{1, k}, {\bf C}_{2, k}, {\bf C}_{3, k})$ is the corresponding value of the $k$-th iteration.

4.1.   Linear disturbance case ($N=0$)

In this example, we will deal with the following stochastic system:

associated with the following auxiliary process:

We choose $\Theta = 0.23$, $\mathfrak{u}_1 = 0.2$, $\mathfrak{u}_3 = 0.19$, $\mathfrak{z} = 0.02$, $\mathfrak{u}_2 = 0.2$, $\mathfrak{m}_1 = 0.1$, $\mathfrak{m}_2 = 0.1$, $\mathfrak{a} = 0.2$, $\mathfrak{c} = 0.02$, $\mathfrak{q}_{10} = 0.11$, $\mathfrak{q}_{20} = 0.112$ and $\mathfrak{q}_{30} = 0.1$. {By setting} $\mathfrak{b} = 0.4$ and considering a large time $T$, we get

From Theorem 2.1, we can infer the disappearance of the illness. Numerically, it can be seen in Figure 1 that the disease disappears after 230 days in the population with a strong permanence of classes 1 and 3. Now, we select $\mathfrak{b} = 0.55$ to switch from the case of extinction to the case of persistence. Then,

In accordance with Theorem 3.1, we infer that the properties of stationarity and ergodicity hold. From Figure 2, we offer a good illustration of these two statistical characteristics. Clearly, in this situation, the continuation of all classes is strongly happening which is depicted in Figure 3.

4.2.   Quadratic disturbance case ($N=1$)

In this example, {we deal with} the following stochastic system:

associated with the following auxiliary process:

For the comparison objective, we keep the same parameter values as the first case and we select $\mathfrak{q}_{11} = 0.022$, $\mathfrak{q}_{21} = 0.013$ and $\mathfrak{q}_{31} = 0.011$. As the above case, we select $\mathfrak{b} = 0.4$. Then, we get

Apparently, the condition of Theorem 2.1 holds. Numerically, the extinction phenomenon of the illness is illustrated in Figure 4. Classes 1 and 3 still persist. Now, we choose $\mathfrak{b} = 0.55$ to move from the case of extinction to the case of persistence. Then,

Based on Theorem 3.1, The last value implies that the system (4.2) admits a stable distribution. From Figure 5, we get an insight into the stationarity of the model (4.2). Furthermore, we offer Figure 6 to clarify the continuation of all classes ${\bf C}_1$, ${\bf C}_2$ and ${\bf C}_3$.

4.3.   Cubic disturbance case ($N=2$)

In this part, we numerically prove that $\mathcal{T}_{\circ}^{\Sigma}$ is the sill of the system (1.4) in the special case of cubic perturbation. So, we firstly introduce this probabilistic model:

associated with the following auxiliary process:

Here, we select $\mathfrak{q}_{12} = 0.014$, $\mathfrak{q}_{22} = 0.0135$, $\mathfrak{q}_{32} = 0.012$ and we keep the other coefficient values as the above two cases. Again, if we select $\mathfrak{b} = 0.4$, then the result is

Theoretically, we have the disappearance of the illness according to Theorem 2.1. It remains to verify it numerically. From Figure 7, the disease will clear up in about 70 days with long-term persistence of categories 1 and 3. Now, we opt $\mathfrak{b} = 0.55$. Then,

From Theorem 3.1, we establish that there is a single stable distribution for (4.3) which is depicted in Figure 8. The persistence of all classes is depicted in Figure 9.

5.   Conclusions

Mathematical formulation plays a major role in understanding epidemics and also in supporting public health decision-making. The regularly used models provide deterministic predictions, that is to say, a strict behavior of the system studied, thus ignoring individual and environmental variations. Actually, we group these unpredictable variations under the name of stochasticity (or randomness), and the present study is devoted to the analysis of an epidemic strewing under heavy stochasticity. The non-linearity and the complexity of the fluctuations pushed us to consider a general form of the probabilistic part. Focusing on these motivations, we have offered an improved generalization of the recent paper ^[9]. In the following, we present our substantial enhancements of the mentioned research.

$\star$ In ^[9], the authors considered a nonlinear prevalence function of the form: $\mathfrak{b}{\bf C}_1(t)g\big({\bf C}_2(t)\big)$. This type of function has certain limitations and is not suitable for covering certain well-known functional responses. By assuming mutual interference between classes ${\bf C}_1$ and ${\bf C}_2$, we have proposed a general function ${\mathcal H}$ which includes all the existing functional incidences.

$\star$ In ^[9], the authors presented two distinct criteria to classify the asymptotic attitude of system (1.4) with the incidence function $\mathfrak{b}{\bf C}_1(t)g\big({\bf C}_2(t)\big)$. Explicitly, they obtained the following results:

$(1)$ The sufficient condition of the disease disappearance is

$(2)$ The sufficient condition of the ergodicity is

However, in this article, we have unified the criterion of the above-mentioned characteristics by providing the following acute threshold value:

Specifically,

$(1)$ The condition of the disease disappearance is $\mathcal{T}_{\circ}^{\Sigma} < 0$.

$(2)$ The condition of the ergodicity is $\mathcal{T}_{\circ}^{\Sigma} > 0$.

As a special case, we mention that the sharp threshold of the perturbed model studied in ^[9] is exactly $\mathcal{R}^C_0$.

Numerically, we have chosen the first three values of $h$, that are, linear, quadratic and cubic perturbations. In all cases, we have confirmed that $\mathcal{T}_{\circ}^{\Sigma}$ is the global sill among disappearance of infection and ergodicity. From Figure 10, we infer that when we raise the perturbation order, the illness disappears swiftly. This indicates that the intense environmental variations have a negative effect on the illness duration. This remark requires further theoretical clarification and explanation. We will deal with this interesting question in our future work.

Some fascinating topics deserve more attention. For example, we can consider our model with fractal-fractional differentiation. This framework is an attractive branch of applied mathematics that deals with derivatives and integrals of non-integer order. Due to its amazing features, it is preferred for describing and simulating real-world problems in various fields such as biological mechanisms, material science, hydrological modeling, economic phenomena. We will address this idea in our future work.

Acknowledgements

This research was sponsored by the Guangzhou Government Project under Grant No. 62216235, also supported by SERB of India (EEQ/2021/001003).

Conflict of interest

The authors declare no conflicts of interest.