There exists the "smartest" movement rate to control the epidemic rather than "city lockdown".

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" κ , which is helpful to control the epidemic. This model is also useful to prevent other infectious diseases.

Introduction

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 was first discovered in late 2019. In common with other infectious diseases (such as influenza, cholera and meningitis), it can infect a large number of people in a short time. The disease symptoms like a cold, a flu or pneumonia, vary from that of uncomplicated upper respiratory tract infection to severe viral pneumonia with multiple organ failure and death. It is mainly transmitted by the air (airborne transmission), when the infected people during activities such as exhaling, coughing, sneezing and talking, aerosols and droplets contaminated with SARS-CoV-2 are formed, leading to infection among susceptible populations [1].

On January 30, 2020, according to the International Health Regulations, the World Health Organization (WHO) listed the SARS-CoV-2 epidemic as a Public Health Emergency of International Importance, which is the highest-level alert in WHO. On March 11, 2020, novel coronavirus’s blood disease was defined as a pandemic by WHO [2], [3], [4]. Creating innovative solutions to fight the epidemic, both inside and outside health institutions around the world, can minimize social and economic impacts [5], [6]. For example, predictive models can use data to predict the future growth and decline of the epidemic, and generate scenarios to guide territorial management agencies in decision-making. Recently, through Refs. [7], [8], we found a method that can directly estimate the epidemic peak, which is beneficial for the prevention and control of the epidemic. A crucial point of studying epidemic phenomena is that we can build models from the data we have to make a reasonable prediction about the future.

Since Hamer and Ross first conducted a quantitative survey of the spread of measles, scientists have been working on developing mathematical models better to investigate the spread of the epidemic [9]. The SIR model is the most classical model in epidemiology and disease control, which was first proposed and investigated by Kermack and McKendrick [10]. Since then, variants of the SIR model (such as SIRS and SEIR) have emerged for theoretical and practical purposes [11], [12]. These models [13], [14], [15], [16], [17] are the mainstream of epidemiological dynamics studies, as they are usually used to predict the peak duration or inflection point of epidemics effectively. However, they have a common weakness: the spread of epidemics in space has been largely ignored [18]. Therefore, the SIR model cannot be applied to every city when simulating the epidemics for the whole country. In addition, it is unreasonable that set models for each city as every city is an open and connected system. Especially with the improvement of transportation technology and the coverage of transport network, the population movement rate 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 steady-state 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 Liu et al. [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 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.

The analysis of the model

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 , and (=1,2) represent the local susceptible, infectious and recovered individuals in area , respectively( is local population size). To make sure the local (i and j) population remains unchanged we assume that the birth rate is equal to the natural mortality rate, the disease mortality rate is equal to the supplemental population rate. , and (=1,2) are recovery, movement, and transmission rates, respectively. () is incubation period, which is a time period. is the infections without symptoms and is infected people with symptoms who are immobile (, represents the infected people in area ).

Let , , , we have , 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 () or die out (). According to Ref. [31], the basic reproduction number is calculated to be . When , the system (2) exists positive equilibrium, we assume that and the proportion of the local susceptible, infectious and recovered individuals are same in the two places ( and ). Therefore, we obtain the equilibrium point , 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 .

Fig. 1

A schematic diagram of the model.

Next, we shall consider the disturbance of environmental noise near the equilibrium pointwhere and represent the Gaussian white noise with power spectral density and , is independent of , is the scale coefficient which is small to ensure the micro disturbance to the system. Note that is small, and it scales to the same order as .

Let The linearized equation of system (3) without Gaussian white noise is

This leads to the characteristic equation are as following:

Due to , , the roots of the characteristic Eq. (5) are the solutions of equations  and .

Assume that (H1) and . whereUnder the condition (H1), it can be obtained by calculation and , thus 0 is not the solution of and

Assume that (H1) holds. Then the equilibrium ofof the system(3)withis asymptotically stable.

When , the Eq. (6) is equivalent to the following equationLet and be two roots of the Eq. (8). Under the condition (H1), thenAs for , we haveLet and be two roots of the . Under the condition (H1), thenThese imply that Eqs. (8) and (9) have no positive roots, hence all the roots of the characteristic Eq. (5) have negative real parts. □

Assume further

(H2) and When , we have the following conclusion.

If (H1) and (H2) hold, the equation  has no purely imaginary roots for all .

When , let () be the root of the equation . Through the work [32], we can obtain that satisfies the following equation:where . Set , then

Let and be the roots of (10), then by (H1) and (H2)This implies when , equation  has no purely imaginary roots. □

When, under the condition (H1) and (H2), the equation  has a pair of purely imaginary roots ( ) for .

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  has a pair of pure imaginary roots. Furthermore, we can solve and , as followwhere

Define , thus is a pair of simple imaginary roots of when .

Assume (H1) and (H2) hold, the equilibrium pointof the system(3)occurs stochastic Hopf bifurcation for.

Due to the sign of and is the same. Differentiating the two sides of the equation  with respect to , assume that (H1) and (H2) hold, we have the Hopf’s transversality conditionWhereAccording to the Hopf bifurcation theory [33], combining the above results and Lemma 2.2, the Eq. (3) has periodic response solution, which implies that Hopf bifurcation arises in this case. □

The reduction and stochastic stability and bifurcation of the SIR model with time delay in the point of Hopf bifurcation

In this section, we will rewrite the system (3) with the Gaussian white noise can be represented as followswhere , . The matrices in Eq. (11) are given in Appendix B.

As , where is the eigenvector of the Eq. (A.1). Furthermore, . Defining that , where , .

With Euler’s formula, we can obtain the expression of , which is given in Appendix C.

According to the symmetry of the matrix A and B, we can calculate that , , ; , , , thus

The relative of and : ()we can solve as following

The elements are shown in Appendix C.

Next, we use the stochastic center manifold theorem to reduce the delay differential Eq. (3) to an ordinary equation [32], [34]. we should notice that and , the space is spanned by the two-dimensional subspace , which is driven by the eigenvalues of the equation  at Hopf bifurcation, plus the infinite-dimensional subspace associated with all the other eigenvalues of , namely . The center manifold is constructed, which is tangent to . Defining the bilinear pairing

The substitution of the inner product matrix into the bilinear pairing, we can obtain the nonsingular matrix

The elements of the matrix is given in Appendix C.

The normalization of to is computed using the formula , following that

Substituting of the new elements of into the bilinear relation (12) as follow:

As for the definition of , where as , then it can be express as and . Thus, the solution of the linearized delay Eq. (3) is the

Considering that the unique solution (where ) of the original nonlinear delay Eq. (2), there must have and . It should be noted that and belong to the space ; and belong to the space . The solution of Eq. (2) through on the center manifold is given by the equation . Under the change of variable , then the first order approximation in , for .

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

Then, on the center manifold , the stochastic ODEs about is obtainedWhere and the parameters of the equation are given in Appendix E.

By the polar coordinates transformationwhere and are the amplitude and phase of the solution, we can obtain the stochastic ODEs with and where .

It is obvious that and are significant only when the order is , thus, the stochastic differential equations including Wong-Zakai correction term for the averaged amplitude equation with the method of integral averaging [35] is shown belowwhere represents the independent unit wiener process, and represent drift coefficient and the diffusion coefficient, respectively, and and are defined below: where

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 . 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 .

Numerical simulation

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).

