Stochastic dynamics of an SIR model for respiratory diseases coupled air pollutant concentration changes

Industrial development has made air pollution increasingly severe, and many respiratory diseases are closely related to air quality in terms of infection and transmission. In this work, we used the classic stochastic susceptible–infectious–recovered (SIR) model to reﬂect the spread of respiratory disease, coupled with the diﬀusion process of air pollutants to the infectious disease model, and we investigated the impact of various environmental noises on the process of disease transmission and air pollutant diﬀusion. The value of this study lies in two aspects. Mathematically, we deﬁne threshold R s 1 for extinction and threshold R s 2 for persistence of the disease in the stochastic model ( R s 2 < R s 1 ) when the parameters are constant, and we show that (i) when R s 1 is less than 1, the disease will go to stochastic extinction; (ii) when R s 2 is larger than 1, the disease will persist almost surely and the model has a unique ergodic stationary distribution; (iii) when R s 1 is larger than 1 and R s 2 is less than 1, the extinction of the disease has randomness, which is demonstrated through numerical experiments. In addition, we derive the exact expression of the probability density function of the stationary distribution by solving the corresponding Fokker–Planck equation under the condition of disease persistence and analyze the eﬀects of random noises on stationary distribution characteristics and the disease extinction. Epidemiologically, the change of the concentration of air pollutants aﬀects the conditions for disease extinction and persistence. The increase in the inﬂow of pollutants and the increase in the clearance rate have negative and positive impacts on the spread of diseases, respectively. We found that an increase in random noise intensity will increase the variance, reduce the kurtosis of distribution, which is not conducive to predicting and controlling the development status of the disease; however, large random noise intensity can also increase the probability of disease extinction and accelerates disease extinction. We further investigate the dynamic of the stochastic model, assuming that the inﬂow rate switches between two levels by numerical experiments. The results show that the random noise has a signiﬁcant impact on disease extinction. The data ﬁtting of the switching model shows that the model can eﬀectively depict the relationship and changes in trends between air pollution and diseases


Introduction
With the development of industrial civilization and human society, people have enjoyed affluent material life and convenient lifestyle, but this also has an adverse effect on the environment.Industrial activities produced large amounts of pollutants emitted into the air [1,2].Coal-dominated energy structure increases the concentration of sulfur dioxide (SO2), nitric oxide (NO), ozone (O3), and carbon monoxide (CO) and suspended particulate matter (PM) in the air [3].The increasing concentration of these air pollutants has raised concerns about their impact on public health.The air quality index (AQI), which is defined based on the concentration of five main air pollutants (CO, SO2, PM10, O3, and NO2), is used to inform the public about the severity of air pollution [4].AQI is a simple and understandable method to measure the impact of air quality on human health [5].The higher the AQI value, the greater the level of air pollution and the greater health concern [6].
Air pollution is an important risk factor for respiratory health.Take common pollutant PM as an example, inhalation of PM into the respiratory tract poses serious harm to the human body.The size of particulate matter determines where it will deposit in the respiratory tract: fine particles less than 10 mm diameter (PM10) mainly deposit in the upper respiratory tract, while fine particles less than 2.5 mm diameter (PM2.5)can reach the alveoli [7].Some nanoparticulates might even be able to disseminate into the systemic circulation after inhalation [8].There is sufficient evidence to confirm that air pollutants can cause cardiovascular and respiratory diseases and even increase the risk of death [9,10].Shortterm exposure to air pollutants is closely related to chronic obstructive pulmonary disease, cough, shortness of breath, wheezing, asthma, respiratory disease.The long-term effects associated with air pollution are chronic asthma, pulmonary insufficiency, cardiovascular diseases, and cardiovascular mortality [11][12][13].Therefore, exploring the associations between air pollution and health is of great importance.
There are many studies on air pollution and linked respiratory diseases.These studies have found that air pollution exposure is associated with an increased risk of respiratory symptoms, hospitalizations, and mortality.The relationship between air pollution and respiratory diseases has been demonstrated by using statistic and epidemiological methods [14,15].Specifically, Wang et al. [16] revealed the analysis of air pollution and the risk of lung cancer through a classification and regression tree model.The risk is often higher for vulnerable populations, such as children, the elderly, and those with pre-existing respiratory conditions.By using age-period-cohort analysis, Zhao et al. [17] explored the key factors including age, sex, and immigration status on cancer mortality.Additionally, they predicted the mortality rates from 2022 to 2030.This method has great inspiration for the study of air pollution and respiratory diseases.Based on the data on concentration of air pollution and daily emergency hospital admissions, Perez et al. [18] revealed that a significant increase of respiratory system problems had been detected as the PM concentrations increased by statistical analysis.Some studies also have found that elevated levels of air pollution increase the risk of mortality and morbidity [19][20][21][22].There are some studies that focus on forecasting models for air quality [23,24].The prediction level of air pollution can provide a basis for the government to take efficient response measures.However, data-based prediction models can only make short-term predictions, and long-term predictions may have significant deviations due to the influence of many fac-tors.From the perspective of dynamics, differential equation models can better interpret the internal mechanisms of respiratory diseases caused by air pollution and thus depict and predict long-term development trends [25][26][27].
In the study of the diseases related to air pollution and based on differential dynamics model, it can be divided into a deterministic model and a stochastic model according to whether stochastic factors are considered.Ganesh et al. predicted the concentrations of indoor pollutants by a physics-based model, which captures the actual scenario of air pollutant flow, and formulated a dynamic optimization problem to calculate the optimal plan for indoor air quality control [28].Tang et al. established a susceptibleexposed-infectious-recovered (SEIR) model with AQI-dependent incidence rate to characterize the seasonal changes in air pollution and their impact on the transmission of related respiratory diseases [29].On the basis of Tang's work, He et al. considered the random factors encountered during the inflow and removal of air pollutants and the process of disease transmission, and they established a corresponding stochastic model with change points, which better fit the data [30].The stochastic model of infectious diseases coupled with air pollution can well describe the relationship between them.To study how changes in the concentration of air pollutants affect the spread of respiratory diseases, it is necessary to explore the dynamic behavior of the coupled model.In order to explore the impact of random noise on the dynamics of disease transmission, He et al. conducted theoretical analysis on a coupled stochastic SIS model and demonstrated the impact of random noise on disease outbreaks and extinctions [31].The method used for analyzing the coupled susceptible-infectious-susceptible (SIS) model is to analyze the properties of one-dimensional stochastic differential equation (SDE) after dimensionality reduction.The model has an assumption that the total population is a constant, but population dynamics should be taken into account when considering long-term impacts.Therefore, further research is needed on the dynamics analysis of more complex coupled infectious disease models.
Considering that some respiratory diseases have immunity after infection, this paper uses the SIR epidemic model to describe the transmission process of the disease.On this basis, we establish a stochastic coupled model driven by air pollution.The study of stochastic epidemic models have attracted the attention of many researchers.Random factors are often not negligible, and the random environment is closer to the real situation.In theory, the addition of random factors makes the dynamics of the system different from the results obtained by the deterministic model [32][33][34][35].Du et al. [36] studied the stochastic SIR equation with general incidence functional responses and assumed that natural death rates and the incidence rate are perturbed by white noises.They provided a threshold number for the stochastic epidemic SIR model, and proved that the solution has a unique invariant measure and obtained the ergodic property when the threshold number is larger than 1.Additionally, there are studies on stochastic models that can obtain the probability density function of a stationary distribution when the disease persists [37][38][39].This plays an important role in exploring the properties of stochastic models and further analyzing the state of disease persistence.Based on the theoretical results of a stochastic SIR model, we explore the dynamic properties of the stochastic coupled model in this paper, hoping to identify key factors controlling the spread of respiratory diseases induced by air pollution.

