The impact of ambient air pollution on an influenza model with partial immunity and vaccination

1.   Introduction

Influenza is an acute respiratory disease caused by influenza viruses, types A, B, C and D, which circulate in all parts of the world. The World Health Organization (WHO) reported that influenza occurs globally with an annual attack rate estimated at $5\%-10\%$ in adults and $20\%-30\%$ in children, and these annual epidemics are estimated to result in about 3,000,000 to 5,000,000 cases of severe illness and about 290,000 to 650,000 deaths ^[1].

In recent years, ambient air pollution (AAP), contamination of the indoor/outdoor environment by any chemical, physical or biological agent that modifies the natural characteristics of the atmosphere, is one of the greatest environmental risks to health and has become a severe problem across the world. Environmental toxins spread at all levels of biological systems, extending from atoms to biospheres, and they have impacts on cells, organs, organisms, populations and the whole ecosystem. For example, fog and haze can be taken as a proof of air pollution resulting in health problems in big cities, containing harmful toxic substances like sulfur dioxide, nitrogen oxides, etc., invading bodies and causing health damage. WHO estimated in 2019 that ambient (outdoor) air pollution caused 4.2 million premature deaths worldwide, and this mortality is due to exposure to fine particulate matter (PM), which causes cardiovascular and respiratory disease, and cancers ^[2]. It is reported that particle droplets are considered to be a primary transmission mode for influenza virus infection due to sneezing and coughing ^[3]. Therefore, contagious PM may cause influenza disease spread ^[4,5].

Mathematical modeling approaches in epidemiology have been a great tool to reveal deep understanding of the mechanisms underlying the spread and thus provide effective methods for the control of influenza disease ^{[6,7,8,9,10,11,12]}. Recently, some epidemiological studies revealed associations of air pollution with influenza infection. There are reports that PM$_{2.5}$, PM$_{10}$, SO$_2$, NO$_2$ and CO are associated with the influenza positivity rate ^[13,14,15]. Huang et al. ^[16] investigated the the adverse health effects of air pollution and the cause of influenza-like illness (ILI). Meng et al. ^[17] observed that air pollution may be associated with the risk of influenza in a broad sense. Lu et al. ^[18] reported a greater short-term impact on childhood pneumonia from PM$_1$ in comparison to PM$_{2.5}$ and PM$_{10}$. Taken together, air pollution could be one of the determinant factors of influenza infection. It should be noted that the results above are based on the statistical analysis. However, it is rare to study the effect of AAP on influenza transmission by a using dynamical model.

Huaian is located in the ranges of $32^{\circ}43'00''-34^{\circ}06'00''$ N and $118^{\circ}12'00''-119^{\circ}36'30''$ E, in the transition region between the southern warm temperate zone and the northern subtropical zone. Huaian, the central area of northern Jiangsu province, China, is composed of four districts, namely, Qingjiangpu, Huaiyin, Huaian and Hongze, and three counties, namely, Lianshui, Xuyi and Jinhu (see Figure 1). The maximum linear distance between east and west is 132 kilometers, and the maximum linear distance between north and south is 150 kilometers, covering an area of 10.03 thousand square kilometers. By the end of 2021, the permanent population of Huaian was 4.5622 million.

Beyond all doubt, influenza is an acute epidemic in Huaian, China. The data of reported cases of ILI and the AAP, including PM$_{2.5}$, PM$_{10}$, SO$_2$, NO$_2$, O$_3$ and CO, from the 1$^{{{\mathrm{st}}}}$ week of 2019 to 52$^{{{\mathrm{nd}}}}$ week of 2019 in Huaian, China, are shown in Figure 2.

Naturally, there comes a question: What is the relation between the incidence of ILI and the AAP?

