Application and analysis of a model with environmental transmission in a periodic environment

1.   Introduction

According to early research, environmental transmission is a crucial driver of infectious diseases, especially enteric and respiratory diseases ^[1]. Some infectious diseases, such as influenza, pertussis and dysentery, show cyclical fluctuations in their transmission dynamics. The cyclical fluctuations of these infectious diseases are mostly strongly associated with seasonal changes and generally consistent with environmental cyclical changes ^[2]. Mathematical models allow the study of disease transmission risks and the dynamic prediction of disease trends. Because the above diseases are highly sensitive to environmental climate fluctuations, autonomous models with constant parameters are no longer applicable, and non-autonomous models are more suitable ^[3,4]. Taking COVID-19 as an example, researchers have proposed many mathematical models and control measures during the first wave of the pandemic to explore the transmission patterns of an early-onset COVID-19 pandemic ^{[5,6,7,8,9,10,11,12,13,14,15]}. Several groups of researchers ^{[16,17,18,19]} have found that the turning point of pandemic fluctuation occurred at the turn of the seasons. Environmental transmission is one of the most critical drivers of the multi-phase outbreak. Other researchers have incorporated viral loads into the environment in mathematical models to study the impact of indirect transmission on pandemics ^[20,21,22]. Musa et al. ^[23] proposed a model of time-varying propagation rates influenced by the environment. However, the periodic oscillations of the environment constitute an essential point that the authors ignore.

Taking into consideration COVID-19, is there a simple and universal mathematical model that can effectively summarize the transmission dynamics of most environmentally transmitted infectious diseases? Several articles ^[24,25,26] show that asymptomatic infected individuals constitute a compartment that cannot be ignored when modeling infectious diseases such as tuberculosis, pertussis and influenza. Others papers ^[2,27,28] show that the transmission rate of some infectious diseases with cyclical fluctuations is closely related to seasonal changes, rather than being static. The study of periodic forcing systems is naturally divided into two main directions: theoretical analysis and application. Researchers focusing on theory have mainly concentrated on system stability and bifurcation analysis. Chithra and Mohamed ^[29] found the existence of strange non-chaotic attractors by using bifurcation and Lyapunov exponents in a single periodic forcing system. M. de Carvalho and Rodrigues ^[30,31] made a breakthrough in the analysis of the bifurcations of periodic forcing systems. They considered the logistic growth in the periodic model and proved the existence of persistent strange attractors. Some researchers ^[32,33] have developed threshold dynamics and defined the next infection operator to analyze the global dynamical behavior of periodic epidemic models. Researchers focusing on model applications have mainly concentrated on the study of the prediction and control of periodic epidemics. Various scholars have set the transmission rates in the model as periodic functions to explore the dynamical behavior of a given disease ^{[28,34,35,36]}. However, few articles have examined a universal model for environmentally transmitted infectious diseases. Based on the work of the above authors, we propose an environmental periodic model to explore the applicability of the model to relevant infectious diseases. The model can inform epidemiological trends in relevant emerging infectious diseases.

The following sections provide model analysis and discussion. Section 2 includes model development, qualitative analysis and stability analysis. Section 3 provides parameter estimation and a comparison of the periodic and non-periodic models for different infectious diseases. Section 4 includes a sensitivity analysis and simulation experiments to explore the impact of changes in the main parameters influencing pandemic transmission dynamics. Section 5 gives an overall discussion of the paper.

2.   Model formulation and qualitative analysis

In contrast to the traditional epidemiological compartmental model, $SEIR$, the periodic environment model is more consistent with the pathogenesis of enteric and respiratory infectious diseases because of the introduction of cyclic environmental changes. Taking COVID-19 as an example, according to a published study recently released by the University of Michigan ^[37], the risk of transmission of the novel coronavirus through aerosols could be thousands of times higher than that of transmission through contact with surfaces. The same is true for the spread of common respiratory infectious diseases such as pertussis, influenza and tuberculosis ^[38]. Enteric infectious diseases like dysentery can spread disease through pathogens in the external environment, such as sewage and feces ^[39]. For this reason, we include the environmental pathogen load, $V(t)$, in the dynamic model.

The model considers two ways of transmission, i.e., direct human-to-human transmission and indirect environment-to-human transmission. Furthermore, five mutually independent epidemiological classes of the host populations are considered: susceptible group $S (t)$, exposed group $E(t)$, asymptomatic infected group $A(t)$, symptomatic infected group $I(t)$ and recovered group $R(t)$ at time $t$. The incubation period is a phase that cannot be ignored for many diseases, and COVID-19 is no exception. In the early stage of infection, some patients do not show symptoms associated with the disease, but they can expel germs or viruses from the body. So, we consider the infectiousness of the exposure period in our model ^[40]. Both asymptomatic and symptomatic patients are confirmed patients. Therefore, we set $E$, $A$ and $I$ as infectious groups. The total population is given as $N(t) = S (t) + E(t) +A(t) + I(t) + R(t)$.