Model derivation
The transmission process of respiratory diseases among people in polluted air environments is closely related to the concentration of air pollutants, and the transmission rate will increase with the increase of pollutant concentration.In this study, we use differential equations to characterize the changes in the number of infected individuals and the concentration of air pollutants.We couple these two processes to establish a transmission system for respiratory diseases in air-polluted environments.We use a classical SIR model to characterize the infection and transmission process of diseases [40].That is, the population is composed of classes of individuals that are susceptible S(t), infectious The birth and death of the population cannot be ignored in the long time scale.Assume that the constant birth rate of the population is and the densitydependent death rate is d.The mortality rate of infected with diseases is α, and the per capita recovery rate of the infected individuals is γ .The concentration change of air pollutants is determined by the inflow and outflow rates of pollutants, and the differential equation reflecting changes in pollutant concentration is given in [31].F(t) denotes the level of air pollution at time t, and the time-dependent inflow and clearance rate of air pollutants are denoted by c(t) and μ(t), respectively.The diagram of the coupled system is shown in Fig. 1.The pollution-dependent incidence rate is bilinear incidence rate β(t)SI.
The probability of transmission is closely related to the severity of air pollution.As mentioned in [31], the transmission probability can be defined as a function of the concentration of air pollutants, i.e., β(t) = f (F(t)), where f (x) ∈ R is a general function, and its specific form can be defined according to the model assumptions.The coupled ordinary differential equation is as follows: Intuitively, the transmission rate is positively correlated with the concentration of air pollutants.Here, we assume that the transmission probability of diseases is proportional to the concentration of air pollutants, and the constant β is a baseline transmission coefficient in the model.We only focus on the changes in the number of infected individuals, and the equation of R(t) is independent of other equations.Under these hypotheses, we Figure 1 Flow diagram of the stochastic epidemic system infected by air pollution consider the following three-dimensional system: Although deterministic models make it easier to understand the differential changes of variables, neither changes in the number of infected cases nor changes in the concentration of air pollutants are smooth, and they are more or less influenced by stochastic factors in the environment.We are interested in investigating the effect of environmental noise on the dynamic.The stochastic differential equation gives us a more realistic way to model the dynamics affected by random perturbations.A stochastic version of the coupled SIR epidemic model can be presented by considering the effects of environmental fluctuations.We assume that the change of each variable is affected by environmental noise, which is proportional to the size of the variable.We put three linear perturbations that depend on the variables S, I, and F in system (2.2).By introducing randomness in the form of white noise into the model, we can establish the corresponding stochastic differential equation model of model (2.2) as follows: where dB i (t) = B i (t + dt) -B i (t) (i = 1, 2, 3) are the increment of a standard Brownian motions and σ i (i = 1, 2, 3) are all real constants; denote the intensity of environmental fluctuations for compartments S(t), I(t) and the pollutant concentration F(t), respectively.In order to comprehensively explore the impact of changes in air pollutants on the spread of respiratory diseases, we first explore the dynamic behavior of the system when the parameters are relatively stable, that is, when the inflow rate and the clearance rate are constant.Next, we focus on the dynamics of the coupled stochastic SIR model driven by air pollution with c(t) = c and μ(t) = μ.