The main goal of this paper is to investigate the effect of the AAP on the incidence of ILI in Huaian, China, and to provide some useful strategies to control the spreading of influenza. The rest of the paper is organized as follows: In the following section, we formulate an environmental AAP-dependent SIPS influenza model with partial immunity and vaccination. In Section 3, we study the dynamics of the proposed model, including the positivity and boundedness of the model, equilibrium analysis with basic reproduction number. In Section 4, based on the statistic data, we estimate the parameters of the model, and study controlling the prevalence of influenza in Huaian, China. The last section includes a general discussion and some concluding remarks.

2.   Model formulation

Suppose that, in the absence of air pollution, the total population $N(t)$ is divided into three parts, the susceptible $S(t)$, the infectious $I(t)$ and the partial immunity $P(t)$, i.e., $N(t) = S(t)+I(t)+P(t)$. Recovered individuals from influenza infection are not permanently immune to the next influenza virus ^[19]. The influenza model incorporating partial immunity of the recovered individuals and vaccination of the susceptible is presented by the following SIPS model:

The biological meanings of the parameters of model (1) are shown in Table 1.

To mechanistically incorporate the AAP into the basic SIPS model (1), we apply some of the existing models in this direction ^{[20,21,22,23]}. Let $E(t)$ be the concentration of the AAP in the environment, and $C(t)$ is the concentration of the AAP in the population, at time $t\geq0$, respectively. Assume that the AAP may be externally introduced into the environment according to some prescribed rate $a(t)$. The AAP in the environment is washed out or broken down with time-dependent rate $\eta(t)$, which may occur if the environment is a lake which on occasion drains into another body of water, or if the AAP is subject to chemical decomposition. Further, we assume that the AAP from the environment is absorbed by the population in direct proportion to its concentration (i.e., $\rho N E$), and the AAP in the population may also be removed from the total environment directly with rate $\xi$. We apply the term $mC(t)$ to describe the impact of the AAP on the transmission rate, and then the transmission rate becomes

The schematic diagram is provided in Figure 3 for clear understanding of the model system. Then, we can obtain the extended model as follows:

The biological meanings of all parameters of model (3) are also shown in Table 1.

It is true that the amount of the AAP consumed by the population is usually very small in quantity, and hence the term $-\rho N E$ in the fifth equation in model (3) can be neglected. Then, the research object of this paper is as follows:

The total population $N(t)$ satisfies the following system:

Then,

The initial values of model (4) are:

3.   Dynamics analysis

3.1.   Positivity and boundedness

In this subsection, we analyze the positivity and boundedness of model (4) to ensure that the model is well-posed.

Theorem 1. All solutions of model (4) that start in $\mathbb{R}_+^5 = \{(S, I, P, C, E): S > 0, I > 0, P > 0, C\geq0, \, E > 0\}$ remain non-negative for all the time.

Proof. Since the right hand side of model (4) is continuous and locally Lipschitzian on $\mathbb{C}$ (space of continuous functions), the solution $(S(t), I(t), P(t))$ of model (4) exists and is unique on $[0, \tau)$, where $0 < \tau \leq +\infty$ ^[24]. From the second equation of model (4), we can obtain that

which implies that when $I(0) > 0$, $I(t) > 0$ for $t > 0$. Note that $I(t)$ may tend to zero when $t\rightarrow \infty$. Thus, we put $I(t)\geq0$. Next, we show that, $S(t) > 0$, $\forall t\in [0, \tau)$. If it does not hold, then $\exists\, t_1\in [0, \tau)$ such that $S(t_1) = 0$, $\dfrac{{{\mathrm{d}}} S}{{{\mathrm{d}}}t}\bigg|_{t = t_1}\leq0$ and $S(t) > 0, \forall\, t\in[0, t_1)$. So, there must be $P(t) > 0, \forall\, t\in[0, t_1).$ Suppose the statement is not true. Then, $\exists\, t_2\in [0, t_1)$ such that $P(t_2) = 0$, $\dfrac{{{\mathrm{d}}} P}{{{\mathrm{d}}}t}\bigg|_{t = t_2}\leq0$ and $P(t) > 0, \forall\, t\in[0, t_2)$. From the third equation of (4), we get