Only the symptomatic compartment $I(t)$ has a disease fatality rate, $d$. $\mu$ represents the natural mortality rate for host compartments. The susceptible populations are increased by recruitment at a constant rate $\Lambda$. Let $\beta_i$, $i = 1, 2, 3$, be the contact transmission rate from $E(t)$, $A(t)$ or $I(t)$ to susceptible $S(t)$. The external environment in which pathogens survive changes periodically and the transmission and decay rates of pathogens are closely related to the fluctuation period of the environment. Therefore, we consider two time-varying functions, $\beta_4(t)$ and $w(t)$, to summarize the cycle variation characteristics, where $\beta_4(t)$ is the function for transmission from the environmental pathogens through cospatial contact or on the surface of objects, and $w(t)$ is the decay function for pathogens. $\beta(t) = \beta_1 \frac{E (t)}{N(t)}+\beta_2 \frac{A (t)}{N(t)}+\beta_3\frac{I (t)}{N(t)}+\beta_4(t) \frac{V (t)}{V(t)+\lambda}$ is the total infectivity of a specific infectious disease through all transmission routes. $\theta$ is the proportion of $E(t)$ converted to $I(t)$, and $1-\theta$ is the proportion of $E(t)$ converted to $A(t)$. The progression rate of class $E(t)$ is denoted by $\kappa$. $c$ is the ratio of $A(t)$ converted to $I(t)$. $\gamma_i$, $i = 1, 2$ represents the recovery rate for $A(t)$ and $I(t)$, respectively. $\eta_i$, $i = 1, 2, 3$ represents the rate of pathogens released into the environment by $E(t), A(t)$ and $I(t)$ through exhalation, talking or some other physical activity such as defecation. Given the actual biological significance of model (2.1), we assume that all parameters in model (2.1) are non-negative. Table 1 details the biological implications of each parameter of the model. Figure 1 gives a more specific flow of the six compartments interacting with each other. The model is given by

where

We introduce the periodic function sin($\omega$t) to describe the effect of periodic changes in the environment on pathogens. $T = \frac{2\pi }{\omega }$ is the period of environmental fluctuations, such as $T = 12$ months, $T = 365$ days, etc. $\beta_4$ is the transmission rate without environmental fluctuations, and $w$ is the decay rate without environmental fluctuations. $a_i \; (i = 1, 2)$ is the magnitude of environmental fluctuations with $0 < a_i < 1$. $a_1 = 0$ indicates that the transmission rate of pathogens is entirely unaffected by environmental changes, and $a_1 = 1$ indicates that changes in the transmission rate are fully consistent with the overall trend of environmental change. These two cases are ideal, so we exclude them. Similarly, $a_2$ has the same interpretation. The initial values for each compartment of model (2.1) are non-negative. The other five compartments, i.e., $S, E, A, I$ and $V,$ are independent of $R$. Model (2.1) can be reduced as follows:

2.1.   Qualitative properties of solutions

Here, we concentrate on the qualitative analysis of system (2.1).

Theorem 1. The solution set $\left \{ S(t), E(t), A(t), I(t), R(t), V(t) \right \}$ of model (2.1) is positive when the initial value for each compartment of the model is non-negative.

Proof. We use the contradiction to prove that $(S(t), E(t), A(t), I(t), R(t), V(t)) \in R_{+}^{6}$ is a solution of system (2.1) for $t \geq 0$. Suppose that $S(t)$ loses positivity for $t_1 > 0$. When $S(t_1) = 0$ holds, we can get $\frac{d S(t_1)}{d t} \leq 0$. However, we know that $\frac{d S(t_1)}{d t} = \Lambda > 0$. It is a contradiction. Therefore, the hypothesis is not valid. The above analysis gives that $S(t) > 0$ for any $t \geq 0$. The second differential equation in model (2.1) can be given by

Solving inequality (2.4) by integration yields

Using the same method, from the third to sixth differential equations of model (2.1), $A(t), I(t), R(t)$ and $V(t)$ are given, respectively, by

Therefore, the above analysis demonstrates that all six state variables remain non-negative at any $t \geq 0$ in model (2.1).

Theorem 2. The solutions of model (2.1) are uniformly bounded in the invariant region $\Omega$, as follows:

Proof. Simplifying model (2.1) yields

Theorem 1 states that each state variable is non-negative at any $t \geq 0$. Therefore, Eq (2.5) is further simplified as follows:

Integrating (2.6) and taking the limit at $t \rightarrow +\infty$ yields