Stochastic extinction and persistence of the disease
Undoubtedly, whether the disease tends to become extinct is a key concern in the research of epidemic models.As given in Theorem A.1, the persistence of the disease for deterministic model (2.2) depends on the asymptotic stability of the positive equilibrium.However, this technique breaks down for random model (2.3), which is because the model does not have the equilibrium.In this section, we use the theory of stochastic differential equation (SDE) to find the threshold conditions that govern disease dynamics for model (2.3), so as to further show the stochastic extinction and stochastic persistence of disease in model (2.3).From now on, w(t) stands for 1 t t 0 w(s) ds.

Stochastic extinction
To demonstrate the stochastic extinction of diseases, we provide the threshold parameter 2 ) . From the threshold expression of disease extinction, it can be seen that the threshold will increase as the pollutant inflow rate becomes larger and will decrease as the clearance rate becomes larger.This indicates that the parameters that cause changes in the concentration of air pollutants can also affect the conditions for disease extinction.By comparing the threshold R s 1 of disease extinction in the stochastic model (2.3) with the basic reproductive number R D 0 in the deterministic model (2.2), it can be found that only the change in noise intensity of I(t) has an impact on the extinction threshold of the disease, while random noise caused by change in the concentration of air pollutants has no impact on the threshold.An increase in noise intensity will lower the threshold for disease extinction.

Stochastic persistence
This subsection is intended to motivate our investigation about stochastic persistence of diseases by defining . Theorem 3.2 When R s 2 > 1, the solution (S(t), I(t), F(t)) for model (2.3) with the initial value (S(0), I(0), F(0)) ∈ R 3 + has the following property: In other words, the infectious individuals will persist if R s 2 > 1 holds.
Proof Define the Lyapunov function: where m 1 and m 2 are positive constants that will be determined later.
Hence, we obtain 2 )(μ + 2 ) 2 )(μ + Then we see that It follows from (3.2) that integrating the above inequality from 0 to t and then dividing by t on both sides may lead to where ϑ(t) := t 0 (-σ 1 dB 1 (s)m 1 σ 2 dB 2 (s)m 2 σ 3 dB 3 (s)).Using Lemma A.3 and the strong law of large numbers for martingale, we obtain that which proves the theorem.
which implies that the inability to control stochastic dynamics with only one threshold parameter is due to the presence of σ 1 and σ 3 ; when which implies that the introduction of noises reduces the value of threshold parameters.
Remark 3. 4 We further realize that the stochastic dynamics of model (2.3) are governed by two threshold parameters R s 1 and R s 2 , which is similar to the dynamics of the stochastic SIS model with media coverage in [41] (see Theorem 3.5 and Theorem 4.5 in [41]).To better understand the dynamics of stochastic SIF model (2.3), we provide a schematic of stochastic extinction and stochastic persistence (as shown in Fig. 2), which is similar to Fig. 8 in [41].
Remark 3.5 Theorem 3.1 sketches a conclusion that under the condition R s 2 < R s 1 < 1, the infected individuals for model (2.3) will disappear with probability one, which means that almost all solutions of model (2.3) will tend to the disease-free equilibrium E 0 (Region I in Fig. 2).It should be noted that as long as the noise intensity σ 2 is too large, the outbreak of the disease will be controlled.Theorem 3.2 indicates that if the noise intensities σ 1 , σ 2 , and σ 3 are relatively small, resulting in R s 1 > R s 2 > 1, then the disease will be widely prevalent (Region II in Fig. 2).Remark 3.6 There exists an example where R s 2 < 1 < R s 1 , which is due to the fact given in Remark 3.3 that R s 2 < R s 1 (see Region III in Fig. 2).However, it is regrettable that the disease dynamics of stochastic model (2.3) failed to be demonstrated in this case and will be further shown through numerical simulations in Sect. 4.
Different from the threshold of disease extinction, the expression of the stochastic persistence threshold we obtained contains three noise intensities, which means that random perturbations of the three variables in the stochastic system will affect the persistent state of disease spread.It can be seen from the expression for R s 2 that an increase in any of the three noises σ 1 , σ 2 , σ 3 will lower the persistence threshold of the disease.This indicates that random noises during variable changes in the system can have adverse effects on disease outbreaks, and large noise can cause disease outbreaks.The two thresholds of disease persistence and extinction illustrate the impact of noise on disease transmission from two different perspectives.Same as the threshold R s 1 , the threshold R s 2 will also increase as the pollutant inflow rate becomes larger and will decrease as the clearance rate becomes larger.However, we cannot draw the conclusion that an increase in the inflow rate of air pollutants is beneficial to the spread of the disease because the dynamics between the two thresholds are still uncertain and will be further explored in numerical experiments.

Stationary distribution
In this subsection, we focus on whether stochastic model (2.3) has a stationary distribution by applying Hasminskii's ergodic theory provided in Lemma A.4.

3) has a unique stationary distribution ν(•).
Proof The diffusion matrix associated with model (2.3) is where X = (S, I, F).
where W 1 is defined in Theorem 3.2, and Here, M > 0 is a sufficiently large constant satisfying It is easy to check that where Then we define a nonnegative function W : R 3 + → R as follows: Simple calculation yields that By Itô's formula, we get Choosing }, we have Similarly, we have Choosing Thus, Next, we construct a compact subset D such that condition (B.2) in Lemma A.4 holds.Define the bounded closed set where is a sufficiently small constant.We choose satisfying the following condition: where For convenience, we divide R 3 + \D into six domains: We will prove that L W (S, I, F) ≤ -1 on R 3 + \D, which is equivalent to showing it on the above six domains.
Case 2: If (S, I, F) ∈ D 2 , then one can obtain that Then we can derive that Then we can derive that Case 6: If (S, I, F) ∈ D 6 , then one can obtain that Hence, one conjectures that Therefore, according to Lemma A.4, one can get that model (2.3) has a unique stationary distribution ν(•).This completes the proof.