which is a contradiction to $\dfrac{{{\mathrm{d}}} P}{{{\mathrm{d}}}t}\bigg|_{t = t_2}\leq0$. So, $V > 0, \forall\, t\in[0, t_1)$. Following from the first equation of (4), we have

which is a contradiction to $\dfrac{{{\mathrm{d}}}S}{{{\mathrm{d}}}t}\bigg|_{t = t_1}\leq0$. It shows that $S(t) > 0$, $\forall t\in [0, \tau)$.

Next, we claim $C(t)\geq0$, $\forall t\in [0, \tau)$. Otherwise, $\exists\, t_3\in [0, \tau)$ such that $C(t_3) = 0$, $\dfrac{{{\mathrm{d}}} C}{{{\mathrm{d}}}t}\bigg|_{t = t_3} < 0$ and $C(t)\geq0, \forall\, t\in[0, t_3)$. Hence, there must be $E(t)\geq0, \forall\, t\in[0, t_3).$ If this statement is not true, then $\exists\, t_4\in [0, t_3)$ such that $E(t_3) = 0$, $\dfrac{{{\mathrm{d}}} E}{{{\mathrm{d}}}t}\bigg|_{t = t_3} < 0$ and $E(t)\geq0, \forall\, t\in[0, t_4)$. From the fifth equation of (4), we get

which is a contradiction to $\dfrac{{{\mathrm{d}}} E}{{{\mathrm{d}}}t}\bigg|_{t = t_2} < 0$. So, $E\geq0, \forall\, t\in[0, t_3)$. Following from the fourth equation of (4), we have

which is a contradiction to $\dfrac{{{\mathrm{d}}}C}{{{\mathrm{d}}}t}\bigg|_{t = t_3} < 0$. It shows that $C(t)\geq0$, $\forall t\in [0, \tau)$. Therefore,

This completes the proof.

It is easy to know that, for model

there is a unique continuous periodic solution $E^*(t)$ which is globally asymptotically stable, where

The study of the dynamics of model (4) requires the introduction of the following set:

where $\hat{E} = \sup\limits_{t\in [0, \omega)}E^*(t), \, \hat{C} = \dfrac{\rho \Lambda \hat{E}}{\mu\xi}$. $\Gamma$ is a positive invariant set for model (4).

3.2.   The basic reproduction number and threshold dynamics

When the disease dies out, we can get the following system from model (4):

Then, model (7) with initial condition $(S(0), P(0), C(0))$ has a unique positive $\omega$-periodic solution $u^*(t) = (S_*, P_*, C_*(t))$, where

Let $X(t) = (I(t), C(t), S(t), P(t), E(t))^{\textbf{T}}$ and then model (4) admits a unique disease-free $\omega$-periodic solution

where $E^*(t)$ is defined as (6), and $C_*(t), S_*(t), P_*(t)$ satisfy (7).

Remark 1. If $a(t) = 0$, i.e., the exogenous input rate of the AAP into the environment is zero, it follows that

which implies that the AAP is eventually washed out of the environment in its totality. This agrees with our intuition. Of course, the time for this occurrence could be very long if the washout rate is small.

The basic reproduction number $\mathcal{R}_0$ is the number of new infective individuals that are generated by a single infectious individual in a completely susceptible population. Now, we define the basic reproduction number $\mathcal{R}_0$. Thanks to ^[25,26], let $\mathcal{F}(t, X)$ be the input rate of newly infected individuals and $\mathcal{{V}}(t, X)$ be the rate of transfer of individuals. Then,

Then,

We can see that the eigenvalues of $-V$ are the diagonal elements and negative.

Let $Y(t, s), t\geq s$ be the evolution operator of the linear $\omega-$periodic system

That is, for each $s\in \mathbb{R},$ the $2\times 2$ matrix $Y(t, x)$ satisfies

