key: cord-1017846-kjud6osg
authors: Wang, Qiubao; Wu, Hao
title: There exists the “smartest” movement rate to control the epidemic rather than “city lockdown”
date: 2022-02-19
journal: Appl Math Model
DOI: 10.1016/j.apm.2022.02.018
sha: 7012ebc72e72daafe010402c2b6a23465acb8bc4
doc_id: 1017846
cord_uid: kjud6osg

The emergency outbreak and spread of coronavirus disease 2019 (COVID-19) has left great damage to individuals over most of the world. Population mobility is the primary reason for the spread of the epidemic. A delayed stochastic epidemic susceptible-infected-recovered (SIR) model with Gaussian white noise is introduced. Compared with traditional models,this model is characterized by time delay, environmental noise and population mobility among municipalities with the convenient transportation network. The stochastic dynamic behavior of the SIR model is analyzed and the existence of the stochastic bifurcation of the system is proved. The effect of time delay and movement rate are investigated. Numerical simulations are performed to support the theoretical results. It is worth mentioning that the movement rate is not as low as possible and appropriate population mobility is conducive to alleviating the epidemic. Through simulation, we demonstrate the existence of the best movement rate named the “smartest” [Formula: see text] , which is helpful to control the epidemic. This model is also useful to prevent other infectious diseases.

(2) The movement rate is not as low as possible, there exists the "smartest"  to control the epidemic more effectively.

(3) The short incubation period is conducive to timely measures to control the epidemic.

(4) We should prevent not only local infections but also external infections.

COVID-19 is the disease caused by a new coronavirus called SARS-CoV-2, where "CO" stands for corona, "VI" for virus and "D" for disease, which between cities has been greatly improved [19, 20] . Besides, it is well known that there has a span for infected people who are not quarantined in time to move from one area to others, namely, the time delay exists in the spread of the epidemic between the two regions.

Time delay has a significant effect on the epidemic model. In the mathematical model with the influence of epidemic controlled by Twitter, Hopf bifurcation occurs when the delay increases in [21] . Liu et al. [22] proposed a periodic SIS epidemic model with the impact of time delay and infection related with transport in the patch environment, then verified the steadystate stability of the model. Here are some epidemic models with time delay on complex networks:the epidemic SIS model with time delay on a scale-free network we have never seen before was given by [23] , and drew a conclusion that the basic reproductive number impacts directly for the extinction and persistence of the epidemic; Kang et al. [24] obtains the global stability of the equilibrium point by studying the delay SIS model on scale-free networks. It would be significant for us to study the important role of time delay in spreading the COVID-19.

It is well known that diffusion mechanisms are more complex than assumed in these models in different regions, subject to various demographic, cultural, and transportation influences. The heterogeneity of the system leads to outbreaks occurring in different times and with diverse modes of transmission in various municipalities. The study of deterministic models is helpful to understand the behavior of epidemics, but the determined results are difficult to describe and include the uncertainty of elements in the real world. Epidemic models driven by white noise with different functional responses are analyzed [25, 26] . Recently, a new type of stochastic epidemic model with time delay has become the focus of research [27, 28] . In nature, the dynamics of epidemic are disturbed by environmental noise to simulate jumps caused by some natural disasters that break the continuity of states behavior.

The goal of our work is to extend the model proposed in [29] to the model with Gaussian white noise perturbation and time delay and take into account the effect of population movement rate in [30] . Different from the existing models, we take time delay and environmental noise into consideration. Furthermore, time delay represents the incubation period of COVID-19. We found that the "smartest" κ was different with the change of time delay. Which can not only alleviate the epidemic situation in the infected areas, but also minimize the number of infections in the two areas. The "smartest" κ would help guide the government to take timely measures to control the epidemic. The harm caused by absolutely "city lockdown" was largely avoided by the "smartest" κ. As is known to all, the harm caused by "city lockdown" involves many aspects, such as society, economy, education, etc. The consequences of long-term absolutely "city lockdown" are unbearable to humans. Finally, numerical simulations are carried out in two regions in China to investigate the impact of movement rate on the spread of the epidemic with the convenient transportation network and find the "smartest" κ.

The article is organized as follows. In section 2, we introduce the model and present the conditions for the Hopf bifurcation to occur. In section 3, we discuss the dynamic behaviors of the Hopf bifurcation and calculate the Itô stochastic differential equation. In section 4, numerical simulations are carried out to support our results and accuracy evaluation of the model with the epidemic data of two areas in China to verify the practicability of the model. Finally, discussions and concluding remarks are made and put forward the prospect of the future work in section 5.