Probability density function
The analysis in the previous subsection makes it clear that if 3) has a unique stationary distribution.From this, we wonder what the exact expression of the probability density function is.Next, we will have a detailed discussion on this issue.
Letting x = ln S, y = ln I, z = ln F, we obtain the equivalent model of model (2.3): 3) admits a quasi-steady state E * = (S 1 , I 1 , F 1 ) = (e x * , e y * , e z * ), where 2 ) βc , , Obviously, E * is consistent with positive equilibrium E * of deterministic model (2.2) when where x * = ln S 1 , y * = ln I 1 , z * = ln F 1 .Then the linearized system of (3.5) is as follows: where In the following, using Lemmas A.5 and A.6, we show the local probability density function of model (2.3) near the quasi-steady equilibrium E * , which is presented in Theorem 3.8 through denoting ) be a solution of (3.7) with any initial value 1 and H = 0 hold, then there is a probability density function (v 1 , v 2 , v 3 ) around the quasi-steady state E * , which can be expressed in the following form: , where and other parameters are provided in the proof of this theorem.
For simplicity, we put the proof of Theorem 3.8 in Appendix B.
Here, we can now rephrase Theorem 3.8 as follows.
Theorem 3.9 Let (S(t), I(t), F(t)) be a solution of (2.3) with any initial value (S(0), 1 and H = 0 hold, then there is a probability density function (S, I, F) around the quasi-steady state E * , which can be expressed in the following form: ) T , and is defined in Theorem 3.8.
Remark 3.10 The marginal density functions of (S, I, F) with respect to S, I, and F are as follows: where 11 , 22 , and 33 are defined in the proof of Theorem 3.8 in Appendix B.
We obtain the explicit form of the probability density function of the stationary distribution of the solution for stochastic model (2.3) when the disease persists.The probability density function is very helpful for us to study the properties of disease transmission when it reaches stability.By analyzing the mean and variance of the distribution, the outbreak level and controllability of the disease epidemic can be determined.The skewness and kurtosis of distribution can be used to determine the extreme conditions of disease development.Based on the form of probability density functions, we will identify the characteristics of the distribution and the impact of different parameters on the distribution by numerical simulations.

Numerical simulations 4.1 Validation of the theoretical results
The results in Theorem 3.1 show a threshold R s 1 for the disease stochastic extinction, i.e., the disease will tend to extinction when R s 1 < 1.The results in Theorem 3.2 indicate that the disease will persist when R s 2 > 1.To verify our analytical findings, we perform the simulation results of stochastic model (2.3) and its corresponding deterministic model by using the EM method under the condition in Theorems 3.1 and 3.2, respectively.As stated in Remarks 3.3 and 3.5, we can see that when R D 0 < 1, the diseases for both deterministic model (2.2) and stochastic model (2.3) will tend to die out.Hence, we mainly focus on the dynamics in the case of R D 0 > 1.Our selection of parameters refers to the parameter estimation results based on the actual data in [30], and some parameter values are set as follows: In this case, it can be calculated that R D 0 = 2.327 > 1, which means that deterministic model (2.3) has a disease-free equilibrium E 0 = (875, 0, 25) and a globally asymptotically stable endemic equilibrium E * = (376, 21.23, 25).As shown in Fig. 2, it can be divided into three regions according to the size relation between R s 1 and R s 2 .Next, we will select parameters satisfying three different regions to investigate the dynamics for model (2.3).Specifically, these experiments are performed for the parameters found in (4.1), except for σ j (j = 1, 2, 3), which are allowed to vary.