where ${\textbf{I}} = {{\mathrm{diag}}} (1, 1)$ is the identity matrix. Following ^[27], let $\phi(s)$ be $\omega$-periodic in $s$ and the initial distribution of infectious individuals, so $F(s)\phi(s)$ is the rate of new infections produced by the infected individuals who are introduced at time $s$. When $t\geq s$, $Y(t, s)F(s)\phi(s)$ gives the distribution of those infected individuals who are newly infected by $\phi(s)$ and remain in the infected compartments at time $t$. Naturally,

is the distribution of accumulative new infections at time $t$ produced by all those infected individuals $\phi(s)$ introduced at time previous to $t$.

Let $C_\omega$ be the ordered Banach space of all $\omega$-periodic functions from $\mathbb{R}$ to $\mathbb{R}^2$, which is equipped with the maximum norm $\|\cdot\|$ and the positive cone $C_\omega^+: = \{\phi\in C_\omega:\phi(t)\geq0, \forall t\in\mathbb{R} \}$. Then, we define a linear operator $L$ which implies that

which is called the next infection operator, and the spectral radius of $L$ is defined as the basic reproduction number

for the periodic epidemic model (4).

Furthermore, we can obtain the global threshold dynamics of model (4).

Theorem 2. If $\mathcal{R}_0 < 1$, model (4) admits a unique disease-free periodic solution $X_*(t) = (0, C_*(t), S_*, V_*(t), E^*(t))$ which is globally asymptotically stable, while if $\mathcal{R}_0 > 1$, it is unstable. If $\mathcal{R}_0 > 1$, there exists a positive constant $\zeta > 0$ such that for any initial values $X(0)\in \Gamma$, the solution of model (4) satisfies

That is, model (4) is uniformly persistent.

The proof of Theorem 2 is standard, and we omit it here. We refer the reader to the proof of Theorem 3.9 in ^[28].

3.3.   The special autonomous case

We now focus on the influenza transmission dynamics of model (4) in a special case of $\beta(t) = \beta, \, a(t) = a, \, \eta(t) = \eta$. In this case, model (4) can be rewritten as follows:

Model (10) always has a disease-free equilibrium (DFE) $X_* = \left(0, C_*, S_*, V_*, E_*\right)$, where

The basic reproduction number $\mathcal{\overline{R}}_0$ of model (10) is

If ${\mathcal{\overline{R}}}_0 > 1$, model (10) has an endemic equilibrium (EE) $X^* = \left(I^*, C^*, S^*, V^*, E^*\right)$, where

Three eigenvalues of the Jacobian matrix of model (10) at the endemic equilibrium $X^*$ are $-\mu, -\xi, -\eta$, and the other two satisfy the following quadratic equation with respect to $\zeta$:

Obviously, the solutions of (12), $\zeta_{1, 2}$, are less than 0, and hence the endemic equilibrium $X^*$ of model (10) is locally asymptotically stable.

Theorem 3. If $\mathcal{R}_0 > 1$, the endemic equilibrium $X^*$ of model (10) is global asymptotically stable.

Proof. In order to study the globally asymptotical stability of $X^*$, we consider

Then, $N, C$ and $E$ approach $N^*, C^*$ and $E^*$ as $t\rightarrow \infty$, respectively. Asymptotically the second and third equations of model (10) can be rewritten as

It follows that model (13) has an interior equilibrium $\overline{X} = (I^*, P^*)$. Two eigenvalues of the Jacobian matrix of model (13) at $\overline{X}$ satisfy (12). Hence $\overline{X}$ is locally asymptotically stable.

We now prove that model (13) has no periodic orbits. We choose the Dulac function

Let

Then,

According to the Bendixson-Dulac criterion, model (13) has no periodic orbits, so the interior equilibrium $\overline{X} = (I^*, P^*)$ is globally asymptotically stable, which completes the proof.