We consider two coupled populations (i and j, i ‰ j ) which are adjacent areas in the geographical and the scale of the population are the same. Introducing the SIR model , a stochastic delayed differential equations(SDDEs) with fractional order. Populations of susceptible, infected and recovered are distribute in the flow network between two adjacent areas, then predict the development of the epidemic.

where S i , I ij and R i (j=1,2) represent the local susceptible, infectious and recovered individuals in area i, respectively(N i " S i`Iij`Ri is local Figure 1 : A schematic diagram of the model. population size). To make sure the local (i and j ) population remains unchanged we assume that the birth rate b is equal to the natural mortality rate, the disease mortality rate α is equal to the supplemental population rate. γ ij , κ and β ij (j=1,2) are recovery, movement, and transmission rates, respectively. τ (τ ą 0) is incubation period, which is a time period. I i1 is the infections without symptoms and I i2 is infected people with symptoms who are immobile (I i1`Ii2 " I i , I i represents the infected people in area i).

Let γ " γ i1

R 0 , the basic reproduction number, is an indicator of the contagiousness or transmissibility of infectious and parasitic agents. It represents the number of new infections estimated to stem from a single case in a population that has never seen the disease before. Besides, it is one of the key values that can predict whether an infectious disease will spread into a population (R 0 ą 1) or die out (R 0 ă 1). According to Ref. [31] , the basic reproduction number is calculated to be R 0 " β γ`αφ`b . When R 0 ą 1, the system (2) exists positive equilibrium, we assume that N i " N j " N and the proportion of the local susceptible, infectious and recovered individuals are same in the two places (i and j). Therefore, we obtain the equilibrium point EpS˚, I˚, R˚q, where

Based on Ref. [29] , we take movement rate,environmental noise, incubation period of COVID-19 , birth rate and so on into consideration. Then according to the transmission mode of the epidemic in Ref. [29] combining with our model (system (2)), a schematic diagram is given in Fig.1 .

Next, we shall consider the disturbance of environmental noise near the equilibrium point

where ξ 1 ptq and ξ 2 ptq represent the Gaussian white noise with power spectral density k 1 and k 2 , ξ 1 ptq is independent of ξ 2 ptq, ε is the scale coefficient which is small to ensure the micro disturbance to the system. Note that β N is small, and it scales to the same order as ε.

This leads to the characteristic equation

F 1 pλq, F 2 pλq are as following:

Due to λ 1 "´b ă 0, λ 2 "´pb`2κq ă 0, the roots of the characteristic equation (5) are the solutions of equations F 1 pλq " 0 and F 2 pλq " 0.

Assume that (H1) γ ą b and k ! 1. where

Under the condition (H1), it can be obtained by calculation F 1 p0q ą 0 and F 2 p0q ą 0, thus 0 is not the solution of F 1 p0q " 0 and F 2 p0q " 0.

Lemma 2.1. Assume that (H1) holds. Then the equilibrium of EpS˚, I˚, R˚q of the system(3) with τ " 0 is asymptotically stable.

Proof When τ " 0, the equation (6) is equivalent to the following equation

Let λ 3 and λ 4 be two roots of the equation (8) . Under the condition (H1), then

As for F 2 pλq, we have

Let λ 5 and λ 6 be two roots of the F 2 pλq| τ "0 " 0. Under the condition (H1), then

These imply that Eq.(8) and (9) have no positive roots, hence all the roots of the characteristic equation(5) have negative real parts.A ssume further (H2) b`αµ`γ ă β ă b`αµ`2γ and b ! 1 When τ ą 0, we have the following conclusion.

Lemma 2.2. If (H1) and (H2) hold, the equation F 2 pλq " 0 has no purely imaginary roots for all τ ą 0.

Proof When τ ą 0, let λ " iω (ω ą 0) be the root of the equation F 2 pλq " 0. Through the work [32] , we can obtain that ω satisfies the following equation:

Let z 1 and z 2 be the roots of (10), then by (H1) and (H2)

This implies when τ ą 0, equation F 2 pλq " 0 has no purely imaginary roots.

Lemma 2.3. When τ ą 0, under the condition (H1) and (H2), the equation F 1 pλq " 0 has a pair of purely imaginary roots˘iω 0 (ω 0 ą 0) for τ " τ j , j " 0, 1, 2¨¨¨.

The method of proof is the same as Lemma 2.2 and the proof of the lemma 2.3 can be found in Appendix A. It demonstrates that the equation F 1 pλq " 0 has a pair of pure imaginary roots. Furthermore, we can solve τ j and ω 0 , as follow

Define τ 0 " mintτ j u, thus˘iω 0 is a pair of simple imaginary roots of F 1 plambdaq " 0 when τ " τ 0 .

Theorem 2.1. Assume (H1) and (H2) hold, the equilibrium point EpS˚, I˚, R˚q of the system (3) occurs stochastic Hopf bifurcation for τ " τ 0 .

Proof Due to the sign of Rep dλ dτ q and Rer dλ dτ s´1 is the same. Differentiating the two sides of the equation F 1 pλq " 0 with respect to τ , assume that (H1) and (H2) hold, we have the Hopfs transversality condition

According to the Hopf bifurcation theory [33] , combining the above results and Lemma 2.2, the equation (3) has periodic response solution, which implies that Hopf bifurcation arises in this case.3

In this section, we will rewrite the system (3) with the Gaussian white noise can be represented as follows

where x " p x 1 x 2 x 3 y 1 y 2 y 3 q T , ξptq " p 0 ξ 1 ptq 0 0 ξ 2 ptq 0 q T . The matrices in Eq.11 are given in Appendix B.

As piω´A 1´A2 e´i ωτ 0 qqp0q " 0, where qp0q is the eigenvector of the equation (A.1). Furthermore, qpθq " qp0qe iωθ . Defining that Φpθq "´φ 1 pθq φ 2 pθq¯T , where φ 1 pθq " Repqpθqq, φ 2 pθq " Impqpθqq.

With Euler's formula, we can obtain the expression of φ ij pi, j " 1, 2q, which is given in Appendix C.

According to the symmetry of the matrix A and B, we can calculate that φ 14 pθq " φ 11 pθq, φ 15 pθq " φ 12 pθq, φ 16 pθq " φ 13 pθq; φ 24 pθq " φ 21 pθq,

The relative of Φpθq and Ψpsq: ( Φpθq P Cpr´τ, 0s,

we can solve Ψpsq as following

Ψpsq "ˆψ 11 psq ψ 12 psq ψ 13 psq ψ 14 psq ψ 15 psq ψ 16 psq

The elements ψ ij psqpi " 1, 2; j " 1, 2, 3, 4, 5, 6q are shown in Appendix C. Next, we use the stochastic center manifold theorem to reduce the delay differential equation (3) to an ordinary equation [32, 34] . we should notice that φ k pθq P Cpr´τ, 0s, R 2 q and ψ j psq PĈpr0, τ s, R 2 q, j, k " 1, 2, the space C is spanned by the two-dimensional subspace P P C, which is driven by the eigenvalues λ 1,2 "˘iω of the equation F 1 pλq " 0 at Hopf bifurcation, plus the infinite-dimensional subspace Q P C associated with all the other eigenvalues of F 1 pλq " 0, namely C " P ' Q. The center manifold M β Ď Cpr´τ, 0s, R n q is constructed, which is tangent to P . Defining the bilinear pairing

The substitution of the inner product matrix p Ψpsq Φpθq q into the bilinear pairing, we can obtain the nonsingular matrix Ψ Φ˘n sg "ˆa 11 a 12 a 21 a 22Ṫ he elements of the matrix p Ψ Φ q nsg is given in Appendix C.

The normalization of Ψpsq toΨpsq P C is computed using the formulā Ψpsq " p Ψ Φ q´1 nsg ψpsq, following that Ψpsq "ˆψ 11 psqψ 12 psqψ 13 psqψ 14 psqψ 15 psqψ 16 psq ψ 21 psqψ 22 psqψ 23 psqψ 24 psqψ 25 psqψ 26 psq˙, 0 ď s ď τ.

Substituting of the new elements p ψ j psq φ k pθq q pj, k " 1, 2q of pΨpsq Φpθq q into the bilinear relation (12) as follow:Ψ

As for the definition of 9

Φpθq " ΦpθqB, where B " p 0 ώ ω 0 q as pΨpSq, Φpθqq id " I, then it can be express as Φpθq " Φp0qe Bθ ,´τ 0 ď θ ď 0 and Ψpsq " Ψp0qe Bs , 0 ď s ď τ . Thus, the solution of the linearized delay Eq.

x 2 pt´τ 0 q "u 1 ptq cos ωτ 0`u2 ptq sin ωτ 0 ; 9 x 2 pt´τ 0 q " 9 u 1 ptq cos ωτ 0`9 u 2 ptq sin ωτ 0 ;

The derivation process of the center manifold is shown in Appendix D.

Then, on the center manifold M f P C, the stochastic ODEs about uptq is obtained 

where Q " β N cos 2 θb 1 cos ωτ 0´s in 2 θb 2 sin ωτ 0`p b 2 cos ωτ 0´b1 sin ωτ 0 q cos θ sin θ pb`β N I˚q 2`ω2

.