The stochastic dynamics in the case
Numerical experiments are initialized from the dynamic behavior of model (2.3) with (σ 1 , σ 2 , σ 3 ) = (0.12, 0.7, 0.12).The parameter values guarantee the condition of Theorem 4.1.At this point, we have which is located in region I in Fig. 2. It corresponds to the case of R s 2 < R s 1 < 1 < R D 0 .The simulated results shown in Fig. 3 are both the solutions of stochastic model (2.3) and the corresponding solutions of deterministic systems (2.2).The red, green, and blue curves represent the solutions of S(t), I(t), and F(t) for stochastic model (2.3) changing with time.The black curves are the corresponding solutions of deterministic model (2.2).We can see that the corresponding solutions of deterministic system stabilize at a constant greater than zero and consequently the disease persists, while the solution I(t) of stochastic system goes to zero after a finite time and consequently the disease will become extinct with probability one.The numerical simulation results are consistent with the conclusion Figure 3 The simulations of stochastic model (2.3) and its corresponding deterministic system under case (Q1) with (σ 1 , σ 2 , σ 3 ) = (0.12, 0.7, 0.12) of Theorem 3.1.In addition, the difference between the solutions of the stochastic model and the deterministic system indicates that random noise may not be conducive to the spread of disease, but rather leads to its extinction.

The stochastic dynamics in the case
Here, our experiment aims to confirm the conclusion in Theorem 3.2.We fix the parameters to satisfy the conditions of case 1 < R s 2 < R s 1 < R D 0 by decreasing the noise intensities σ j (j = 1, 2, 3).Taking (σ 1 , σ 2 , σ 3 ) = (0.01, 0.01, 0.01), then the threshold parameters which is located in region II in Fig. 2. Furthermore, we have σ 2 1 2 = 0.09 = 0.According to Theorem 3.7, Theorem 3.9, and Remark 3.10, we know that model (2.3) admits a unique stationary distribution and the distribution has a probability density function around the quasi-steady state E * = (376.In Fig. 4(a-c), we show a single stochastic realization of stochastic model (2.3) together with the solution of its corresponding deterministic model.The simulated results shown in Fig. 4 reveal that the corresponding solutions of deterministic system stabilize at a constant greater than zero and the solutions of stochastic model have a positive stationary distribution, which reveals that the diseases will be almost surely prevalent when R s 2 > 1. Figure 4(d-f ) shows the frequency histogram fitting curves of (S(t), I(t), F(t)) at the iteration times equal to 200,000 and the theoretical marginal density functions ϒ(S), ϒ(I), and ϒ(F), respectively.These results indicate that the fitting curves almost coincide with the corresponding theoretical marginal density functions.

The stochastic dynamics in the case
To investigate the dynamic behavior of model (2.3) in the gap between the extinction and persistence conditions we theoretically provide, we consider the case that the parameters satisfy R s 2 < 1 < R s 1 , which is in region III of Fig. 2. When parameters other than noise intensity are fixed, the value of R s 1 is determined by σ 2 , and the value of R s 2 is determined by three noise intensities σ 1 , σ 2 , and σ 3 .To better reflect the changes in R s 1 and R s 2 using noise intensity, we fixed parameter σ 1 and parameter σ 3 for the dynamic analysis of the stochastic model, respectively.Firstly, we investigate the impact of noise intensity σ 2 on disease transmission.We choose a parameter interval of [0.05, 0.1] for σ 1 to ensure that the value of R s 1 is greater than 1.The interval of parameter σ 2 is set to [0.001, 0.2] so that the value of R s 2 is less than 1.The surface in Fig. 5(a) shows the variation of R s 1 and R s 2 with changes in noise intensity.Figure 5(b) shows the persistence and extinction states of the disease when the parameters satisfy the condition R s 2 < 1 < R s 1 .The green region represents the parameter range of disease extinction, while the magenta region represents the parameter range of disease persistence.
From Fig. 5(b), we observe that when R s 2 < 1 < R s 1 , the disease will become extinct under some noise values and will persist at others.From the color distribution in the figure, it (b) When parameter values of σ 1 and σ 3 are chosen from the green point, the solution of I(t) will converge to 0. When parameter values are chosen from the magenta point, the solution of I(t) will stabilize at a nonzero state can be seen that the space of disease persistence and extinction in region III is irregular.
Similarly, we investigate the impact of noise intensity σ 3 on disease transmission.The value of σ 1 is fixed and the interval of parameter σ 3 is set to [0.6, 1.5] to ensure that R s 2 is less than 1. Figure 6(a) shows the variation of R s 1 and R s 2 with changes in noise intensity σ 2 and σ 3 .Figure 6(b) shows the persistence and extinction states of the disease when σ 2 and σ 3 change.The color distribution in Fig. 6(b) also shows that the conditions of disease persistence and extinction in region III are not clearly defined but chaotic.

The effect of noise intensities on the stochastic dynamics
The above theory has proved that when the parameters satisfy certain conditions, the number of infected cases for the disease will reach a stationary distribution after a period of time.There are three random noise parameters in the stochastic model (2.3), we explore the impact of different noise parameters on the stationary distribution in this section.Figure 7(a-c) shows simulations of ϒ(I) as three noise intensities σ 1 , σ 2 , and σ 3 change, respectively.The other parameters are fixed as R s 2 > 1.We compare the changes with noise by calculating the four characteristics of stationary distribution: mean, variance, kurtosis, and skewness kurtosis.In Fig. 7(a), we can see that the impact of noise intensity σ 1 changes on the mean is minimal, while the variance increases with the increase of noise.After calculation, it is found that the kurtosis decreases from 15.95 to 2.73 as σ 1 increases, and the skewness decreases from 3.71 to 1.1.It indicates that the increase of random noise of S(t) will make the stationary distribution for I(t) more diffuse.In Fig. 7(b), both the mean and variance of I(t) change little as noise intensity σ 2 increases, and the kurtosis decreases from 4.67 to 4.51 as σ 2 increases, and the skewness decreases from 1.73 to 1.67.It can be seen from Fig. 7(c) that, similar to the effect of noise intensity σ 2 , noise intensity σ 3 has little effect on the stationary distribution for I(t).
When diseases can become extinct is a meaningful question that naturally comes to our mind.Based on the fact revealed in Remark 3. and the average extinction time is 9.0.Thus, we assert that the increase in noise intensity σ 2 triggers an exciting result, which is that the extinction time is significantly shortened.
In the previous numerical simulation results, we found that when R s 1 is greater than 1 and R s 2 is less than 1, the persistence and extinction of the disease are uncertain.But as can be seen from Fig. 5, the distribution of extinction is not completely random.We want to see how the extinction probability of diseases changes as R s 2 changes between 0 and 1 when R s 1 satisfied is larger than 1.The variation of R s 2 > 1 is mainly determined by the noise intensity σ 1 and σ 3 when other parameters are fixed.Therefore, we separately inves- tigated the extinction probability of diseases under changes in random noise intensities σ 1 and σ 3 from small to large.The probability of disease extinction is calculated by the frequency at which the solution I(t) of stochastic model (2.3) tends to 0 in 1000 simulations.The results in Fig. 9(a) show that as the noise intensity σ 1 increases, the probability of disease extinction increases.The same increase in noise intensity σ 3 will also increase the probability of extinction.This indicates that large random noise is beneficial for the elimination of diseases.

Stochastic environments with varying pollutant inflow rates
In the previous section, we discussed the dynamics of stochastic model (2.1) when the inflow rate is constant, and we found that when the system satisfied the conditions that R s 2 > 1, the disease will persist and its stable state obeys a stationary distribution.In real life, the levels of air pollution are not consistent throughout the year but show seasonal changes.For the rate of inflow of pollutants, we know that pollutants are not permanent and that their concentrations change with seasonal factors.For example, the emission of pollutants will significantly increase in winter due to burning fuels for heating systems.To further investigate the impact of air quality differences in different seasons of the year on disease transmission, we assume that the inflow rate in model (2.1) has the following form: -1) , T (2n) ), n = 1, 2, 3, . . ., (5.1) where T i is the change point of time.Without losing generality, we assume that c 1 > c 2 .This means that the inflow rate constantly switches between c 1 and c 2 for a long time, switching once in an interval [T (2n-2) , T (2n) ], n = 1, 2, 3, . . . .Within the interval [T (2n-2) , T (2n-1) ) and [T (2n-1) , T (2n) ), the inflow rate remains constant, indicating that its dynamic behavior is consistent with the situation we analyzed in the previous section.Therefore, we will explore through numerical simulation experiments how the disease spreads under such a switching pattern.According to the relationship between c 1 and c 2 , there exists a basic reproductive number corresponding to the deterministic model and two thresholds for the stochastic model in each switching time interval.We use superscripts c 1 and c 2 to distinguish the parameter values of two switching intervals.From the expression of basic reproductive number R 0 and the thresholds R s 1 and R s 2 , it is easy to determine that R Dc There are many possibilities for the relationship between these six parameters and the size of 1.We list various cases of their size relationship in Table 1.For the sake of brevity and focus, certain intermediary scenarios have been omitted from our consideration.We conduct 1000 stochastic simulation experiments for each scenario to observe whether the disease becomes extinct after prolonged transmission.The results of the stochastic model and the corresponding deterministic model are also presented in Table 1.From the results, it can be seen that the final state of the deterministic model and the stochastic model is consistent only when the selected parameters guarantee all values are greater than 1.However, when the threshold for stochastic extinction within the lower range of pollution flows is less than 1, the disease may persist or become extinct.When the threshold for stochastic persistence and extinction within the lower range of pollution flows is less than 1, the disease will become extinct.This indicates that switching air pollution levels can help control diseases.