4.   Numerical results via influenza disease dynamics

4.1.   Data collection

In China, influenza has been classified as a class Ⅲ statutory reportable disease in the National Infectious Disease Reporting System. Hospitals and clinics collect nasopharyngeal swabs for each confirmed case, then send them to a designated laboratory for virus isolation and further identification and submits the result online within 24 h ^[17,29]. The data of ILI cases in Huaian city during the 1$^{{{\mathrm{st}}}}$ week of 2019 (the data of 2019 is used for model simulation) to the 12$^{{{\mathrm{th}}}}$ week of 2020 (the data of 2020 is used to test the model prediction) were derived from the Huaian Center for Disease Control and Prevention, which included information about ILI case number, vaccination number and the demographic data ^[30]. The ILI case is defined as a person with a sudden onset of fever $(\geq 38^{\circ} C)$, chills, cough and/or sore throat, a generalized feeling of weakness and pain in the muscles, together with varying degrees of soreness in the head and abdomen. The information regarding the AAP was gathered from the China Online Air Quality Monitoring and Analysis Platform ^[31]. The details can be found in Figure 2.

4.2.   The selection of $E(t)$ in the model

In model (4), $E(t)$ is the concentration of the AAP in the environment. For studying the relations between ILI and the AAP (e.g., PM$_{2.5}$, PM$_{10}$, SO$_2$, NO$_2$, O$_3$ and CO), we compute the Spearman's rank correlation coefficients between them.

Spearman's correlation coefficient, the strength of the rank correlation between variables, is given by

where $\sum_{i = 1}^n {d_i^2}$ represents the sum of the squared differences between $x = (x_1, x_2, \cdots, x_n)$ and $y = (y_1, y_2, \cdots, y_n)$ variable ranks, and $n$ is the sample size.

Form Table 2, we can see that the Spearman's rank correlation coefficients between PM$_{2.5}$, PM$_{10}$, SO$_2$, NO$_2$, O$_3$, CO and ILI are 0.34, 0.40, 0.07, 0.02, -0.17, 0.56, respectively, which indicates that PM$_{2.5}$, PM$_{10}$, SO$_2$, NO$_2$ and CO are positively significantly correlated with the incidence of ILI, while O$_3$ is negatively significantly correlated with the incidence of ILI. Noting that $\rho_s$ of CO, PM$_{10}$ and PM$_{2.5}$ are the three largest Spearman's rank correlation coefficients, we choose CO, PM$_{10}$ and PM$_{2.5}$ as $E(t)$ and focus on the effects of these factors on the prevalence of influenza in Huaian, China, respectively.

4.3.   The values of parameters

From the statistical data ^[30,31,32], we adopt

and the initial values are taken as $S(0) = 49,885$, $I(0 = 623$, $V(0) = 0$, $C(0) = 100$.

To parameterize the model, according to the statistical data of CO, PM$_{10}$ and PM$_{2.5}$ in Huaian, China ^[31], we use the Matlab function fminsearch with Nelder-Mead algorithm and estimate $a(t)$ and $\eta(t)$ via model (5).

1) With $E_{{{\mathrm{CO}}}}(0) = 1.04$, we can obtain

2) With $E_{{{\mathrm{PM}}}_{10}}(0) = 107$, we can obtain

3) With $E_{{{\mathrm{PM}}}_{2.5}}(0) = 76$, we can obtain

Substitute the values above into model (4), respectively, and by using the Markov Chain Monte Carlo (MCMC) method, we can get

1) $\beta\, _{{{\mathrm{CO}}}}(t) = 1.6131\sin\bigg(0.0088\, t+0.9976\bigg).$

2) $\beta\, _{{{\mathrm{PM}}}_{10}}(t) = 1.6135\sin\bigg(0.0088\, t+0.9972\bigg).$

3) $\beta\, _{{{\mathrm{PM}}}_{2.5}}(t) = 1.6135\sin\bigg(0.0088\, t+0.9972\bigg).$