All of the possible solutions of model (2.1) at any $t > 0$ initiating in $R_{+}^{6}$ are confined in the following region: $\Omega = \left\{(S, E, A, I, R, V) \in R_{+}^{6}: 0 < N(t) \leq \frac{\Lambda}{\mu}, 0 \leq V(t) \leq \frac{(\eta_1 + \eta_2 + \eta_3) \Lambda}{\mu w (1-a_2)}\right\}$.

2.2.   Equilibrium points and basic reproduction number

2.2.1.   Disease-free equilibrium point

When E = A = I = V = 0, there is no disease in the system and the disease-free equilibrium, $X^0$, can be obtained as follows:

It shows that the system (2.1) has a unique disease-free equilibrium.

2.2.2.   The basic reproduction number

The basic reproduction number, $R_0$, is an essential epidemiological parameter. It indicates the severity of infectious diseases and can be used as a reference indicator for the degree of prevention and control ^[41]. $R_0$ represents the average number of secondary infections in susceptible individuals through direct contact with a patient or environmental pathogens ^[42].

$R_0$ of system (2.1) is now solved by using the method in ^[43]. The matrix $f$ represents the introduction of the new-infection part, and $v$ represents the generation of the transition part, which are respectively given by

Calculating the Jacobian matrices of (2.7) and then substituting the disease-free equilibrium point yields

Then, $R_0$ of time-averaged autonomous systems is defined as $[R_0] = \rho([F] [V^{-1}])$, where $\rho([F] [V^{-1}])$ is the spectral radius of the next-generation matrix $[F] [V^{-1}]$. $[R_0]$ ^[43] is given by

where

[$R_E$] is the secondary infection caused by an exposed person. In other words, [$R_E$] is the number of susceptible persons who become infected through direct contact with an exposed person. Similarly, [$R_A$] is a secondary infection caused by an asymptomatic infected individual, and [$R_I$] is a secondary infection caused by a symptomatic infected individual. [$R_V$] is a secondary infection caused by pathogens from the environment. In other words, [$R_V$] is the number of susceptible persons infected through indirect contact with pathogens, such as aerosols and waterborne transmission.

However, $[R_0]$ may overestimate or underestimate the risk of disease prevalence. This makes the proposed control measures either too relaxed or too restrictive. Wang and Zhao extend the classical framework described in ^[44] and introduced the next infection operator ^[32,33] in a periodic environment. $R_0$ ^[32] is given by

where $F(s)$ is the newly introduced infected class at time $s$, and $\phi \left (s \right)$ is the distribution of the initialization state of the infected class. Then, $F(s)\phi \left (s \right)$ is a function that descirbes the distribution of new infections generated by the newly introduced infected class at time $s$. $Y(t, s), t \geq s$ is a developmental operator for linear $\omega$-periodic systems $\frac{d y}{d t} = -V(t) y$. The $4 \times 4$ matrix $Y (t, s)$ ^[34], of model (2.1) satisfies

where

with the initial conditions $y_{ii}(0, 0) = 1, i = 1, 2, 3, 4,$ and $y_{ij}(0, 0) = 0, i\ne j.$

$R_0$ of the periodic environment system is defined by the spectral radius of the next-infection operator $L$:

Lemma 1. (see ^[32], Theorem 2.2) $R_{0}$ is strongly correlated with the spectral radius of the next-infection operator, and it satisfies that

(1) $R_{0} = 1$ if and only if $\rho\left(\Phi_{F-V}(\omega)\right) = 1$.

(2) $R_{0} < 1$ if and only if $\rho\left(\Phi_{F-V}(\omega)\right) < 1$.

(3) $R_{0} > 1$ if and only if $\rho\left(\Phi_{F-V}(\omega)\right) > 1$.

If $R_{0}$ satisfies condition $(2)$, it follows that $X^0$ is locally asymptotically stable; if $R_{0}$ satisfies condition $(3)$, it follows that $X^0$ is unstable.

Lemma 2. (see ^[33], Lemma 2.1) Let $r = \frac{1}{\omega} \ln \rho\left(\Phi_{F-V}(\omega)\right)$; then, there exists a positive $\omega$-periodic function $u(t)$ such that $e^{r t}u(t)$ is a solution to $\frac{d z}{d t} = [F(t)-V(t)] z$. (The proof is given in Appendix A.)

Theorem 3. $X^0$ is globally asymptotically stable in the invariant region $\Omega$ when and only when $R_{0} < 1$.

Proof. Lemma 1 states that $X^0$ is locally asymptotically stable when $R_{0} < 1$. Therefore, the next step is to state that $X^0$ is globally attractive. Given that $R_0 < 1$, we have that $\rho\left (\Phi_{F-V}(\omega) \right) < 1$. Now, restrict $\varepsilon$ to a positive number small enough such that $\rho\left(\Phi_{F_\varepsilon-V}(\omega)\right) < 1$, where