Table 1 The extinction or persistence of the disease with inflow rate switching
The size relationship of different thresholds The final state of the disease for deterministic model (2.1) The final state of the disease for stochastic model (2.3) Persistence Extinction We utilize the relevant disease infection dataset and air quality data from Xi'an City, Shaanxi Province, China, to identify the model.The disease infection data comes from influenza-like illness (ILI) hospitalizations data originally from the surveillance system at the Shaanxi Center for Disease Control and Prevention, Xi'an city Shaanxi province.The data for the air quality index (AQI) are obtained from the web of www.tianqihoubao.com.The time series spans November 15th, 2010 to November 14th, 2016.Parameter estimation is carried out by the least squares method, and the fitting results are shown in Fig. 10.We can see that when the inflow rate of air pollutants switches between two values, the stochastic model can present the trend of periodic changes in air pollution levels and also drive the periodic changes in the number of infected cases.

Conclusion and discussion
The gradual deterioration of air quality has led to an increase in the number of respiratory infections.The transmission of respiratory disease and diffusion of air pollutants are subject to many random factors such as individual migration, behavior changes of humans, and various climate factors.To investigate the impact of random factors on the spread of disease, we assumed that the transmission rate of the disease depends on the concentration of air pollutants F(t) and developed a coupled stochastic SIR epidemic model with multiple noise influence.We analyze fundamental properties of the stochastic model and its corresponding deterministic model when the parameters do not change over time.By applying Itô's formula and constructing the Lyapunov function, we found two thresholds R s 1 and R s 2 for disease persistence and extinction.We prove that the disease will go to stochastic extinction when threshold R s 1 is less than 1, and the model has a unique ergodic stationary distribution when R s 2 is larger than 1, which reveals that the disease will persist almost surely.In addition, we derive the exact expression of the probability density function of the stationary distribution by solving the corresponding Fokker-Planck equation.It can be determined that R s 1 > R s 2 , but we cannot theoretically show whether the disease tends to become extinct when R s 1 > 1 > R s 2 .By comparing thresholds of disease persistence and extinction in the stochastic model with the basic reproduction number of the deterministic model, the threshold of disease extinction in the stochastic model is smaller than that in the deterministic model.This means that under the same conditions, there will be a situation where the disease will break out without considering environmental noises, but the disease will become extinct if the impact of environmental noise is considered.We have further verified the theoretical results by numerical experiments.We refer to the parameter estimation results in [31] for parameter settings.Firstly, we selected a parameter that satisfies R s 1 < 1 and R D 0 > 1.From the solutions shown in Fig. 3, it was observed that the infection is eventually eliminated from the population in the stochastic model, while the disease persists in the deterministic model under this condition.When the parameter satisfies R s 2 > 1, the persistence of the disease occurred to the model in the long-term behavior.
Although the dynamics between the two thresholds cannot be determined theoretically, we have discovered interesting phenomena from numerical experiments.When R s 1 > 1 > R s 2 , we plotted the persistence and extinction of the disease as the noise intensities change (also representing changes in R s 1 and R s 2 ).According to the mixed space distribution in Figs. 5 and 6, it can be explained that the persistence or extinction of the disease is random between two thresholds, which is different from the basic reproductive number of the deterministic system.These numerical results indicate that the dynamics of the stochastic model are more complex considering the influence of noise.On the one hand, we investigated the impact of different random noises on the stationary distribution, and the results show that the noise intensity σ 1 that affects the variation of S(t) has a significant impact on the variance, kurtosis, and skewness of the stationary distribution.The larger the noise, the greater the variance and the smaller the kurtosis of distribution.It indicates that changes in the number of cases make it more difficult to control the disease.On the other hand, to explore the impact of noise on disease extinction, we calculated the change of the average extinction time of disease with the change of noise intensity σ 2 , and the results show that large noise intensity of σ 2 will accelerate the extinction of disease.When the system parameters satisfy R s 1 > 1 > R s 2 , an increase in the random noise intensities σ 1 and σ 3 would increase the probability of disease extinction.
We simulated the situation when the parameters change over time using numerical experiments.We considered the inflow rate of air pollutants switching between two different levels over time.The inflow rate is constant within every two switching time intervals, and the dynamics of the model can refer to the case where the parameters remain unchanged.Based on the relationship between the threshold conditions for persistence and extinction of diseases at two switching levels, there are six scenarios for the switched stochastic system.We investigated the persistence or extinction states of diseases in each of these six scenarios, and the results are listed in Table 1.By comparing the results with the deterministic model, the disease is persistent in the deterministic model under six scenarios, it is found that the disease in the stochastic model persist when the threshold of stochastic persistence is greater than 1 with a low switching level.When the threshold of stochastic extinction is less than 1 and the threshold persistence is greater than 1 with a low switching level, both the persistence and extinction of the disease are possible.The disease will be extinct in other scenarios.This indicates that random noise has a significant impact on the extinction of diseases.To demonstrate the applicability of the model, we use the stochastic model with inflow rate switching to fit the actual data, and the fitting results show that the switched model can well reflect the characteristics of periodic changes in air quality and disease transmission.
In summary, the study of this paper is to couple the diffusion of air pollutants in the classic SIR model, which characterize the transmission of air pollution-induced respira-tory diseases.Based on the theoretical analysis method of stochastic epidemic models, we analyzed the dynamics of this stochastic coupled model, focusing on the impact of multiple random noises on disease outbreaks and extinctions.In reality, the spread of diseases is affected by many factors.The establishment of coupling models can better reflect the interrelationships at multiple levels.Therefore, it is necessary to develop stochastic coupling epidemic model analysis.