Simulating dynamics of the SIR model

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.09 to 1.12; refs.4, 31), (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 (; ref.34) that . And and the rate of movement is . According to the values given above, we can calculate that the equilibrium point and the stochastic Hopf bifurcation occurs at the critical value . When , the equilibrium point is asymptotically stable. When , the periodic solutions of the system (3) occur. These are illustrated in Fig. 2 .

Fig. 2

Effect of time delay on the infected cases (for a more intuitive display, the ratio of the number of infected people to the total population is taken as the vertical ordinate): (a); (b).

In Figs. 3 and 4 , the top right plane of two figures show the marginal probability density of and the marginal probability density of in the top left. And the joint probability density of with , see bottom left plane of Figs. 3 and 4, the projection of the joint probability density are shown in the bottom right plane. From the perspective of probability density, we can see that from Figs. 3 to 4, the properties of the probability density have changed. It demonstrated that the stochastic D bifurcation happened.

Fig. 3

Effect of time delay on susceptible(S) and infected(I), respectively and joint with .

Fig. 4

Effect of time delay on susceptible(S) and infected(I), respectively and joint with .

The transmission of COVID-19 between two regions

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.

Fig. 5

Comparison of time series of infected cases for actual data and the Monte carlo simulation with SIR model.

“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.

Fig. 6

The actual probability density curves of BTHUA and TNP.

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 .

Fig. 7

Time evolution of the infectious cases for 50 pairs of adjacent areas.

Fig. 8

Time evolution of the infectious cases for two pairs of adjacent areas.

Note ( and , ) is time interval, which represents the time it takes for the epidemic peak from one area to another in group . It should be noted that in Fig. 8, the delay in Group 1 denotes and in group 2 is . Obviously, (, ) 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 Wei et al. [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 , and change the movement rate , we investigate the influence of the epidemic from the infected area to the uninfected area in Figs. 9 and 10 and the areas of and are adjacent. In Fig. 9, the movement rate is an absolutely idealized situation as it requires strict enforcement by the government. Besides, in the infected area 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 () displayed in Fig. 10.

Fig. 9

Ideal state: The influence of infected area on uninfected area with movement rate .

Fig. 10

The influence of infected area on uninfected area with different movement rate .

By comparing Figs. 9 with 10(a), (b), (c) and (d), we find that proper population mobility can alleviate the number of infections in the infected area . In terms of time, with appropriate movement rate uninfected people in the infected area 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 large-scale 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 has a negative effect on , which exacerbates the epidemic in region . 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, 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 , due to the “smartest” movement rate is helpful to control the epidemic of area . Besides, the total number of infections is lower than other cases and the trend of the epidemic in and is almost the same with , which is conducive to epidemic control. As the prevention and control policy in area can take 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 is a long time, the epidemic peak of the uninfected area before is seriously higher than the infected area . 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.

Conclusion and discussions

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 Lvy noise. Many works have been done with Lvy 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 Al-Darabsah [45], D’Antonio et al. [46], Mingione et al. [47], Ghosh and Bhattacharya [48], Jithesh [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.