4.4.   Model fitting and prediction performance

In numerical experiments, the weekly cumulative number of influenza cases and the values of CO, PM$_{10}$ and PM$_{2.5}$ of Huaian are respectively used to fit the number of infected cases. The model fits the reported cases in Huaian well generally (see Figure 4). Also, we find a decreasing trend in the projection results for the first twelve weeks in 2020 (see Figure 4). It should be noted that the model does not consider any interventions.

To quantify the model's simulation and prediction performance, denoted by the error between the statistical data $A = (a_1, a_2, \cdots, a_m)$ and the model values $B = (b_1, b_2, \cdots, b_m)$, we employ four statistical indices:

1) the correlation coefficient (CC):

where $\bar{a} = \frac{1}{m}\sum\limits_{i = 1}^{m}a_i,$ $\bar{b} = \frac{1}{m}\sum\limits_{i = 1}^{m}b_i.$

2) absolute error (AE):

3) root mean square error (RMSE):

4) the distance between indices of simulation and observation (DISO) ^[33,34,35]:

The performance is evaluated by the data from the 1$^{{{\mathrm{st}}}}$ week of 2019 to the 12$^{{{\mathrm{th}}}}$ week of 2020, and the results of CC, AE, RMSE and DISO are displayed in Tables 3–5, respectively.

Figure 4 and Tables 3–5 show that the fitting/predicting results of our model with CO, PM$_{10}$ and PM$_{2.5}$ have a certain reliability and rationality; hence, model (4) can well capture the history of influenza transmissions in Huaian, China.

Furthermore, from Table 2, we can also know that PM$_{2.5}$ is significantly related to CO and PM$_{10}$, and the Spearman's rank correlation coefficients are 0.73 and 0.86, respectively. Hence, for simplicity, in the remainder, as an example, we only focus on the effect of PM$_{2.5}$ on the spreading of influenza and the corresponding control measures.

4.5.   Controlling the prevalence of influenza

In this subsection, we will focus on the factors of controlling the prevalence of influenza. As an example, in Figure 5, we show the relations between the cumulative solutions of model (4) and parameters $\phi$, $\sigma$ and $\alpha$. In Figure 5a, if we decrease $\phi^* = 0.365$ by 5%, i.e., $\phi$ decreases from $\phi^*$ to $0.95\phi^* = 0.34675$, the solutions of model (4) increase; meanwhile, if we increase $\phi^*$ by 5%, i.e., $\phi$ increases from $\phi^*$ to $1.05\phi^* = 0.38325$, the solutions of model (4) decrease. In Figure 5b, if we decrease $\sigma^* = 0.485$ by 5%, i.e., $\sigma$ decreases from $\phi^*$ to $0.95\sigma^* = 0.46075$, the solutions of model (4) increase; meanwhile, if we increase $\sigma^*$ by 5%, i.e., $\sigma$ increases from $\sigma^*$ to $1.05\sigma^* = 0.50925$, the solutions of model (4) decrease. In Figure 5c, if we decrease $\alpha^* = 0.15$ by 10%, i.e., $\alpha$ decreases from $\alpha^*$ to $0.9\alpha^* = 0.135$, the solutions of model (4) decrease; meanwhile, if we increase $\phi^*$ by 10%, i.e., $\alpha$ increases from $\alpha^*$ to $1.1\alpha^* = 0.165$, the solutions of model (4) increase. Hence, in order to control the prevalence of influenza, we must increase the vaccination rate $\phi$ and the recovery rate $\sigma$ and decrease the rate of the vaccine wearing off $\alpha$.

4.5.1.   The effect of PM$_{2.5}$