Appendix A: Preliminaries
Let us briefly provide the dynamics of the model without considering random factors, i.e., deterministic model (2.2).Denote the basic reproduction number which governs the dynamics of SIF model (2.2).
It is easy to verify that model (2.2) has the following equilibrium: (i) Disease-free equilibrium: E 0 = ( d , 0, c μ ); (ii) Endemic equilibrium: E * = (S * , I * , F * ), here R D 0 > 1 and Next, we mainly focus on the stability of the equilibria for model (2.2).This proof is similar to the proof of Theorem 3.1 in [42] and Theorem 2.2 in [43], and we present it in Appendix B.
Let ( , F, {F t } t≥0 , P) be a complete filtered probability space with the filtration {F t } t≥0 satisfying the usual condition, i.e., it is increasing and right continuous, while F 0 contains all P-null sets.
Here, we first give the following theorem about the existence and uniqueness of global positive solution.
The proof is straightforward and goes back to the work of [31].Next, some useful lemmas for stochastic threshold dynamics are given as follows.
Lemma A.3 Let (S(t), I(t), F(t)) be the solution of model (2.3) with the initial value (S(0), I(0), F(0)) ∈ R This proof is standard, and hence we omit it.Furthermore, we introduce the theory of Hasminskii [44], which has been applied to show the existence of stationary distribution.
Let X(t) be a regular time-homogeneous Markov process in R l described by the following stochastic differential equation: The diffusion matrix is defined as follows: where X(t) is nonsingular.
Lemma A.4 [44] If there exists a bounded open domain U ⊂ R l with regular boundary , having the following properties: (B.1)In the domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix (X) is bounded away from zero; (B.2) If X ∈ R l \ U, then the mean time τ at which a path issuing from X reaches the set U is finite, and sup X∈K E X τ < ∞ for every compact subset K ⊂ R l , then the Markov process X(t) has a stationary distribution.
Finally, we provide some lemmas that are helpful to derive the expression of the probability density function.The proofs can be referred to in [45].

