
Citation: Giuseppina Albano. Detecting time-changes in PM10 during Covid pandemic by means of an Ornstein Uhlenbeck type process[J]. Mathematical Biosciences and Engineering, 2021, 18(1): 888-903. doi: 10.3934/mbe.2021047
[1] | Marte K.R. Kjøllesdal, Gerd Holmboe-Ottesen . Dietary Patterns and Birth Weight—a Review. AIMS Public Health, 2014, 1(4): 211-225. doi: 10.3934/Publichealth.2014.4.211 |
[2] | Carla Gonçalves, Helena Moreira, Ricardo Santos . Systematic review of mediterranean diet interventions in menopausal women. AIMS Public Health, 2024, 11(1): 110-129. doi: 10.3934/publichealth.2024005 |
[3] | Hanan E. Badr, Travis Saunders, Omar Bayoumy, Angelie Carter, Laura Reyes Castillo, Marilyn Barrett . Reversal for metabolic syndrome criteria following the CHANGE program: What are the driving forces? Results from an intervention community-based study. AIMS Public Health, 2025, 12(1): 162-184. doi: 10.3934/publichealth.2025011 |
[4] | Truong Tuyet Mai, Tran Thu Trang, Tran Thi Hai . Effectiveness of germinated brown rice on metabolic syndrome: A randomized control trial in Vietnam. AIMS Public Health, 2020, 7(1): 33-43. doi: 10.3934/publichealth.2020005 |
[5] | Jürgen Vormann . Magnesium: Nutrition and Homoeostasis. AIMS Public Health, 2016, 3(2): 329-340. doi: 10.3934/publichealth.2016.2.329 |
[6] | Yu Li, Zhehui Luo, Claudia Holzman, Hui Liu, Claire E. Margerison . Paternal race/ethnicity and risk of adverse birth outcomes in the United States, 1989–2013. AIMS Public Health, 2018, 5(3): 312-323. doi: 10.3934/publichealth.2018.3.312 |
[7] | Jennifer L Lemacks, Laurie S Abbott, Cali Navarro, Stephanie McCoy, Tammy Greer, Sermin Aras, Michael B Madson, Jacqueline Reese-Smith, Chelsey Lawrick, June Gipson, Byron K Buck, Marcus Johnson . Passive recruitment reach of a lifestyle management program to address obesity in the deep south during the COVID-19 pandemic. AIMS Public Health, 2023, 10(1): 116-128. doi: 10.3934/publichealth.2023010 |
[8] | Trong Hung Nguyen, Thi Hang Nga Nguyen, Hung Le Xuan, Phuong Thao Nguyen, Kim Cuong Nguyen, Tuyet Nhung Le Thi . Nutritional status and dietary intake before hospital admission of pulmonary tuberculosis patients. AIMS Public Health, 2023, 10(2): 443-455. doi: 10.3934/publichealth.2023031 |
[9] | Ragnar Rylander . Treatment with Magnesium in Pregnancy. AIMS Public Health, 2015, 2(4): 804-809. doi: 10.3934/publichealth.2015.4.804 |
[10] | John P. Elliott, John C. Morrison, James A. Bofill . Risks and Benefits of Magnesium Sulfate Tocolysis in Preterm Labor (PTL). AIMS Public Health, 2016, 3(2): 348-356. doi: 10.3934/publichealth.2016.2.348 |
Coronavirus disease 2019 (COVID-19) remains an on-going global pandemic at present. The World Health Organization declared COVID-19 as a Public Health Emergency of International Concern (PHEIC) on January 30, 2020. According to data released by Johns Hopkins University, there are 37, 213, 592 confirmed cases and 1, 072, 959 deaths in 188 countries and regions around the world by October 11, 2020 [1]. The disease is caused by a novel coronavirus named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [2]. Although most SARS-CoV-2-infected cases have asymptomatic or mild-to-moderate diseases, around 10% of those infected may develop severe pneumonia and other associated organ malfunctions [3].
In recent investigations, many scholars studied different epidemic models of COVID-19. In Zhang's study [4], a new mathematical model (SEIRD) was proposed, which is constructed with five classes including susceptible, exposed, infected, recovered and deaths to describe the possibility of transmission in a given general population. Bhadauria et al. [5] studied the SIQ model by using the stability theory of nonlinear ordinary differential equations. Li et al. [6] developed a numerical method preserving positivity for a stochastic SIQS epidemic model. In Ref. [7], Higazy proved the existence of a stable solution for the fractional order COVID-19 SIDARTHE model. A θ−SEIHQRD model which is more relevant to COVID-19 was developed by Ramos et al. [8]. Nisar et al. [9] constructed the SIRD model and verified the correctness. Batistela et al. [10] proposed an SIRSi model of COVID-19. Especially, Paré et al. [11] proposed traditional group models, continuous-time and discrete-time versions of the models with non-trivial networks on simple SIR-based models, and supported the need for networked models through presenting a set of simulations.
In the propagation process of COVID-19, we consider that there is a time-delay from infection to recovery. At present, there are some researches on epidemic models with time-delay. Zhu et al. [12] constructed a time delay reaction-diffusion model that is closer to the actual spread of the COVID-19 epidemic with considering the time delay effect of infected persons during the spread of the epidemic. In Ref. [13], the mathematical modeling of COVID-19 fatality trends was constructed by Scheiner et al., which was based on infection-to-death delay rule. Wei et al. [14] considered the time delay from susceptible individuals to infected individuals, thus, proposed a new SVEIR epidemic disease model with time delay, and analyzed the dynamic behavior of the model under pulse vaccination. In Ref. [15], Mukherjee considered an S-I epidemic model with time delay, the time delay is the immune period and the incubation period, and gave an estimate on the length of delays for system which is stable in the absence of delays remains stable.
At the same time, many researchers analyzed the stability of the COVID-19 model. In Ref. [16], Araz dealt with a mathematical model about COVID-19 spread, and analyzed global and local stability for the considered model. Annas et al. [17] carried out the stability analysis and numerical simulation of the SEIR model on the spread of COVID-19 in the research. Besides, some scholars conducted mathematical analysis on the principle of COVID-19 infection. Almocera et al. [18] studied an in-host model, and the stability of a unique positive equilibrium point, with viral load V∗, suggested that the virus may replicate fast enough to overcome T cell response and cause infection. In Ref. [19], Samanta analyzed the stability of the proposed model to control the epidemic.
And some scholars predicted the epidemic situation. León et al. [20] proposed an SEIARD mathematical model and attempted to forecast the evolution of the outbreak. Youssef et al. [21] carried out numerical verification and predictions of the proposed SEIR model, and compared the results with the real data due to the spreading of the COVID-19 in Saudi Arabia.
In addition, some scholars used the models with several compartments to put forward analysis and opinions on epidemic prevention and control. In Ref. [22], Carli et al. proposed a multi-region SIRQTHE model and an optimal control approach, which supported governments in defining the most effective strategies to be adopted during post-lockdown mitigation phases in a multi-region scenario. Giordano et al. [23] established a SIDARTHE model that predicted the course of the epidemic to help plan an effective control strategy. On the basis of the expected utility theory, Odagaki [24] carried out a theoretical framework to find out an optimum strategy for minimizing the maximum number of infeced and for controlling the outbreak of pandemic. In Ref. [25], Saldaña et al. developed a compartmental epidemic model of the COVID-19 epidemic outbreak to evaluate the theoretical impact of plausible control interventions.
In this paper, we establish the COVID-19 epidemic model based on the SIQR model. However, these large deviations between model predictions and the actually recorded numbers stem from the uncertainty of the underlying SIQR model parameters: they may not be sufficiently well known for the novel COVID-19 pandemic yet. Generally, considering the latent period gives rise to models with the incorporation of delays to solve problems. Therefore, time-delay has important biologic meaning in epidemic models. Differently, by describing the time-delay from the infection period to the recovery period, we propose a new delay SIQR epidemic model with horizontal transmission in our paper, and obtain the critical treatment time through calculation and stability analysis, which can contribute greatly to medicine.
The rest of the paper is organized as follows: In Section 2, considering the treatment time of epidemic, we present a delayed differential equation for COVID-19. In Section 3, we analyze the existence and stability of equilibria and the existence Hopf bifurcation of the system. In Section 4, we derive the normal form of Hopf bifurcation for above model, then determine the direction of Hopf bifurcation and the stability of bifurcating periodic solutions. In Section 5, we present the simulated results to verify the correctness of the theoretical analysis, and show the effect of treatment time on epidemic. Finally, some conclusions are shown in Section 6.
COVID-19 is similar to other infectious diseases. The susceptible population can be infected by COVID-19 carriers. Some of the infected people are quarantined, and some are killed by COVID-19. Both infected and quarantined people may be cured after treatment. A feature of COVID-19 can be obtained by combining a large number of cases, that is, cured infected people are difficult to be infected again. In Ref. [26], Liu et al. proposed an infectious disease model as follows:
{dSdt=Λ−μS−βSIN,dIdt=βSIN−(μ+γ+δ+α)I,dQdt=δI−(μ+ϵ+α)Q,dRdt=γI+ϵQ−μR, | (2.1) |
where S, I, Q, R denote the numbers of susceptible, infective, quarantined and removed, N=S+I+Q+R is the number of total population individuals. The parameter Λ is the recruitment rate of S corresponding to births and immigration; β denotes the average number of adequate contacts; μ is the natural death rate; γ and ϵ denote the recover rates from group I, Q to R, respectively; δ denotes the removal rate from I; α is the disease-caused death rate of I and Q. The parameters involved in system (2.1) are all positive constants.
On the basis of reference [26], we define β as the contact rate, which eliminates the calculation of the total population number N and simplifies the form of the equation. It is a good way to reduce the error of the model results by reducing the parameters. Combined with practical reports, it is known that some quarantined patients will receive treatment, and the mortality of treated patients will be greatly reduced. Therefore, unlike the model in Ref.[26], we believe that the COVID-19 mortality of infected and quarantined person is different. Meanwhile, the uncertainty of the model parameters is also the main reason for the difference between the predicted results and the actual. Here, we consider the time-delay (τ) from infection to recovery process on the basis of the model in Ref.[26], and the critical treatment time obtained by the calculation with time-delay is helpful to the treatment of COVID-19, which is rarely studied before and should be paid attention to in severe epidemic areas. In particular, the model can change the parameters according to the situation of different regions to get the corresponding critical treatment time, and take reasonable measures to control the epidemic effectively. Finally, we present a new COVID-19 epidemic model, and the relationships between the four populations (susceptible population(S), infected population(I), quarantined population(Q), recovered population(R)) are obtained, as shown in Figure 1.
Thus, we construct the following COVID-19 epidemic model:
{˙S=Λ−μS−βSI,˙I=βSI−δI−μI−α1I−γI(t−τ),˙Q=δI−ϵQ−α2Q−μQ,˙R=γI(t−τ)+ϵQ−μR, | (2.2) |
where Λ,μ,β,δ,γ,ϵ,α1,α2 are parameters; S,I,Q,R are control variables and τ is the time-delay. The specific definitions are given in the Table 1.
Symbol | Definition |
S | Number of susceptible people |
I | Number of infected people |
Q | Number of quarantined people |
R | Number of recovered people |
Λ | Natural increase of population |
β | Transition rate from S to I |
δ | Transition rate from I to Q |
γ | Transition rate from I to R, the cure rate of infected persons |
ϵ | Transition rate from Q to R, the cure rate of quarantined persons |
α1 | COVID-19 mortality rate of infected persons |
α2 | COVID-19 mortality rate of quarantined persons |
μ | Natural death of population |
τ | The time-delay from infection to recovery process |
In this section, system (2.2) is considered. Obviously, system (2.2) has two equilibria:
E1=(S∗1,I∗1,Q∗1,R∗1),E2=(S∗2,I∗2,Q∗2,R∗2), |
where
S∗1=Λμ, I∗1=0, Q∗1=0, R∗1=0, S∗2 = γ+α1+δ+μβ, I∗2 = Λβ−μ(γ+α1+δ+μ)β(γ+α1+δ+μ),
Q∗2 = δ[Λβ−μ(γ+α1+δ+μ)]β(γ+α1+δ+μ)(α2+ϵ+μ), R∗2 = [Λβ−μ(γ+α1+δ+μ)][γ(ϵ+α2+μ)+ϵδ]βμ(γ+α1+δ+μ)(α2+ϵ+μ).
Firstly, we consider the first equilibrium E1=(S∗1,I∗1,Q∗1,R∗1)=(Λμ,0,0,0). Transferring the equilibrium to the origin and linearizing the system around it, we obtain the characteristic equation of the linearized system as follows:
(λ+μ)2(λ+ϵ+α2+μ)(λ−βS∗1+δ+α1+μ+γe−λτ)=0, | (3.1) |
where S∗1 = Λμ.
When τ=0, Eq. (3.1) becomes
(λ+μ)2(λ+ϵ+α2+μ)(λ−βS∗1+δ+α1+μ+γ)=0. | (3.2) |
Eq. (3.2) has four roots: λ1=λ2=−μ,λ3=−ϵ−α2−μ,λ4=βS∗1−δ−α1−μ−γ, due to μ>0,ϵ>0,α2>0, so in the actual situation, λ1<0, λ2<0, λ3<0.
We consider the following assumption:
(H1) βΛμ<μ+δ+α1+γ.
When (H1) holds, all the roots of Eq. (3.2) have negative real parts, and the equilibrium E1 is locally asymptotically stable when τ=0.
When τ>0, let λ=iω(ω>0) be a root of Eq. (3.1). Due to λ1<0, λ2<0, λ3<0, actually, we only need to discuss the following equation:
λ−βS∗1+δ+α1+μ+γe−λτ=0. | (3.3) |
Substituting λ=iω(ω>0) into Eq. (3.3) and separating the real and imaginary parts, we obtain
{βS∗1−(δ+α1+μ)=γcos(ωτ),ω=γsin(ωτ). | (3.4) |
Eq. (3.4) leads to
{cos(ωτ)=βS∗1−(δ+α1+μ)γ,sin(ωτ)=ωγ. | (3.5) |
Adding the square of two equations of Eq. (3.5), we have
h(ω)=ω2+(βS∗1−δ−α1−μ)2−γ2=0. | (3.6) |
Therefore, we give the following assumption:
(H2) (βS∗1−δ−α1−μ)2−γ2<0.
If (H2) holds, then Eq. (3.6) has one positive root ω0=√γ2−(βS∗1−δ−α1−μ)2. Substituting ω0 into Eq. (3.5), we get
τ(j)1={1ω0[arccos(P0)+2jπ],Q0≥0,1ω0[2π−arccos(P0)+2jπ],Q0<0,j=0,1,2,⋯, | (3.7) |
where
Q0=sin(ω0τ(j)1)=ω0γ,P0=cos(ω0τ(j)1)=βS∗1−(δ+α1+μ)γ. |
Lemma 3.1. If (H2) holds, when τ=τ(j)1(j=0,1,2,⋯), then Eq. (3.1) has a pair of pure imaginary roots ±iω0, and all the other roots of Eq. (3.1) have nonzero real parts.
Furthermore, let λ(τ)=α(τ)+iω(τ) be the root of Eq. (3.1) satisfying α(τ(j)1)=0, ω(τ(j)1)=ω0(j=0,1,2,⋯).
Lemma 3.2. If (H2) holds, we have the following transversality conclusions:
Re(dτdλ)|τ=τ(j)1=Re(dλdτ)−1|τ=τ(j)1=1γ2>0, where j=0,1,2,⋯.
Secondly, for the other equilibrium E2 = (S∗2,I∗2,Q∗2,R∗2) of the system (2.2), similarly, transferring the equilibrium to the origin and linearizing the system around it, we obtain the characteristic equation of the linearized system as follows:
(λ+μ)(λ+ϵ+α2+μ)[(λ+μ+βI∗2)(λ−γ+γe−λτ)+β2I∗2S∗2]=0, | (3.8) |
where S∗2 = γ+α1+δ+μβ, I∗2 = Λβ−μ(γ+α1+δ+μ)β(γ+α1+δ+μ).
When τ=0, Eq. (3.8) becomes
(λ+μ)(λ+ϵ+α2+μ)[λ2+(μ+βI∗2)λ+β2I∗2S∗2]=0. | (3.9) |
Eq. (3.9) has four roots: λ1=−μ, λ2=−ϵ−α2−μ, λ3=−μ−βS∗2+√(μ+βS∗2)2−4β2S∗2I∗224, λ4=−μ−βS∗2−√(μ+βS∗2)2−4β2S∗2I∗22. Due to μ>0, ϵ>0, α2>0, in the actual situation, λ1<0, λ2<0.
Note that β2S∗2I∗2>0, we consider the following assumption:
(H3) βΛμ>μ+δ+α1+γ.
Therefore, under the assumption (H3), all the roots of Eq. (3.9) have negative real parts, and the equilibrium E2=(S∗2,I∗2,Q∗2,R∗2) is locally asymptotically stable when τ=0.
When τ>0, let λ=iω(ω>0) be a root of Eq. (3.8). Due to μ>0, ϵ>0, α2>0, actually, we only need to consider the equation:
(λ+μ+βI∗2)(λ−γ+γe−λτ)+β2I∗2S∗2=0. | (3.10) |
Substituting λ=iω(ω>0) into Eq. (3.10) and separating the real and imaginary parts, we have
{ω2+γμ+γβI∗2−β2S∗2I∗2=(γμ+γβI∗2)cosωτ+γωsinωτ,ω(μ−γ+βI∗2)=−γωcosωτ+(γμ+γβI∗2)sinωτ. | (3.11) |
Eq. (3.11) leads to
{cosωτ=(ω2+γμ+γβI∗2−β2S∗2I∗2)(γμ+γβI∗2)−γω2(μ−γ+βI∗2)(γμ+γβI∗2)2+γ2ω2,sinωτ=γω(μ+βI∗2)(μ−γ+βI∗2)+γω(ω2+γμ+γβI∗2−β2S∗2I∗2)(γμ+γβI∗2)2+γ2ω2. | (3.12) |
Adding the square of two equations of Eq. (3.12), let z=ω2, then
h(z)=z2+c1z+c0=0, | (3.13) |
where c1=(μ+βI∗2)2−2β2S∗2I∗2, c0=(γμ+γβI∗2−β2S∗2I∗2)2−(γμ+γβI∗2)2. Therefore, we give the following assumptions:
(H4) c0<0.
(H5) c21−4c0>0, c1<0,c0>0.
If (H4) holds, then Eq. (3.13) has only one positive real root z1. If (H5) holds, then Eq. (3.13) has two positive real roots z2 and z3. Substituting ωk=√zk(k=1,2,3) into Eq. (3.12), we get
τ(j)2,k={1ωk[arccos(Pk)+2jπ],Qk≥0,1ωk[2π−arccos(Pk)+2jπ],Qk<0,k=1,2,3,j=0,1,2,⋯, | (3.14) |
where
Qk=sin(ωkτ(j)2,k)=γωk(μ+βI∗2)(μ−γ+βI∗2)+γωk(ω2k+γμ+γβI∗2−β2S∗2I∗2)(γμ+γβI∗2)2+γ2ω2k,Pk=cos(ωkτ(j)2,k)=(ω2k+γμ+γβI∗2−β2S∗2I∗2)(γμ+γβI∗2)−γω2k(μ−γ+βI∗2)(γμ+γβI∗2)2+γ2ω2k. |
Lemma 3.3. If (H4) or (H5) holds, when τ=τ(j)2,k(k=1,2,3;j=0,1,2,⋯), then Eq. (3.8) has a pair of pure imaginary roots ±iωk, and all the other roots of Eq. (3.8) have nonzero real parts.
Furthermore, let λ(τ)=α(τ)+iω(τ) be the root of Eq. (3.8) satisfying α(τ(j)2,k)=0, ω(τ(j)2,k)=ωk(k=1,2,3;j=0,1,2,⋯).
Lemma 3.4. If (H4) or (H5) holds, and zk=ω2k, h′(zk)≠0, then we have the following transversality conclusions:
Re(dτdλ)|τ=τ(j)2,k=Re(dλdτ)−1|τ=τ(j)2,k=h′(zk)γ2[(μ+βI∗2)2+ω2k]≠0, k=1,2,3,j=0,1,2,⋯.
Theorem 3.1. We show the conclusion associated with two equilibria of the system (2.2).
(1) If the assumptions (H1) and (H2) hold, the equilibrium E1 of the system (2.2) undergoes Hopf bifurcation at τ=τ(j)1(j=0,1,2,⋯), where τ(j)1 is given by Eq. (3.7), and we have: when τ∈[0,τ(0)1), the equilibrium E1 is locally asymptotically stable, and the equilibrium E1 is unstable when τ>τ(0)1.
(2) If the assumptions (H4) or (H5) holds, the equilibrium E2 of the system (2.2) undergoes Hopf bifurcation at τ=τ(j)2,k(k=1,2,3;j=0,1,2,⋯), where τ(j)2,k is given by Eq. (3.14), and
(a) If the assumptions (H3) and (H4) hold, h(z) has one positive root z1, then when τ∈[0,τ(0)2,1), the equilibrium E2 is locally asymptotically stable, and the equilibrium E2 is unstable when τ>τ(0)2,1.
(b) If the assumptions (H3) and (H5) hold, h(z) has two positive roots z2 and z3, we suppose z2<z3, then h′(z2)<0,h′(z3)>0, note that τ(0)2,2>τ(0)2,3. Then ∃ m∈N makes 0<τ(0)2,3<τ(0)2,2<τ(1)2,3<τ(1)2,2<⋯<τ(m−1)2,2<τ(m)2,3<τ(m+1)2,3. When τ∈[0,τ(0)2,3)∪m⋃l=1(τ(l−1)2,2,τ(l)2,3), the equilibrium E2 of the system (2.2) is locally asymptotically stable, and when τ∈m−1⋃l=0(τ(l)2,3,τ(l)2,2)∪(τ(m)2,3,+∞), the equilibrium E2 of the system (2.2) is unstable.
In this section, we discuss the normal form of Hopf bifurcation for the system (2.2) by using the multiple time scales method. Combining with actual situation, we concern about the impact of treatment time on epidemic control. Therefore, we consider the time-delay τ as a bifurcation parameter, let τ=τc+ετε, where τc is the critical value of Hopf bifurcation given in Eq. (3.7) or Eq. (3.14) respectively, τε is the disturbance parameter, and ε is the dimensionless scale parameter. We suppose the characteristic Eq. (3.1) and Eq. (3.8) have eigenvalue λ=iω(k)(k=1,2), where ω(1)=ω0, ω(2)=ω1, ω2 or ω3, at which system (2.2) undergoes a Hopf bifurcation at equilibrium Ek=(S∗k,I∗k,Q∗k,R∗k), k=1,2, respectively.
Then system (2.2) can be written as
˙X(t)=AX(t)+BX(t−τ)+F[X(t),X(t−τ)], | (4.1) |
where X(t)=(Sk,Ik,Qk,Rk)T, X(t−τ)=(Sk(t−τ),Ik(t−τ),Qk(t−τ),Rk(t−τ))T,
A=(−μ−βI∗k−βS∗k00βI∗kβS∗k−(δ+μ+α1)000δ−(ϵ+α2+μ)000−ϵ−μ),B=(00000−γ0000000−γ00), |
F(X(t),X(t−τ))=(FSFIFQFR)=(Λ−μS∗k−βS∗kI∗k−βSkIkβSkIk+βS∗kI∗k−(δ+μ+α1+γ)I∗kδI∗k−(ϵ+α2+μ)Q∗kγI∗k−ϵQ∗k−μR∗k). |
We suppose hk, h∗k(k=1,2) are the eigenvector of the corresponding eigenvalue λ=iω(k), λ=−iω(k) respectively of Eq. (4.1) for equilibrium Ek, and satisfies ⟨h∗k,hk⟩=¯h∗kThk=1. By simple calculation, we can get:
hk=(hk1,hk2,hk3,hk4)T=(1,−λ+μ+βI∗kβS∗k,−δ(λ+μ+βI∗k)(λ+ϵ+α2+μ)βS∗k,δϵ(λ+μ+βI∗k)−γ(λ+μ+βI∗k)(λ+ϵ+α2+μ)e−λτ(λ+μ)(λ+ϵ+α2+μ)βS∗k)T,h∗k=(h∗k1,h∗k2,h∗k3,h∗k4)T=dk(1,−λ+μ+βI∗kβI∗k,0,0)T, | (4.2) |
where dk=β2S∗kI∗kβ2S∗kI∗k−(−λ+μ+βI∗k)2,k=1,2.
We suppose the solution of Eq. (4.1) as follows:
X(t)=X(T0,T1,T2,⋯)=∞∑k=1εkXk(T0,T1,T2,⋯), | (4.3) |
where X(T0,T1,T2,⋯)=[S(T0,T1,T2,⋯),I(T0,T1,T2,⋯),Q(T0,T1,T2,⋯),R(T0,T1,T2,⋯)]T,Xk(T0,T1,T2,⋯)=[Sk(T0,T1,T2,⋯),Ik(T0,T1,T2,⋯),Qk(T0,T1,T2,⋯),Rk(T0,T1,T2,⋯)]T.
The derivative with respect to t is transformed into:
ddt=∂∂T0+ε∂∂T1+ε2∂∂T2+⋯=D0+εD1+ε2D2+⋯, |
where Di=∂∂Ti, i=0,1,2,⋯.
We make Xj=(Sj,Ij,Qj,Rj)T=Xj(T0,T1,T2,⋯), Xj,τc=(Sj,τc,Ij,τc,Qj,τc,Rj,τc)T=Xj(T0−τc,T1,T2,⋯), j=1,2,3,⋯.
From Eq. (4.3) we can get:
˙X(t)=εD0X1+ε2D1X1+ε3D2X1+ε2D0X2+ε3D0X3+⋯. | (4.4) |
The Taylor expansion of X(t−τ) is carried out:
X(t−τ)=εX1,τc+ε2X2,τc+ε3X3,τc−ε2τεD0X1,τc−ε3τεD0X2,τc−ε2τcD1X1,τc−ε3τεD1X1,τc−ε3τcD2X1,τc−ε3τcD1X2,τc+⋯, | (4.5) |
where Xj,τc=Xj(T0−τc,T1,T2,⋯),j=1,2,3,⋯.
Substituting Eqs. (4.3) ∼ (4.5) into Eq. (4.1), and balancing the coefficients before ε on both sides of the equation, the following expression is obtained:
D0Sk1+μSk1+βI∗kSk1+βS∗kIk1=0,D0Ik1−βI∗kSk1−βS∗kIk1+(δ+μ+α1)Ik1+γIk1,τc=0,D0Qk1−δIk1+(ϵ+α2+μ)Qk1=0,D0Rk1−γIk1,τc+ϵQk1+μRk1=0,k=1,2. | (4.6) |
Thus Eq. (4.6) has the following solution form:
X1(T1,T2,T3,⋯)=G(T1,T2,T3,⋯)eiω(k)T0hk+ˉG(T1,T2,T3,⋯)e−iω(k)T0ˉhk,k=1,2. | (4.7) |
The expression of the coefficient before ε2 is as follows:
D0Sk2+μSk2+βI∗kSk2+βS∗kIk2=−D1Sk1−βSk1Ik1,D0Ik2−βI∗kSk2−βS∗kIk2+(δ+μ+α1)Ik2+γIk2,τc=−D1Ik1+βSk1Ik1+γ(τεD0Ik1,τc+τcD1Ik1,τc),D0Qk2−δIk2+(ϵ+α2+μ)Qk2=−D1Qk1,D0Rk2−γIk2,τc+ϵQk2+μRk2=−D1Rk1−γ(τεD0Ik1,τc+τcD1Ik1,τc),k=1,2. | (4.8) |
Substituting Eq. (4.7) into the right hand side of Eq. (4.8), and the coefficient vector of eiω(k)T0 is denoted by m1. According to the solvable condition ⟨h∗k,m1⟩=0, the expression of ∂G∂T1 can be obtained as follows:
∂G∂T1=NkτεG, | (4.9) |
where Nk=(iω(k)+μ+βI∗k)2γ⋅iω(k)e−iω(k)τc(iω(k)+μ+βI∗k)2−β2S∗kI∗k−(iω(k)+μ+βI∗k)2γτce−iω(k)τc,k=1,2.
Since τε is a disturbance parameter, we only consider its effect on the linear part. It has little effect on the high order, so it can be ignored. Therefore, we ignore the part containing τε in the higher order. We suppose:
Sk2=gk1e2iω(k)T0G2+ˉgk1e−2iω(k)T0ˉG2+lk1GˉG,Ik2=gk2e2iω(k)T0G2+ˉgk2e−2iω(k)T0ˉG2+lk2GˉG,Qk2=gk3e2iω(k)T0G2+ˉgk3e−2iω(k)T0ˉG2+lk3GˉG,Rk2=gk4e2iω(k)T0G2+ˉgk4e−2iω(k)T0ˉG2+lk4GˉG, | (4.10) |
where
gk1=−βhk2(2iω(k)+δ+μ+α1+γe−2iω(k)τc)J,gk2=(2iω(k)+μ)βhk2J,gk3=(2iω(k)+μ)δβhk2(2iω(k)+μ+α2+ϵ)J,gk4=(2iω(k)+μ)βhk2[ϵ(2iω(k)+μ+α2+ϵ)−δϵ](2iω(k)+μ+α2+ϵ)(2iω(k)+μ)J,lk1=−β(hk2+ˉhk2)(δ+μ+α1+γ)V,lk2=βμ(hk2+ˉhk2)V,lk3=δβμ(hk2+ˉhk2)(μ+α2+ϵ)V,lk4=γβμ(hk2+ˉhk2)(μ+α2+ϵ)−ϵδβμ(hk2+ˉhk2)μ(μ+α2+ϵ)V,J=(2iω(k)+μ+βI∗k)(2iω(k)+δ+μ+α1+γe−2iω(k)τc)−βS∗k(2iω(k)+μ),V=(μ+βI∗k)(δ+μ+α1+γ)+βμS∗k,hk2=−iω(k)+μ+βI∗kβS∗k,k=1,2. | (4.11) |
The expression of the coefficient before ε3 is:
D0Sk3+μSk3+βI∗kSk3+βS∗kIk3=−D2Sk1−D1Sk2−βSk2Ik1−βSk1Ik2,D0Ik3−βI∗kSk3−βS∗kIk3+(δ+μ+α1)Ik3+γIk3,τc=−D2Ik1−D1Ik2+βSk2Ik1+βSk1Ik2+γ(τcD1Ik2,τc+τcD2Ik1,τc),D0Qk3−δIk3+(ϵ+α2+μ)Qk3=−D2Qk1−D1Qk2,D0Rk3−γIk3,τc+ϵQk3+μRk3=−D2Rk1−D1Rk2−γ(τcD1Ik2,τc+τcD2Ik1,τc),k=1,2. | (4.12) |
Substituting Eq. (4.7), Eq. (4.10) and Eq. (4.11) into the right hand side of Eq. (4.12), and the coefficient vector of eiω(k)T0 is denoted by m2. According to the solvable condition ⟨h∗k,m2⟩=0, the expression of ∂G∂T2 can be obtained as follows:
∂G∂T2=ξkG2ˉG, | (4.13) |
where
ξk=β(lk2+gk2+lk1hk2+gk1ˉhk2)(iω(k)+μ)βI∗k−(iω(k)+μ+βI∗k)(γτce−iω(k)τc−1)hk2,k=1,2. |
Let G→G/ε, then the deduce to third-order normal form of Hopf bifurcation of system (2.2) is:
˙G=NkτεG+ξkG2ˉG, | (4.14) |
where Nk is given in Eq. (4.9), and ξk is given in Eq. (4.13), k=1,2.
By substituting G=reiθ into Eq. (4.14), the following normal form of Hopf bifurcation in polar coordinates can be obtained:
{˙r=Re(Nk)τεr+Re(ξk)r3,˙θ=Im(Nk)τε+Im(ξk)r2, | (4.15) |
where Nk is given in Eq. (4.9), and ξk is given in Eq. (4.13), k=1,2.
According to the normal form of Hopf bifurcation in polar coordinates, we just need to consider the first equation in system (4.15). Thus, there is the following theorem:
Theorem 4.1. For the system (4.15), when Re(Nk)τεRe(ξk)<0(k=1,2), there is a nontrivial fixed point r=√−Re(Nk)τεRe(ξk), and system (2.2) has periodic solution:
(1) If Re(Nk)τε<0, then the periodic solution reduced on the center manifold is unstable.
(2) If Re(Nk)τε>0, then the periodic solution reduced on the center manifold is stable.
In this section, according to the data presented in Refs. [27] and [28], we let Λ=10.48, β=0.0004, δ=0.7, α1=0.0484, α2=0.022778, γ=0.42386, ϵ=0.475, μ=0.00714. Then, system (2.2) only has one nonnegative equilibrium E1=(S∗1,I∗1,Q∗1,R∗1)≈(1467.8,0,0,0). Obviously, the assumption (H1) holds, thus the equilibrium E1 is locally asymptotically stable when τ = 0.
Substituding these parameter values into the Eqs. (3.5)∼(3.7), using MATLAB, we can obtain ω0≈0.3900, Q0≈0.9201, P0≈−0.4284, τ(0)1≈5.1629. According to the Theorem 3.1, the equilibrium E1 is locally asymptotically stable at τ∈[0,τ(0)1), and when τ=τ(0)1, Hopf bifurcation occurs near the equilibrium E1. Then, we obtain Re(N1)>0, Re(ξ1)>0 from Eqs. (4.9)∼(4.13), thus according to the Theorem 4.1, the system (2.2) has backward periodic solution and the periodic solution is unstable when τε<0.
When τ=0, we choose the initial value (1473,1,2,1) and the corresponding locally asymptotically stable equilibrium is shown in Figure 2.
When τ=4.2∈[0,τ(0)1), we choose the initial value (1473,1,2,1), the equilibrium E1 is also locally asymptotically stable (see Figure 3). Thus, the epidemic will be controlled eventually.
When τ=5.5∈(τ(0)1,+∞), we choose the initial value (1473,1,2,1), and the equilibrium E1 is unstable shown in Figure 4.
It can be seen from Figure 2∼ Figure 4, when τ∈[0,τ(0)1), the equilibrium E1 of system (2.2) is locally asymptotically stable (see Fig. 2 and Fig. 3). When τ∈(τ(0)1,+∞), the equilibrium E1 of system (2.2) is unstable, and accompanied by the fluctuation with increased amplitude (see Fig. 4). There is no periodic solution, so we believe the theoretical analysis is correct.
Remark 1: Through numerical simulations, it can be found that when the treatment time τ<τ(0)1, the shorter treatment time is, the faster the epidemic tends to be stable. If τ>τ(0)1, the epidemic cannot be controlled effectively, and accompanied by the fluctuation with increased amplitude with time. Thus, we obtain the critical treatment time τ(0)1 for the epidemic.
According to the data presented in Refs. [27] and [28], we choose another set of parameters, Λ=2.48, β=0.86834, δ=0.3, α1=0.0484, α2=0.022778, γ=0.42386, ϵ=0.475, μ=0.00714. We obtain E1=(S∗1,I∗1,Q∗1,R∗1)≈(347.3,0,0,0), E2=(S∗2,I∗2,Q∗2,R∗2)≈(0.8976,3.1737,1.8857,313.8526). Apparently, the assumption (H1) does not hold, but the assumption (H3) holds, therefore, the equilibrium E1 is unstable and the equilibrium E2 is locally asymptotically stable when τ=0.
Substituding these parameter values into the Eq. (3.13), we get c1≈3.3383, c0≈−0.4174. It satisfies the assumption (H4), therefore, Eq. (3.13) has only one positive root z1≈0.1204. Substituting these parameter values into the Eqs. (3.12)∼(3.14), using MATLAB, we can obtain ω1≈0.3474, Q1≈0.5926, P1≈−0.8055, τ(0)2,1≈7.2175. According to the Theorem 3.1, the equilibrium E2 is locally asymptotically stable at τ∈[0,τ(0)2,1), and when τ=τ(0)2,1, Hopf bifurcation occurs near the equilibrium E2. Then we obtain Re(N2)>0, Re(ξ2)<0 by Eqs. (4.9)∼(4.13), according to the Theorem 4.1, the system (2.2) has stable periodic solution near E2, and the periodic solution reduced on the center manifold is stable.
When τ = 0, we choose the initial value (0.85,3,1.9,300) and the corresponding locally asymptotically stable equilibrium is shown in Figure 5.
When τ=4∈[0,τ(0)2,1), we choose the initial value (0.85,3,1.9,300), and the periodic solution of this model is locally asymptotically stable shown in Figure 6.
When τ=5∈[0,τ(0)2,1), we choose the initial value (0.85,3,1.9,300), the equilibrium E2 has a stable equilibrium is shown in Figure 7.
It can be seen from Figure 5∼Figure 7, when τ∈[0,τ(0)2,1), the equilibrium E2 of system (2.2) is locally asymptotically stable. Especially, if τ is larger, the time required for the system (2.2) to be controlled stable is longer.
When τ=7.23>τ(0)2,1=7.2175, we choose the initial value (0.85,3,1.9,300), the model has stable periodic solution shown in Figure 8.
When τ=7.5>7.23, we choose the initial value (0.85,3,1.9,300), and the equilibrium E2 is unstable shown in Figure 9.
When τ=7.9>7.5, we choose the initial value (0.85,3,1.9,300), and the equilibrium E2 is unstable shown in Figure 10.
It can be seen from the Figure 8∼Figure 10, when τ∈(τ(0)2,1,+∞), the equilibrium E2 of system (2.2) is unstable. When τ approaches the critical value τ(0)2,1, the equilibrium E2 of system (2.2) exhibits periodic fluctuation and bifurcates stable periodic solutions (see Fig. 8). As τ increasing, the fluctuation phenomenon lasts for a period of time, and shows increasing volatility trends (see Fig. 9) and disappears eventually (see Fig. 10). According to the Theorem 3.1 and Theorem 4.1, we know that the theoretical analysis is correct.
Remark 2: Through numerical simulations, it can be found that when the treatment time τ<τ(0)2,1, the epidemic can be controlled effectively, and the smaller the time-delay τ is the faster the epidemic is controlled. When the treatment time τ∈(τ(0)2,1,+∞) is close to τ(0)2,1, the epidemic will repeatedly occur periodically, otherwise, the epidemic cannot be controlled. Since the small variation in treatment time will lead to the epidemic polarization, therefore, whether to grasp the treatment time is essential to control the epidemic.
Remark 3: In Ref. [26], Liu et al. used the Milstein's Higher Order Method mentioned to illustrate their main results. And they obtained that the disease is persistent. Similarly, in our simulation, the epidemic is persistent and tends to stable with time.
In Ref. [29], the study carried out numerical simulation. First of all, they changed the value of parameters to get the "without control'' and "with control'' graphs, from which we can see that different populations tend to be stable in the end. At the same time, through comparison, it is concluded that the control technique is useful in controlling the disease. Our simulation results show that the system is locally asymptotically stable in the critical treatment time, and the disease will be effectively controlled. Then they fit real data with the infected class of their model and found that the number of cases grows exponentially with time, the disease would be serious without applying the proper optimal control strategies. And we conclude that if not treat within the critical time, that is to say, without timely and effective treatment, the epidemic will be uncontrollable.
In a word, although the model studied is different from Ref. [29], and we focus on the impact of time delay on the epidemic situation, but the overall idea and conclusion are consistent, so we can still verify the correctness of our simulation results. Moreover, we emphasize the time-delay of treatment, so we also obtain the critical treatment time, which can take more targeted measures to control the epidemic. To effectively control the epidemic situation, we need to take measures not only in the external environment, but also in the internal medical care, and now there is little research on the treatment direction. Therefore, the critical treatment time we get is helpful for the alleviation of the disease.
In this paper, according to the propagation characteristics of COVID-19, we have constructed the SIRQ epidemic model with time delay for the COVID-19. We have also analyzed the stability of the equilibria and the existence of Hopf bifurcation associated with both equilibria. Then, we have used the multiple time scales method to derive the normal form of Hopf bifurcation for above COVID-19 epidemic model. Finally, We have chosen two groups of parameter values according to the data presented in Refs. [27] and [28] for numerical simulations to verify the correctness of the theoretical analysis. Compared with the numerical simulation results of Refs. [26] and [29], we have got the conclusion that our numerical simulation results were accuracy.
The numerical simulations showed that the small change of time-delay τ leads to the epidemic from stable to uncontrollable: the smaller τ is the better the epidemic controlled, but with τ increasing, the epidemic will occur repeatedly and outbreak eventually. Therefore, it is significant that find the critical treatment time to control the epidemic. In addition, we would predict the critical treatment time of epidemic in different regions according to the relevant parameters, so as to effectively control the epidemic in treatment, and realize the major breakthrough of medicine in the epidemic situation.
As for the part of numerical simulation, according to the opinions of reviewers, we will use more accurate data to get more valuable conclusions in further research.
The authors are extremely grateful to the anonymous referees and editor for their careful reading, valuable comments and helpful suggestions, which have helped us to improve the presentation of this work significantly. This study was funded by Fundamental Research Funds for the Central Universities of China (Grant No. 2572019BC14), the Heilongjiang Provincial Natural Science Foundation of China (Grant No. LH2019A001) and College Students Innovations Special Project funded by Northeast Forestry University of China (No.202010225035).
The Authors declare that this work has no conflict of interest.
[1] | G. He, Y. Pan, T. Tanaka, The short-term impacts of COVID-19 lockdown on urban air pollution in China, Nat. Sustain., 3 (2020). |
[2] |
X. Wu, R. C. Nethery, B. M. Sabath, D. Braun, F. Dominici, Air pollution and COVID-19 mortality in the United States: Strengths and limitations of an ecological regression analysis, Sci. Adv., 6 (2020), 1005–1011. doi: 10.1126/sciadv.abb1005
![]() |
[3] | Y. Han, J. C. Lam, V. O. Li, P. Guo, Q. Zhang, A. Wang, et al., The effects of outdoor air pollution concentrations and lockdowns on Covid-19 infections in Wuhan and other provincial capitals in China, preprints, (2020), 2020030364. Available from: https://www.preprints.org/manuscript/202003.0364/v1 |
[4] |
G. Albano, V. Giorno, Inferring time non-homogeneous Ornstein Uhlenbeck type stochastic process, Comput. Stat. Data Anal., 150 (2020), 107008–107008. doi: 10.1016/j.csda.2020.107008
![]() |
[5] |
G. Albano, V. Giorno, On short-term loan interest rate models: a first passage time approach, Mathematics, 6 (2018), 70. doi: 10.3390/math6050070
![]() |
[6] | V. Linetsky, Computing hitting time densities for CIR and OU diffusions: Applications to mean reverting models, J. Comput. Financ., 4 (2004), 1–22. |
[7] |
S. Ditlevsen, P. Lánský, Estimation of the input parameters in the Ornstein-Uhlenbeck neuronal model, Phys. Rev. E, 71 (2005), 011907. doi: 10.1103/PhysRevE.71.011907
![]() |
[8] |
H. Tuckwell, F. Wan, J. P. Rospars, A spatial stochastic neuronal model with Ornstein–Uhlenbeck input current, Biol. Cybern., 86 (2002), 137–145. doi: 10.1007/s004220100283
![]() |
[9] | G. Albano, V. Giorno, P. Román-Pomán, F. Torres-Ruiz, On a non-homogeneous Gompertz-type diffusion process: inference and first passage time, in LNCS 10672 (eds. R. Moreno-Díaz et al.), Springer (2018), 47–54. |
[10] |
A. Buonocore, L. Caputo, A. G. Nobile, E. Pirozzi, Restricted ornstein uhlenbeck process and applications in neuronal models with periodic input signals, J. Comput. Appl. Math., 285 (2015), 59–71. doi: 10.1016/j.cam.2015.01.042
![]() |
[11] |
V. Giorno, S. Spina, On the return process with refractoriness for non-homogeneous Ornstein-Uhlenbeck neuronal model, Math. Biosci. En., 11 (2014), 285–302. doi: 10.3934/mbe.2014.11.285
![]() |
[12] |
R. Gutiérrez, R. Gutiérrez-Sánchez, A. Nafidi, A. Pascual, Detection, modelling and estimation of non linear trends by using a non-homogeneous Vasicek stochastic diffusion. Application to CO2 emissions in Morocco, Stoch. Environ. Res. Risk Assess., 26 (2012), 533–543. doi: 10.1007/s00477-011-0499-z
![]() |
[13] |
G. Albano, V. Giorno, Inference on the effect of non homogeneous inputs in Ornstein Uhlenbeck neuronal modeling, Math. Biosci. Eng., 17 (2020), 328–348. doi: 10.3934/mbe.2020018
![]() |
[14] |
D. W. K. Andrews, R. C. Fair, Inference in nonlinear econometric models with structural change, Rev. Econ. Stud., 55 (1988), 615–639. doi: 10.2307/2297408
![]() |
[15] | E. Hansen, Approximate asymptotic p values for structural-change tests, J. Bus. Econ. Stat., 15 (1997), 60–67. |
[16] |
M. L. Parrella, G. Albano, M. La Rocca, C. Perna, Reconstructing missing data sequences in multivariate time series: an application to environmental data, Stat. Methods Appl., 28 (2019), 359–383. doi: 10.1007/s10260-018-00435-9
![]() |
[17] |
G. Albano, V. Giorno, P. Román-Pomán, S. Román-Pomán, J. J. Serrano-Pérez, F. Torres-Ruiz, Inference on an heteroscedastic Gompertz tumor growth model, Math. Biosci., 328 (2020), 108428. doi: 10.1016/j.mbs.2020.108428
![]() |
1. | Yong Li, Huixin Li, Xiangwen Chen, Yingying Wang, Lei Chen, Man Yang, Jingguang Tang, 2023, Review of Interruption Control Technologies for Cascading Faults in Power Systems, 979-8-3503-0389-6, 426, 10.1109/ACFPE59335.2023.10455331 | |
2. | Jinlong Fu, Yan Song, Yike Feng, Rumor Spreading Model Considering the Roles of Online Social Networks and Information Overload, 2023, 11, 2169-3536, 123947, 10.1109/ACCESS.2023.3328396 | |
3. | Yanchao Liu, Pengzhou Zhang, Lei Shi, Junpeng Gong, A Survey of Information Dissemination Model, Datasets, and Insight, 2023, 11, 2227-7390, 3707, 10.3390/math11173707 | |
4. | Junnan Jiang, 2024, Topic Prediction and Trend Analysis in Social Media Data Mining, 9798400710148, 165, 10.1145/3702879.3702908 |
Symbol | Definition |
S | Number of susceptible people |
I | Number of infected people |
Q | Number of quarantined people |
R | Number of recovered people |
Λ | Natural increase of population |
β | Transition rate from S to I |
δ | Transition rate from I to Q |
γ | Transition rate from I to R, the cure rate of infected persons |
ϵ | Transition rate from Q to R, the cure rate of quarantined persons |
α1 | COVID-19 mortality rate of infected persons |
α2 | COVID-19 mortality rate of quarantined persons |
μ | Natural death of population |
τ | The time-delay from infection to recovery process |
Symbol | Definition |
S | Number of susceptible people |
I | Number of infected people |
Q | Number of quarantined people |
R | Number of recovered people |
Λ | Natural increase of population |
β | Transition rate from S to I |
δ | Transition rate from I to Q |
γ | Transition rate from I to R, the cure rate of infected persons |
ϵ | Transition rate from Q to R, the cure rate of quarantined persons |
α1 | COVID-19 mortality rate of infected persons |
α2 | COVID-19 mortality rate of quarantined persons |
μ | Natural death of population |
τ | The time-delay from infection to recovery process |