In the fourth equation of model (4), parameters $\rho$ and $\xi$ are closely related to the effect of the AAP (e.g., PM$_{2.5}$) on the incidence of ILI. The numerical results are shown in Figure 6. In Figure 6a, if we decrease $\rho^* = 1$ by 50% and 75%, respectively, the solutions of model (4) decrease. In Figure 6b, if we decrease $\xi^* = 0.5$ by 50%, i.e., $\alpha$ decreases from $\xi^*$ to $0.5\alpha^* = 0.5$, the solutions of model (4) increase; meanwhile, if we increase $\xi^*$ by 50%, i.e., $\xi$ increases from $\xi^*$ to $1.5\xi^* = 1.5$, the solutions of model (4) decrease. These results show that the decrease of the uptake coefficient $\rho$ and the increase of the depletion rate $\xi$ are useful to control the spread of influenza.

On the other hand, in model (2), the term $mC(t)$ describes the impact of the AAP (e.g., PM$_{2.5}$) on the transmission rate. In Figure 6c, we show the relations between the solutions of model (4) with $m$. Obviously, $m > 0$ increases the values of transmission rate $\lambda(C)$, which results in increases in the solutions of model (4). In addition, when increasing $m^* = 10^{-10}$ by 50%, i.e., $m$ increases from $m^*$ to $1.5\times10^{-10}$, the solutions of model (4) increase. Hence, the decrease of PM$_{2.5}$ (measured by $m$) can efficiently suppress the influenza outbreak.

4.5.2.   The effect of the baseline rate

It is widely known that the transmission rate $\beta(t) = b_1\sin(b_2 t+b_3)$ plays an important role in epidemic model (4), where $b_1$ is called the baseline rate. In Figure 7, we give the relations between the cumulative solutions of model (4) with the baseline rate $b_1$, which shows that the cumulative influenza cases increase with the increase of $b_1$, and the decrease of the baseline rate $b_1$ is beneficial to controlling the incidence of influenza.

5.   Conclusions and discussions

It is now widely believed that ambient particulate air pollution is a key factor in the incidence of influenza. More different from existing research ^{[3,4,13,14,15,16,17,18]}, in this work, we have derived an AAP-dependent dynamical model that incorporates partial immunity and vaccination, and the infection rate with an AAP-dependent linear effect, and we proved that the basic reproduction number $\mathcal{R}_0$ can be used to govern the threshold dynamics of model (4): If $\mathcal{R}_0 < 1$, the influenza disease will go extinct, while if $\mathcal{R}_0 > 1$, model (4) is uniformly persistent, i.e., the influenza disease is always endemic. Based on the statistical data, via numerical simulations, we find that, in order to control the prevalence of influenza in Huaian, China, we must raise the vaccination rate $\phi$ (see Figure 5(a)) and the recovery rate $\sigma$ (see Figure 5(b)) and reduce the rate of the vaccine wearing off $\alpha$ (see Figure 5(c)), the uptake coefficient $\rho$ (see Figure 6(a)), the depletion rate $\xi$ (see Figure 6(b)), the effect coefficient of the AAP on transmission rate $m$ (see Figure 6(c)) and the baseline rate $b_1$ (see Figure 7). The theoretical and numerical results show that the dynamics of model (4) provides a profile of the prevalence of influenza in Huaian, China and gives us some useful information of controlling the influenza.

It is worthy to note that, decreasing/increasing the values of the parameters in model (4) or model (10) is complex system engineering, involving the government's intervention strategies and human behavior. In reality, what we can do is change our behavior, especially, changing our traveling plan and staying at home to reduce the contact rate $\beta(t)$ and thus decrease the incidence rate $\lambda(C)$, or increase the close-contact distance and wear protective masks to reduce the influence of the AAP (measured by $m$) on the influenza transmission.

Acknowledgments

The authors would like to thank the anonymous referees for helpful suggestions and comments which led to improvements of the original manuscript. This research was supported by the National Natural Science Foundation of China (Grant numbers 12071173 and 12171192) and Huaian Key Laboratory for Infectious Diseases Control and Prevention (HAP201704).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.