Lemma A.5 Consider the algebraic equation G
If a 1 > 0, a 3 > 0, and a 1 a 2a 3 > 0, then the matrix 0 is positive definite.In this case, the matrix of the form A 0 is called the standard R 1 matrix.

Appendix B
Proof of Theorem A.1 (i) In the view of model (2.2), we know that the solution of (2.2) has the following properties: From the first equation, it is easy to check that Hence, when R D 0 < 1, there exists ε > 0 such that dI(t) dt ≤ -εI(t).
According the Lyapunov-LaSalle theorem, E 0 is globally asymptotically stable.
(ii) Through summing the first two equations of model (2.2), we find that the total population N(t) = S(t) + I(t) satisfies the following equation: Thus, model (2.2) is equivalent to the following model: where .
Then the derivative of V along the solution of model (B.2) is given by Consequently, we can conjecture that where the sufficiently small constant ε and the appropriate p satisfy Then we have dV dt < 0. According to the Lyapunov-LaSalle theorem, the equilibrium Ẽ * of model (B.2) is globally asymptotically stable, which implies that the equilibrium E * of model (2.2) is globally asymptotically stable.
This ends the proof.
Proof of Theorem 3.8 We rewrite model (3.7) in the following form: Following Roozen [46], the density function (V ) = (v 1 , v 2 , v 3 ) of the quasi-stationary distribution of model (3.7) at the origin can be described by the three-dimensional Fokker-Planck equation as follows: which can be approximated by a Gaussian distribution where V * = (0, 0, 0), and Q is a real symmetric matrix, which is described by In the light of the finite independent superposition principle [47], (B.3) can be replaced by the sum of the following three equations: where To verify the positive definiteness of , we consider the corresponding characteristic equation of A as follows: where Furthermore, we can verify that a 1 a 2a 3 > 0.
Next, we focus on the positive definiteness of ( = 1 + 2 + 3 ) by the following three steps.

Figure 2
Figure 2Schematic of stochastic extinction and stochastic persistence for stochastic model(2.3)

Case 1 :
If (S, I, F) ∈ D 1 , then one can obtain that

Case 4 :
If (S, I, F) ∈ D 4 , then one can obtain that

Case 5 :
If (S, I, F) ∈ D 5 , then one can obtain that

Figure 5 Figure 6
Figure 5The space of disease extinction to persistence in region III with σ 1 and σ 2 change.(a) The gray surface is a plane with a value of 1, which is for the convenience of judging the relationship between the other two surfaces and 1. Above the gray surface is the value of R s 1 , below the gray surface is the value of R s 2 .(b)When parameter values of σ 1 and σ 2 are chosen from the green point, the solution of I(t) will converge to 0. When parameter values are chosen from the magenta point, the solution of I(t) will stabilize at a nonzero state 5, we attempt to investigate the impact of σ 2 on extinction time of the infected individuals I(t) by fixing the noise intensities σ 1 = σ 3 = 0.12.Here, the results of the extinction time are shown in the box plot of Fig. 8.When σ 2 = 0.7, the median and confidence interval of the extinction time for I(t) of model (2.3) is 34.2 (95% CI: [14.2,46.7]),and the average extinction time is 35.0; when σ 2 = 1.1, the median and confidence interval of the extinction time for I(t) of model (2.3) is 15.8 (95% CI: [4.6, 30.0]), and the average extinction time is 15.9; when σ 2 = 1.5, the median and confidence interval of the extinction time for I(t) of model (2.3) is 8.3 (95% CI: [2.5, 18.3]),

Figure 7 Figure 8 12 Figure 9
Figure 7 Probability density distribution of I(t) with three noise variations

Figure 10
Figure 10 Data fitting of for the ILI cases and AQI using stochastic model (2.3) with inflow rate switching

Theorem A. 1
(i) If R D 0 < 1, then the disease-free equilibrium E 0 is globally asymptotically stable.(ii) If R D 0 > 1, then the endemic equilibrium E * is globally asymptotically stable.
follows from the third equation of the model that lim t→∞ F(t) = c μ .