System (2.1) has non-negative solutions $(S(t), E(t), A(t), I(t), R(t), V(t))$. Considering the non-disease compartments, it is obvious that there is only the susceptible class, $S$. The differential equation for the susceptible class can be reduced as follows:

Integrating (2.11) yields $S(t)\le \frac{\Lambda }{\mu } +(S(0)-\frac{\Lambda }{\mu })e^{-\mu t}$. It implies that $\lim_{t \to \infty} S(t) = \frac{\Lambda}{\mu } = S^{0}$. Using the standard comparison principle, it follows that, for any $\varepsilon > 0$, there is a $T > 0$ such that $S(t)\le S^{0}+\varepsilon$ for $t > T$. Thus, the following differential inequalities are obtained:

The corresponding auxiliary model can be given by

Lemma 2 states that there exists a $\omega$-periodic function $u(t) = (u_1(t), u_2(t), u_3(t), u_4(t), u_5(t))^T$ such that $(\bar{E}, \bar{A}, \bar{I}, \bar{R}, \bar{V}) = e^{r t} u(t)$ is the solution to model (2.13), where $r = \frac{1}{\omega} \ln \rho\left(\Phi_{F_{\varepsilon}-V}(\omega)\right)$. When $R_{0} < 1$, we have that $\rho\left(\Phi_{F-V}(\omega)\right) < 1$. If $\varepsilon$ is a sufficiently small positive number, we can get that $\rho\left(\Phi_{F_\varepsilon-V}(\omega)\right) < 1$, and then $r < 0$. Therefore, there is a $T > 0$ such that, for any non-negative initial value $x^0$ and sufficiently small $\alpha$, we can get

Following the comparison theorem, for $t > T$, we have

We can get

The above analysis indicates that $\lim_{t \to \infty} S(t) = \frac{\Lambda }{\mu }$, which proves that $X^0$ is globally attractive in $\Omega$ only if $R_0 < 1$. If $R_0 < 1$, the above theorems and analysis show that $X^0$ is globally asymptotically stable.

2.3.   Disease persistence

Definition 1. ^[33] Let $\rho\left(\Phi_{F-V}(\omega)\right) > 1$; $u(t, X_1)$ is a positive $\omega$-periodic solution of model (2.1) if it satisfies the condition that, for any $\varepsilon > 0$, there exists $\delta _0 > 0$ such that, for any $X^1\in X \backslash \left \{ \left (S, 0, 0, 0, 0, 0 \right) \right \}$ with $\left\|X^{1}-X^{0}\right\| \leq \delta_{0}$, we have

Lemma 3. ^[34] The solutions of model (2.1) are uniformly persistent if there exist some positive numbers $\varphi$ such that

The initial values of the six state variables are all positive.

Lemma 4. (See ^[45], Theorem 1.3.1) If $f: X \longrightarrow X$ satisfies the following:

(A) $f(X_0) \subset X_0$ and $f$ has a global attractor $A$;

(B) The maximal compact invariant set $A_\partial = A\cap M_\partial$ of f in $\partial X_0$, possibly empty, admits a Morse decomposition $\left \{ M_1, \cdots, M_k \right \}$ with the following properties:

(a) $M_i$ is isolated in $X$;

(b) $W^s\left \{ M_i \right \} \cap X_0 = \emptyset$ for each $1\le i\le k$,

then $f: X \longrightarrow X$ is uniformly persistent with respect to $(X_0, \partial X_0)$.

Theorem 4. The solutions of model (2.1) are uniformly persistent, and there is at least one positive $\omega$-periodic solution when $R_{0} > 1$.

Proof. Define

Assume that model (2.1) has a unique solution with the initial condition $X^0 = (S^0, E^0, A^0, I^0, R^0, V^0)$, set as $u(t, X^0)$. Let $P: X \longrightarrow X$ be the Poincaré map for model (2.1), that is, $P(X^0) = u(\omega, X^0), \forall X^0 \in X$, where $u(0, X^0) = X^0, P^n(X^0) = u(n\omega, X^0), \forall n > 0$.

First, it is known by definition that $X$ and $X_0$ are positive invariant sets, and that $\partial X_0$ is a relatively closed set in $X$. It has been shown that model (2.1) is always eventually bounded. Therefore, $P$ has a global attractor in $X$.