It is obvious that R and ϕ are significant only when the order is O 1{ε , thus, the Itô stochastic differential equations including Wong-Zakai correction term for the averaged amplitude equation with the method of integral averaging [35] is shown below

where Bptq represents the independent unit wiener process, mpRq and σpRq represent drift coefficient and the diffusion coefficient, respectively, and mpRq and σpRq are defined below:

In the sequel of this section, we shall consider the stochastic stability and bifurcation of the SDDEs. With the Eq.(15), the local stability of the stochastic systems is analyzed by using the method of maximum Lyapunov exponent on the basis of Oseledec multiplicative ergodic theory [36] .Obviously, the system (15) is a geometric Brownian motion and the solution is Rptq " Rp0qexpppµ 1`µ 2 2 qt`?µ 2 Bptqq. Then the Lyapunov exponent of the Eq.(15) can be obtained

According to the Ref. [37] , we can infer that the stochastic D-bifurcation occurs when λ " 0.

In this section, numerical simulations were made to reveal the complex dynamics of the system (3). With the reveal of the analytical results in previous sections, numerical simulations are employed to analyze the effect of time delay on the nonlinear stochastic delay system. It can be found that time delay takes an active part in the system (3) as a time change function, in other words, the variation of the time delay value as a switch of the system (3) controls the stability of the system. Then, to better verify the validity and rationality of the model, two adjacent regions in China are taken as the research objects to carry out the simulation. One of the regions consists of Heilongjiang province, Jilin province and Liaoning province, and the three provinces are interrelated, known as Three northeastern provinces (TNP). Another is Beijing, Tianjin and Hebei province, marked as Beijing-Tianjin-Hebei Urban Agglomeration (BTHUA).

The parameters of the SIR model we used were within the range of their empirical estimates [4, [38] [39] [40] [41] . Under the conditions((H1), (H2) and (H3)) which are given before, we choose the parameters as β " 0.57d´1(0.09 to 1.12; refs.4, 31) , γ " 0.32d´1 (0.07 to 0. 29; refs.4,31,33) . The value of the fraction of infections without symptoms µ is 0.87, (ref. 34) . To achieve the best estimate of infection mortality rate (p1´µqα{pα`γq " 0.6%; ref.34) that α " 0.015d´1. And b " 0.01 and the rate of movement is κ " 0.006d´1. According to the values given above, we can calculate that the equilibrium point EpS˚, I˚, R˚q " p3053070.18, 87175.96, 1859753.86q and the stochastic Hopf bifurcation occurs at the critical value τ 0 " 4.53. When τ " 4.38 ă τ 0 , the equilibrium point EpS˚, I˚, R˚q is asymptotically stable.When τ ą τ 0 , the periodic solutions of the system (3) occur. These are illustrated in Fig.2 . 

Infectious disease can cause a large number of deaths and injuries in a short period of time in severe cases [9] . Therefore, it is necessary to take effective measures in time. For infectious diseases, the most important thing is to cut off the source of infection to prevent the spread of the disease.

It spreads to other cities of China as pneumonia cases from unknown origin. Now we apply the model to the two regions of China. Geographically, BTHUA and TNP are adjacent to each other and the transportation is convenient. Besides, the total population of these two areas is also similar to meet the model assumptions. Data Source: National Health Commission of the People's Republic of China with high precision.

In Fig.5 , a comparison between the Monte Carlo simulated curve and the real data for the TNP and BTHUA are displayed, almost all the points of the actual data stand on the Monte Carlo simulated curve in the first 180 days, the effectiveness of SIR model is revealed. Then the Monte Carlo simulated curve starts to break away from the real data points. The actual data is asymptotically stable and approach zero, while the direction of the estimated curve shows an obvious upward trend. That's because on January 23, 2020, the Chinese government took measures of "lockdown of cities" in Hubei and decreased population mobility in China.

"Lockdown of cities" means that residents can't travel between cities, counties, and communities in Hubei. Besides, on January 28, 2020, the traffic of Xiangyang (a city of Henan) has been cut off. At the same time, due to traffic control, media publicity, extended holidays and other factors, the scale of population flow between cities in other regions has been greatly reduced, to contain the rapid spread of the COVID-19 epidemic.

The probability density of the infected cases with the real data is shown in Fig.6 . The decline of the probability density curve with the two regions-BTHUA and TNP indicates that the effectiveness of decreasing population mobility. Comparing Fig.6(a) and (b) , obviously, different from the BTHUA has only one peak of probability density, TNP has two. The reason is that after the epidemic was basically under the control of China, and companies resume work and production across China, the TNP sent an alert for imported infectious cases broke out in the border city of China and Russia, followed by a local case where the source of the infection was unknown. On may 10, 2020, Shulan, a city in Jilin, was declared to be a high-risk area and entered "wartime state" . Moreover, Jilin cut off external traffic, limited the entry and exit of people. Which made Jilin, with a population of 4.4 million, has become the first major city in China to be shut down since the resumption of work and production in China.

Above all, as is known to all under the condition that emergency event occurs in one area the probability density peak of the two regions is different. To find without emergency event happened whether the epidemic peak in two adjacent regions is the same or not. We simulated 50 pairs of areas with adjacent geographical locations and convenient transportation. Assume there has 500 infected cases in every area at first, then simulating the curve of the infected cases with 300 days as shown in Fig.7 . Obviously, the epidemic peaks of the 100 curves in Fig.7 are different from each other, that's mainly because of the differences in movement rate κ among areas and the existence of the incubation period τ of the coronavirus. For explanation in detail, numerical simulation on two groups of these was carried out displayed in Fig.8 .

Note T i (T i ą 0 and T i P N , i " 1, 2) is time interval, which represents the time it takes for the epidemic peak from one area to another in group i. It should be noted that in Fig.8 , the delay in Group 1 denotes τ 1 " 6.8 and in group 2 is τ 2 " 10.6. Obviously, T 1 ă T 2 (T 1 « 9, T 2 « 15) and the long duration of a large outbreak in group 2 is harmful to healthy people. It can be explained by when the epidemic of the infected city spreads rapidly to another city, it can immediately attract the government's attention and quickly take measures to control the epidemic at the first time. From [30] we knew that the faster we took the measures, the more beneficial it is to reduce the number of infections and control the epidemic. Therefore, we can draw a conclusion that a large delay τ is not conducive to the control of the epidemic. Fixed the parameters α " 0.015, γ " 0.32, ε " 0.2, β " 0.57, µ " 0.87, b " 0.01, τ " 10.686, and change the movement rate κ P t0, 0.001, 0.0045, 0.0058, 0.01, 0.086, 0.092u, we investigate the influence of the epidemic from the infected area i to the uninfected area j in Fig.9 and Fig.10 and the areas of i and j are adjacent. In Fig.9 , the movement rate κ " 0 is an absolutely idealized situation as it requires strict enforcement by the government. Besides, in the infected area i which is locked down, the economy is bound to be traumatized by destruction and the basic civil rights of the healthy people will be deprived. Considering the situation that the population movement (κ ‰ 0) displayed in Fig.10 .

By comparing Fig.9 with Fig.10(a) , (b), (c) and (d), we find that proper population mobility can alleviate the number of infections in the infected area i. In terms of time, with appropriate movement rate uninfected people in the infected area i can be transferred to safe areas to separate the uninfected from the infected, bought time to control the epidemic. From the perspective of space, population density has been alleviated to avoid largescale outbreaks. As the frequency of population mobility in the two coupled regions in Fig.10(e) and (f), the trend of the epidemic in the two regions is synchronous. The almost simultaneous outbreak in two areas makes it difficult to control, and the high movement rate enhances the spread of the epidemic(see Fig.10(a), (b), (c) and (d) ). Region j has a negative effect on i, which exacerbates the epidemic in region i. In Fig.10(f) , the number of infections in both regions reaches a maximum, unfortunately.

The results in Fig.10 indicates that instead of decreasing monotonously as κ increases, T i first decreases then has a minimum in Fig.10(c) , after increases and finally decreases. Which is interesting and subverts our habitual thinking. We named the value of κ in Fig.10(c) is the "smartest" κ with τ " 10.686, due to the "smartest" movement rate is helpful to control the epidemic of area i. Besides, the total number of infections is lower than other cases and the trend of the epidemic in i and j is almost the same with T 3 « 12, which is conducive to epidemic control. As the prevention and control policy in area j can take i as a reference in large extent. This finding highlights the need to find the smartest κ under different values τ to control not only epidemic but also any pandemic.

Another interesting finding in Fig.10(a) and (b) is that even if the movement rate is low and the T is a long time, the epidemic peak of the uninfected area j before is seriously higher than the infected area i. Which indicates it is not absolutely safe as long as there are infected people, as the old saying goes-sparks can be bonfires and ruin everything. Such as on December, 2020, an infected person in Shenyang, China, who led to many people being infected and the entire city of Shenyang has been put on emergency quarantine. Moreover, due to only one infected person, Shijiazhuang, China took measures-"city lockdown " on January, 2021. Although, it has been more than one year since the outbreak in Wuhan, China. It warns us strengthen epidemic prevention and control in key regions. Precise and taking differentiated epidemic control strategies is crucial. Besides, epidemic monitoring should be strength, under normal conditions, once the outbreak occurs, the first-level public health emergency response can be activated immediately.

In this paper, we investigated the dynamic behavior of a delayed stochastic SIR epidemic model with Gaussian white noise in two coupled populations. Different from other models, our model takes time delay and environment noise into account. We find that the short incubation period is helpful to control the epidemic. As the small time delay with the short time interval, thus the epidemic can be detected at the first time, and prevention and control measures can be taken as soon as possible. Moreover, the population flow among municipalities with the convenient transportation network also considered. However, the spread of the epidemic is a complex progress, due to the country has high-speed and large-capacity transport facilities, so that the population can move quickly, flexibly and massively between cities. Complete physical isolation is impossible to achieve in real life, thus the movement rate should as low as possible? No. Interestingly, we find there exists an optimal population movement rate, which is so called "smart" κ. This finding could provide effective guidance for global epidemic prevention and control.

This study also demonstrated that we should prevent not only local infections but also external infections. In particular, preventing local to international flows and transmission is essential for the control of outbreaks and future pandemics. Curbing population flow with the "smartest" κ is the best way to contain the outbreak.

In this study, from the perspective of mathematical theory, assuming that the two objects in the research have the same population, it is beneficial to analysis the dynamic behaviors of the model and simplify the calculation. However, such cases are rare in real life. Nevertheless, this study provides a reasonable basis for future studies with different population sizes.

The results of the study only focus on the spread of the epidemic with adjacent cities of a country. The spread of the epidemic between multinational cities and countries will serve as the content of expanding the model into a chain system. Based on this model, we can add nodes as needed to carry out our research on a multi-area chain system, like the circle area surrounding the sea or a chained area that is geographically presented in a narrow form. It is possible to embed the system of time-delay stochastic differential equations with noise into the neural network, then we can establish mathematical models around the world to predict and control infectious diseases. Moreover, contributing to the health of people all over the world and researchers who are interested and have the same goal as us are also welcome to join us. At present, the international situation with COVID-19 remains serious. While the outbreak is almost under control in some countries, if the outbreak is not controlled globally, the epidemic still risks spreading across borders. And we must prevent secondary transmission from other regions. For a sudden case that may occur at any time, such environmental noise with jumping behavior can be described by Lévy noise. Many works have been done with Lévy noise to avoid and control the epidemics [42] [43] [44] .

In addition, the model of this article was built in the early stage of the outbreak. However, with the spread of the epidemic, this model is too simple to be applied to the study of COVID-19. This is because it does not consider some of the major processes of this epidemic. For instance, the effects of SARS-CoV-2 variants and vaccines on COVID-19 is not included. These factors complicate the study of epidemics, so one type of approach cannot solve these problems. Motivated by [45] [46] [47] [48] [49] the next related research should not use only one method of differential equation modeling, but a combination of various methods. Such as visualization and analysis of Spatio-temporal geographic data, cellular-automata method, gene sequencing techniques in biology, etc. Although most of people have been vaccinated with different types of COVID-19 vaccine, unfortunately, the virus is constantly mutating, which makes fighting the epidemic a long-term task. People need to communicate, and the world needs to recover normalcy, not a long-term lockdown. Therefore, the discovery of the "smartest" movement rate for the long-term battle is more meaningful. As is known to all, all models can only be infinitely close to the truth. But our work has taken us to a meaningful step toward the truth. Proof Assume iωpω ą 0q is a root of the equation F 1 pλq " 0. The form of F 1 piωq " 0 is as followś

With Eulers formula, we obtain the equations

then, we have

where h "´pb`γ`αµ`βq and 0 ă h ă γ Set z " ω 2 , then we obtain

Let z 3 and z 4 be the roots of (A.1), then by (H1) and (H2)

It demonstrates that the equation F 1 pλq " 0 has a pair of pure imaginary roots.2

The matrices in Eq.11 are shown as following:

Note, we record the elements in matrix A 2 as a ij . the matrix F can be shown as´ε β N x 1 x 2 pt´τ 0 q`ετ β N S˚x 1 2 pt´τ 0 q ε β N x 1 x 2 pt´τ 0 q`ετ p β N S˚´γqx 1 2 pt´τ 0 q 0 ε β N y 1 y 2 pt´τ 0 q`ετ β N S˚y 1 2 pt´τ 0 q ε β N y 1 y 2 pt´τ 0 q`ετ p β N S˚´γqy The elements φ ij pθqpi, j " 1, 2q in the matrix φpθq are shown below φ 11 pθq " 1 pb`β N qI˚`ω 2 rαµpb`β N I˚q cos ωθ`ωαµ sin ωθ´pb`β N I˚q β N S˚cos ωpτ 0´θ q`ω β N S˚sin ωpτ 0´θ qqs; φ 12 pθq " cos ωθ; φ 13 pθq " γ b 2`ω2 rb cos ωpτ 0´θ q`ω sin ωpτ 0´θ qs; φ 21 pθq " 1 pb`β N qI˚`ω 2 rαµpb`β N I˚q sin ωθ´ωαµ cos ωθ`pb`β N I˚q β N S˚sin ωpτ 0´θ q`ω β N S˚cos ωpτ 0´θ qqs; φ 22 pθq " sin ωθ; φ 23 pθq " γ b 2`ω2 rb sin ωpτ 0´θ q´ω cos ωpτ 0´θ qs;

The elements ψ ij psqpi " 1, 2; j " 1, 2, 3q of the matrix ψpsq are here ψ 11 psq " 1 pb`β N qI˚`ω 2 rαµpb`β N I˚q cos ωs`ωαµ sin ωs´pb`β N I˚q β N S˚cos ωpτ 0´s q`ω β N S˚sin ωpτ 0´s qqs; ψ 12 psq " cos ωs; ψ 13 psq " γ b 2`ω2 rb cos ωpτ 0´s q`ω sin ωpτ 0´s qs; ψ 21 psq " 1 pb`β N qI˚`ω 2 rαµpb`β N I˚q sin ωs´ωαµ cos ωs`pb`β N I˚q β N S˚sin ωpτ 0´s q`ω β N S˚cos ωpτ 0´s qqs; ψ 22 psq " sin ωs; ψ 23 psq " γ b 2`ω2 rb sin ωpτ 0´s q´ω cos ωpτ 0´s qs; ψ 14 psq " ψ 11 psq, ψ 15 psq " ψ 12 psq, ψ 16 psq " ψ 13 psq ψ 24 psq " ψ 21 psq, ψ 25 psq " ψ 22 psq, ψ 26 psq " ψ 23 psq

The elements of the matrix p Ψ Φ q nsg , following that a 11 "2r pb`β N I˚q 2 pαµ´β N S˚cos ωτ 0 q 2 ppb`β N I˚q 2`ω2 q 2`1`γ 2 pb 2`ω2 q 2 pb cos ωτ 0`ω sin ωτ 0 q 2 ś β N S˚1 pb`β N I˚q 2`ω2 rαµpb`β N I˚q pcos 2ωτ 0´1 q´ωpτ 0`s in 2ωτ 0 qq. a 22 "2r 1 ppb`β N I˚q 2`ω2 q 2 p´ωαµ`pb`β N I˚q β N S˚sin ωτ 0`ω β N S˚cos ωτ 0 q 2γ 2 pb 2`ω2 q 2 pb sin ωτ 0´ω cos ωτ 0 q 2 s´β N S˚1 pb`β N I˚q 2`ω2 ppb`β N I˚qαµ pτ 0 cos ωτ 0´1 ω sin ωτ 0 q`ωαµτ 0 sin ωτ 0´2 pb`β N I˚q β N S˚pτ 0`1 2ω sin 2ωτ 0 qq´1 2 β N S˚p1´cos 2ωτ 0 q´p β N S˚´γqp 1 ω sin ωτ 0´τ0 cos ωτ 0 q γ 2 b 2`ω2 pbpτ 0´1 2ω sin 2ωτ 0 q´1 2 p1´cos 2ωτ 0 qq.

where G is G " pb`β N I˚q 2 1 ppb`β N I˚qppαµ´β SN cos ωτ 0 q`ω βSN sin ωτ 0 qp´ωαµp b`β N I˚q βSN sin ωτ 0`ω βSN cos ωτ 0 q`p b 2`ω2 q 2 γ 2 pb cos ωτ 0`ω sin ωτ 0 q pb sin ωτ 0´ω cos ωτ 0 q Appendix E. Description of the parameters

The parameters in Eq.(13) are shown below b 1 "pb`β N I˚qpαµ´β N S˚cos ωτ 0 q`ω β N S˚sin ωτ 0 . b 2 "´ωαµ`pb`β N I˚q β N S˚sin ωτ 0`ω β N S˚cos ωτ 0 q. J 11 "ψ 11 p0q " a 22 b 1´a12 b 2 pb`β N I˚q 2`ω2 . J 12 "ψ 12 p0q " a 22 .

J 21 "ψ 21 p0q "´a 21 b 1´a11 b 2 pb`β N I˚q 2`ω2 . J 22 "ψ 22 p0q "´a 21 .

α 1 " pu 1 b 1`u2 b 2 qppu 1 cos ωτ 0´u2 sin ωτ 0 q pb`β N I˚q 2`ω2 .

A "p1´α 2 cos ωτ 0 qα 3´p 1´α 3 cos ωτ 0 qα 2 . B "p1´α 2 sin ωτ 0 qα 3´p 1´α 3 cos ωτ 0 qα 2 sin ωτ 0 .

Whereτ " τ´τ 0 ε

Aerosol transmission of infectious disease

Coronavirus disease 2019 (covid-19): situation report

Impact of international travel and border control measures on the global spread of the novel 2019 coronavirus outbreak

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (sars-cov-2)

Explicit formulae for the peak time of an epidemic from the sir model. which approximant to use?

Explicit formulae for the peak time of an epidemic from the sir model

Infectious diseases of humans: dynamics and control

A contribution to the mathematical theory of epidemics

Mathematical biology II: spatial models and biomedical applications

Small world and scale free model of transmission of sars

The inflection point about covid-19 may have passed

Effect of non-pharmaceutical interventions to contain covid-19 in china

Transmissibility of 2009 pandemic influenza a (h1n1) in new zealand: effective reproduction number and influence of age, ethnicity and importations

Novel coronavirus 2019-ncov: early estimation of epidemiological parameters and epidemic predictions

Nowcasting and forecasting the potential domestic and international spread of the 2019-ncov outbreak originating in wuhan, china: a modelling study

Quarantine in a multi-species epidemic model with spatial dynamics

Exploring the roles of high-speed train, air and coach services in the spread of covid-19 in china

Dynamical behaviour of an epidemic on complex networks with population mobility

Dynamics of the impact of twitter with time delay on the spread of infectious diseases

A periodic two-patch sis model with time delay and transport-related infection

The analysis of an epidemic model with time delay on scale-free networks

Global stability analysis of an sir epidemic model with demographics and time delay on networks

A threshold of a delayed stochastic epidemic model with crowly-martin functional response and vaccination

Dynamics of a stochastic viral infection model with immune response

Threshold dynamics of a time-delayed seirs model with pulse vaccination

An seir epidemic model with constant latency time and infectious period

The interplay of movement and spatiotemporal variation in transmission degrades pandemic control

Spread of covid-19 in china: analysis from a city-based epidemic and mobility model

Dynamics for an seirs epidemic model with time delay on a scale-free network

Weak signal detection based on mathieuduffing oscillator with time-delay feedback and multiplicative noise

Elements of applied bifurcation theory

Delay dynamical systems and applications to nonlinear machine-tool chatter

On stochastic processes defined by differential equations with a small parameter

Lyapunov exponents of linear stochastic functional differential equations driven by semimartingales. part i: the multiplicative ergodic theory

The Hopf bifurcation and its applications

Inferring change points in the spread of covid-19 reveals the effectiveness of interventions

The challenges of modeling and forecasting the spread of covid-19

Spread and dynamics of the covid-19 epidemic in italy: Effects of emergency containment measures

The implications of silent transmission for the control of covid-19 outbreaks

A stochastic sirs epidemic model incorporating media coverage and driven by lévy noise

Stochastic sirs model driven by lévy noise

Dynamics of a stochastic sis model with double epidemic diseases driven by lévy jumps

A time-delayed sveir model for imperfect vaccine with a generalized nonmonotone incidence and application to measles

Sars-cov-2 susceptibility and covid-19 disease severity are associated with genetic variants affecting gene expression in a variety of tissues

Spatio-temporal modelling of covid-19 incident cases using richards curve: an application to the italian regions

A data-driven understanding of covid-19 dynamics using sequential genetic algorithm based probabilistic cellular automata

A model based on cellular automata for investigating the impact of lockdown, migration and vaccination on covid-19 dynamics

The center manifold M f " tφ P C|φ " Φu`hpuq, h P Su P C, and the solution of Eq.(11) on the center manifold as followwhere Xptq " px 1 , x 2 , x 3 , y 1 , y 2 , y 3 q T and θ P r´τ 0 , 0s. The equation (11) can be expressed as an extended equation of X t pθq to compute the central manifoldwhere Lrx t pθqs is a linear operator of the equation (11) and F rx t pθqs`ε 1{2 Grx t pθqs is a part of nonlinear of Eq.(11). By solving the Eq.(D.1) and Eq.(D.2), we obtain rΦpθq`D u hpθ, uptqqs 9 uptq "Due to xΨpsq, hpθ, uptqqy " 0, we can obtain the ODEs about uptq, namely 9 uptq " Buptq`Ψp0qF rΨpθquptq`hpθ, uptqqs`ε 1{2 GrΦpθquptq`hpθ, uptqqswhere uptq " pu 1 ptq, u 2 ptqq. Furthermore, we can expand it on the central manifold.