Next, we focus on proving that $D_{\partial} = \tilde{D_{\partial}}$. It is evident that $\tilde{D_{\partial}} \subseteq D_{\partial}$; then, we only need to prove that $D_{\partial} \subseteq \tilde{D_{\partial}}$. We consider any initial values $(S(0), E(0), A(0), I(0), R(0), V(0)) \in \partial X_{0} \backslash\ \tilde{D_{\partial}}$, which is equivalent to $P^{n}(S(0), E(0), A(0), I(0), R(0), V(0)) \notin \partial X_{0}$. The initial value of the system is discussed in the following 13 cases for three categories (only one compartment's initial value is greater than 0 (1 – 4), two compartments' initial values are greater than 0 (5 – 10) and three compartments' initial values are greater than 0 (11 – 13)):

(1) $V(0) > 0$, and the other three infected compartments $E(0) = 0, I(0) = 0, A(0) = 0$ hold;

(2) $A(0) > 0$, and the other three infected compartments $E(0) = 0, I(0) = 0, V(0) = 0$ hold;

(3) $I(0) > 0$, and the other three infected compartments $E(0) = 0, A(0) = 0, V(0) = 0$ hold;

(4) $E(0) > 0$, and the other three infected compartments $I(0) = 0, A(0) = 0, V(0) = 0$ hold;

(5) There are two infected compartments $E(0) > 0, I(0) > 0$, and the other two infected compartments $A(0) = 0, V(0) = 0$ hold;

(6) There are two infected compartments $E(0) > 0, A(0) > 0$, and the other two infected compartments $I(0) = 0, V(0) = 0$ hold;

(7) There are two infected compartments $E(0) > 0, V(0) > 0$, and the other two infected compartments $I(0) = 0, A(0) = 0$ hold;

(8) There are two infected compartments $I(0) > 0, A(0) > 0$, and the other two infected compartments $E(0) = 0, V(0) = 0$ hold;

(9) There are two infected compartments $I(0) > 0, V(0) > 0$, and the other two infected compartments $E(0) = 0, A(0) = 0$ hold;

(10) There are two infected compartments $A(0) > 0, V(0) > 0$, and the other two infected compartments $E(0) = 0, I(0) = 0$ hold;

(11) There are three infected compartments $E(0) > 0, I(0) > 0, A(0) > 0$, leaving one compartment $V(0) = 0$;

(12) There are three infected compartments $E(0) > 0, I(0) > 0, V(0) > 0$, leaving one infected compartment $A(0) = 0$;

(13) There are three infected compartments $I(0) > 0, A(0) > 0, V(0) > 0$, leaving one infected compartment $E(0) = 0$;

When case (1) holds, we can get

Integrating (2.14) yields that $E(t) > 0$ for any $0 < t \ll 1$. Solving for the remaining four state variables at $0 < t \ll 1$ yields

The above results imply that $\left (S(t), E(t), A(t), I(t), R(t), V(t) \right) \notin \partial X_0$ for any $0 < t \ll 1$. The proofs of cases (5) and (11) in the other two categories are given in Appendix B, respectively, and the proofs of the other cases in the category can be obtained similarly. It shows that

In other words, if there exists $(S(t), E(t), A(t), I(t), R(t), V(t))$, which is a solution of the system from $D_{\partial }$, then the limit of the solution must satisfy that

Therefore, $X^{0}$ is an isolated invariant set. The above analysis proves that $D_{\partial} = \tilde{D_{\partial}}$.

Finally, we prove that $W^{s}\left(X^{0}\right) \cap X_{0} = \emptyset$, where $W^{s}\left(X^{0}\right) = \left\{X^{0}: P^{n}\left(X^{0}\right) \rightarrow X^{0}, n \rightarrow +\infty\right\}$. The continuous dependence of the solution of the differential equation on the initial value shows that, for any $\varepsilon > 0$, there exists $\delta _0 > 0$; for any $X^1\in X_0$, with $\left\|X^{1}-X^{0}\right\| \leq \delta_{0}$, we have

We claim that

Proving (2.15) by contradiction yields

Let us assume that

Thus,

Furthermore, for any $t > 0$, we have that $t = n \omega+t^{\prime}, t^{\prime} \in[0, \omega], n = \left[\frac{t}{\omega}\right]$.

Then, the following equation can be obtained:

Let $(S(t), E(t), A(t), I(t), R(t), V(t)) = u(t, X^1)$, where

Then, we have

Now, restrict $\varepsilon$ of model (2.16) to be sufficiently small and consider the following auxiliary system:

According to Lemma 2, there exists a $\omega$-periodic function $p(t) = (p_1(t), p_2(t), p_3(t), p_4(t), p_5(t))$ such that $(\bar{E}, \bar{A}, \bar{I}, \bar{R}, \bar{V}) = e^{r_{1} t} p(t)$ is the solution to system (2.17), where $r_{1} = \frac{1}{\omega} \ln \rho\left(\Phi_{F_{\varepsilon}-V}(\omega)\right.).$ According to Lemma 1, we have that $\rho\left(\Phi_{F-V}(\omega)\right.) > 1$ when $R_{0} > 1$. Now, restricting $\varepsilon$ to be a sufficiently small positive number yileds that $\rho\left(\Phi_{F_{\varepsilon}-V}(\omega)\right.) > 1$; then, $r_1 > 0.$ Hence, for $T_1 > 0$ and $\alpha$ small enough, we have

Following the comparison theorem yields

It shows that $\lim_{t \rightarrow \infty} E(t) = +\infty, \; \lim_{t \rightarrow \infty} A(t) = +\infty, \; \lim_{t \rightarrow \infty} I(t) = +\infty, \; \lim_{t \rightarrow \infty} R(t) = +\infty$ and $\lim_{t \rightarrow \infty} V(t) = +\infty$. It contradicts what was mentioned earlier, which proves that $W^{s}\left(X^{0}\right) \cap X_{0} = \emptyset$. Hence, $X^0$ is acyclic in $D_\partial$ and $P$ is uniformly persistent with respect to $(X_0, \partial X_0)$. It shows that the solutions of model (2.1) are uniformly persistent. The Poincaré map $P$ has a stationary point $\tilde{X}(0) = (\tilde{S}(0), \tilde{E}(0), \tilde{A}(0), \tilde{I}(0), \tilde{R}(0), \tilde{V}(0)) \in X_{0}$ with $\tilde{S}(0) > 0$. Thus, $\tilde{X}(0)\in Int (R_+^6)$ and $\tilde{X}(t) = u(t, \tilde{X}(0))$ is a positive $\omega$-periodic solution of the model.

To visualize the uniform persistence of the solutions when $R_0 > 1$, Figure 2 is presented to simulate the spread of the pandemic through the use of numerical examples. Taking the values of the parameters in Table 2 and setting $w = \frac{\pi}{180}, a_1 = a_2 = 0.8$, it is observed that all six compartments eventually oscillate periodically with time and persist within a specific range. Two three-dimensional phase diagrams of the virus-containing compartments are given in Figure 3. The solutions of model (2.1) eventually stabilize in a range of cycles with different initial conditions.

3.   Model applicability

We selected several typical infectious diseases for data analysis to illustrate the applicability of the proposed periodic model for respiratory or intestinal infectious diseases. For COVID-19, we calibrated model (2.1) by fitting actual case data for Shanghai from March 1, 2022 to August 31, 2022 ^[46]. Daily new and cumulative cases were taken, and the cumulative asymptomatic cases and symptomatic cases are denoted as follows:

The least squares method was used to fit model (2.1). We fit all parameters and $S(0)$ of model (2.1), except $\mu$ and $\lambda$. $\lambda$ and $\mu$ were obtained from published literature ^[23] and a local government website ^[46], respectively. $E (0)$, $A (0)$, $I (0)$ and $R (0)$ were taken from the number of reported cases from Shanghai's outbreak statistics on March 1, 2022. Thus, we set $E (0) = 80$, $A (0) = 1$, $I (0) = 1$, $R (0) = 0$ and assumed that $V (0) = 0$. The climate is different throughout the year, and there is also a certain gap in the survivability of the virus in the environment; hence, we set the oscillation period to one year ($\omega = \frac{\pi }{180}$). After continuous program debugging of the data and model (2.1), we selected a set of optimal parameters for the environmental oscillation magnitude, that is, $a_1 = 0.009$, $a_2 = 0.01$. Table 2 shows the exact fitting values of model (2.1) parameters.

In Figure 4, we fit model (2.1) to confirmed case data in Shanghai. The red line represents the fitted curve for the infected class in model (2.1). The blue points represent reported cases. Figure 4(a), (b) show that the fitted curve of $A(t)$ broadly matches the actual epidemiological trend of the outbreak. But, it may ultimately underestimate the cumulative cases. Figure 4(c), (d) show a certain gap in the fitting of daily new cases for $I(t)$, even though the fitting trend of cumulative cases is consistent. To further validate the applicability of model (2.1) to COVID-19, the residuals in statistics are used to assess the merits that

where $D(t)$ is the daily predicted value for the asymptomatic or symptomatic infected populations, and $d(t)$ is the daily observed value for the asymptomatic or symptomatic infected populations. As shown in Figure 5, the residuals are small and randomly distributed. Therefore, we can conclude that the proposed periodic model is applicable to COVID-19.

Next, to further illustrate the advantages of the periodic environment model compared to the autonomous model, we propose a non-periodic model (3.1) based on model (2.1):

Model (3.1) is fitted to actual epidemiological data and compared with model (2.1). Figure 6 shows that model (2.1) fits the actual epidemic data much better, while model (3.1) underestimates the epidemic peak. The sizes of the viral carriage compartments, i.e., $(E(t), A(t), I(t))$, were all underestimated based on the overestimation of $S(t)$. In contrast, model (3.1) predicts a much higher viral load on the environment than model (2.1). The similarity between the two systems is that there are some small fluctuations after the end of the first wave of the pandemic, suggesting that the trend in the prevalence of COVID-19 is affected by environmental fluctuations. For this reason, it is timely and reasonable to consider the periodic environment model.

To observe the epidemic trend over a long period, Figure 7 was constructed to show the time-response diagram for system (2.1) over 8 years. Due to the periodic environmental changes, the epidemic shows periodic fluctuations, with a decreasing trend in wave peaks and an increasing trend in wavelengths. As the epidemic continues to develop, the natural immunity of humans will increase to a certain extent and eventually develop to the point where the virus coexists with humans, as the message conveys in Figure 7.

To further validate the efficiency of model (2.1) in capturing the cyclical changes in the environment, we selected seasonal highly prevalent infectious diseases (influenza, dysentery and pertussis) and perennial diseases (tuberculosis) for our experiments. Since some infectious diseases, such as tuberculosis, are not contagious during incubation ^[47], we set $\beta_{1} = 0$. For influenza, we selected China's case data from August 2018 to October 2021 ^[40] and fitted it with model (2.1) and model (3.1). Figure 8(a) shows that model (2.1) can better capture the dynamics of influenza over time, and that the model forecasts high disease periods that are entirely consistent with reality. Due to the variability of influenza viruses, there may be discrepancies between the predicted and actual case data. The performance of model (3.1) in terms of epidemiological trends and accuracy is quite disappointing. For dysentery, we used cumulative infections in China from October 2017 to April 2020 for our analysis ^[40]. As can be seen in Figure 8(b), model (2.1) fully captures the seasonal high-prevalence characteristics of dysentery, and the data fitting error is small enough to be negligible. As expected, model (3.1) failed to reflect the seasonal trend of dysentery incidence. For pertussis, two years of case data from October 2017 to September 2019 in China were selected ^[40]. Figure 8(c) shows that the periodic and non-periodic models predict the same trend, but the periodic model fits the data relatively better. For perennial tuberculosis, 24 months of reported cases from January 2021 to December 2022 were selected for China ^[40]. Figure 8(d) shows that, for infectious diseases with no seasonal characteristics, the sensitivity of the periodic and non-periodic models to the data is basically the same. The periodicity advantage described by model (3.1) does not apply. In addition, noise may be more favorable if the effects of environmental oscillations on perennial infectious diseases, or the stochastic interference of certain uncertainties in contagious diseases, are to be studied ^[48,49].

4.   Sensitivity analysis and numerical simulation

This section explores the parameters with critical drivers in model (2.1). Sensitivity analysis has been used to obtain the parameters that significantly influence $R_0$ of model (2.1). Therefore, some intervention strategies can be proposed to help reduce the impact of infection peaks and environmental fluctuations. On the other hand, local officials can use the sensitivity analysis results to develop effective interventions to reduce the impact of infectious diseases for which environmental transmission is the main driver. Now, the normalized forward sensitivity index of $R_0$ for the parameters is defined by

where $p$ represents the various parameters in $R_0$. If the sensitivity index of a parameter is positive, it indicates that the parameter has to promote an effect on $R_0$. Similarly, parameters with a negative sensitivity index have an inhibitory effect on $R_0$.

Using COVID-19 as the subject of the study, for computational convenience, we use the fitted values in Table 2 to examine the effect of each parameter on $R_0$; we also give the sensitivity index values for each parameter. Figure 9 clearly shows that the rate of transmission via indirect contact with the virus $(\beta_4)$ has a much more significant impact on $R_0$ than does direct contact on $R_0$, and that the viral decay rate $w$ has a great inhibitory effect on $R_0$. Hence, we need to be highly concerned about the viral load on the environment. This result is the same as in the literature ^[20,21,22], which once again demonstrates the applicability of periodic environment model (2.1) to COVID-19. Based on the above analysis, we can determine that $\beta_{4}$ and $w$ significantly impact $R_0$. For COVID-19, the focus should be on the virus in the environment.

Therefore, numerical simulation describes the influence of these critical parameters on the system. First, we focus on the impact of changes in $\beta_{4}$ on the system. As shown in Figure 10, the compartments respond sharply during the initial phase of infection, and the infected population increases with increasing $\beta_{4}$. The increase of $\beta_{4}$ shortens the outbreak period and causes it to move into the coexistence phase earlier. This may be due to herd immunity. The red dashed line, shown in Figure 10, represents the lowest level of $\beta_{4}$. The second wave is the latest and has the longest wavelength. If we reduce the values of $\beta_1$, $\beta_2$, $\beta_3$ and $\beta_4$ by the same degree, as shown in Figure 11, the system that changes $\beta_{4}$ is the first to respond. Reducing $\beta_{4}$ will significantly reduce the number of infected individuals relative to the original system. It will push back the second wave oscillation, reducing the size of the second wave. In particular, $\beta_{4}$ represents the virus's transmission rate. We can reduce the number of infections via indirect contact in the post-epidemic era by washing our hands and ventilating regularly. Finally, the influence of the decay rate on the system should not be ignored. As shown in Figure 12, although decreasing the value of $w$ will increase the size of the first wave of infection and make the second wave arrive earlier, it will shorten the time to reach stability and reach a state of coexistence as soon as possible. $w$ represents the decay rate of the virus, and we can reduce the impact of infection by implementing environmental decontamination.

5.   Conclusions

Periodic oscillations of the environment are often considered in dynamic models of respiratory and enteric infectious diseases. Mitigating the impact of COVID-19 is an expectation of people worldwide. Inspired by the mechanism of transmission of COVID-19, a new environmental periodic model, SEAIRV, has been proposed based on the previous SEIR model, considering the infectivity of exposed and asymptomatic patients and the influence of environmental cycle oscillation on pathogens. It was our goal to explore the applicability of the periodic model. First, the disease-free equilibrium point and $R_0$ of the system were calculated. The global asymptotic stability of the disease-free equilibrium point and the uniform persistence of the disease have been demonstrated. Second, the applicability of model (2.1) to related infectious diseases was explored. We used the least squares method to fit model (2.1) with COVID-19-related reported data from Shanghai. The results show that model (2.1) better reflects the epidemiological trends in Shanghai than model (3.1). Using residuals and model comparisons powerfully demonstrates that the periodic environment model generalizes the transmission characteristics of COVID-19 better than the constant-parameter model. Meanwhile, we validated the applicability of model (2.1) by using case data for four other infectious diseases with similar pathogeneses to COVID-19. The results show that, for seasonal high-incidence infectious diseases, such as influenza, dysentery and pertussis, model (2.1) can well reflect the seasonal characteristics of infectious diseases. In contrast, the deterministic model (3.1) cannot do this. For perennial infectious diseases, such as tuberculosis, model (2.1) is about as effective as model (3.1). Finally, a sensitivity analysis of the system was carried out by using the fitted parameter values for COVID-19. The results show that the parameters related to pathogens in the environment greatly influence the system. Increasing the transmission rate of the virus or the decay rate of the virus can shorten the time for the system to stabilize and causing the system to reach equilibrium as quickly as possible. It reveals that, in the post-pandemic era, we need to focus on viral loads on the environment. With humankind's combined efforts, we will finally overcome this epidemic.

Use of AI tools declaration

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

Conflict of interest

All authors declare no conflict of interest that may affect the publication of this paper.

Appendix A. Proof of Lemma 2

Proof. It follows from Lemma 2.1 of ^[33] that $F(t)-V(t)$ is a continuous, cooperative, irreducible and $\omega$-periodic $k\times k$ matrix function. $\Phi _{\left (F\left (\cdot \right) -V\left (\cdot \right) \right) }\left (t \right)$ is the fundamental solution matrix of the linear system $\frac{d z}{d t} = [F(t)-V(t)] z$, and $\rho\left(\Phi_{F\left (\cdot \right)-V\left (\cdot \right)}(\omega)\right)$ is the spectral radius of $\Phi _{\left (F\left (\cdot \right) -V\left (\cdot \right) \right) }\left (\omega \right)$. Assume that $u^* \gg 0$ is the eigenvector corresponding to the principal eigenvalue $\rho\left(\Phi_{F\left (\cdot \right)-V\left (\cdot \right)}(\omega)\right)$. It is known that $z(t) = e^{r t}u(t)$, the derivation of which yields

Define $u \left (t \right) : = \Phi _{\left (F\left (\cdot \right) -V\left (\cdot \right) -rI \right) }\left (t \right)u^*$ to be a positive solution of the above system. It can be obtained that $e^{rt}\Phi _{\left (F\left (\cdot \right) -V\left (\cdot \right) -rI \right) }\left (t \right) = \Phi _{\left (F\left (\cdot \right) -V\left (\cdot \right) \right) }\left (t \right)$. In addition,

Hence, $u(t)$ is a positive $\omega$-periodic solution of the system, which proves Lemma 2.

Appendix B. Proof of cases (5) and (11)

When case (5) holds, we can get

Integrating the above equations gives $A(t) > 0, V(t) > 0$ for any $0 < t \ll 1$. Solving for the remaining four state variables at $0 < t \ll 1$ gives

The above results imply that, when case (5) holds, $\left (S(t), E(t), A(t), I(t), R(t), V(t) \right) \notin \partial X_0$ for any $0 < t \ll 1$.

When case (11) holds, we can get

Integrating the above equation gives $V(t) > 0$ for any $0 < t \ll 1$. Solving for the remaining four state variables at $0 < t \ll 1$ gives

The above results imply that, when case (11) holds, $\left (S(t), E(t), A(t), I(t), R(t), V(t) \right) \notin \partial X_0$ for any $0 < t \ll 1$.