Urban spatial epidemic simulation model: A case study of the second COVID-19 outbreak in Beijing, China.

The second COVID-19 outbreak in Beijing was controlled by non-pharmaceutical interventions, which avoided a second pandemic. Until mass vaccination achieves herd immunity, cities are at risk of similar outbreaks. It is vital to quantify and simulate Beijing's non-pharmaceutical interventions to find effective intervention policies for the second outbreak. Few models have achieved accurate intra-city spatio-temporal epidemic spread simulation, and most modeling studies focused on the initial pandemic. We built a dynamic module of infected case movement within the city, and established an urban spatially epidemic simulation model (USESM), using mobile phone signaling data to create scenarios to assess the impact of interventions. We found that: (1) USESM simulated the transmission process of the epidemic within Beijing; (2) USESM showed the epidemic curve and presented the spatial distribution of epidemic spread on a map; and (3) to balance resources, interventions, and economic development, nucleic acid testing intensity could be increased and restrictions on human mobility in non-epidemic areas eased.

INTRODUCTION

A series of atypical pneumonia cases first detected in Wuhan, China, in December 2019, has developed into a global pandemic (Huang et al., 2020; Tang et al., 2020; Wang et al., 2020; World Health Organization, 2020). The pathogen responsible is the severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2), and the disease was named coronavirus disease 2019 (COVID‐19) by the World Health Organization (Bai et al., 2020; Wu, Leung, & Leung, 2020). China's government implemented a series of non‐pharmaceutical interventions that effectively controlled the first wave (Ferretti et al., 2020; Zhou et al., 2020). After 56 days without a new case, new cases were reported in Beijing on June 11. From June 11 to July 5, 335 cases were reported (Wu, Wang, et al., 2020). The source of the outbreak was the Xinfadi Wholesale Market (“Xinfadi market”). Through targeted large‐scale testing and precise human mobility restrictions, the second outbreak in Beijing was stopped, and a second pandemic was avoided. Until mass vaccination achieves herd immunity, many cities are at risk of similar outbreaks. Therefore, we used the second COVID‐19 outbreak in Beijing as an example to develop an accurate spatial epidemic simulation model to quantify and simulate Beijing’s non‐pharmaceutical interventions.

Previous modeling studies have used publicly available population migration data (e.g., Tencent location request data and Baidu migration data) to develop a series of SEIR (susceptible, exposed, infected, and removed)‐based models for predicting the size of the pandemic (Chinazzi et al., 2020; Gilbert et al., 2020; Wu, Leung, et al., 2020), evaluating the effectiveness of non‐pharmaceutical interventions (Ferretti et al., 2020; Jia et al., 2020; Kraemer et al., 2020; Tian et al., 2020), and inferring transmission characteristics (Kissler, Tedijanto, Goldstein, Grad, & Lipsitch, 2020; Li et al., 2020). However, these models tend to minimize spatial heterogeneity and ignore that epidemics have different infectious effects due to different locations (Arino, Brauer, van den Driessche, Watmough, & Wu, 2007). Therefore, they produce statistical values or large‐scale infection spatial distribution maps, and it is difficult to change the research scale to the inner city.

To simulate the epidemic in the city, researchers used some fine‐grained data to highlight the spatial heterogeneity of the intra‐city epidemic spread in the model. For example, Zhou et al. (2020) used mobile phone signal data to estimate human mobility and the force of infection within the city. Based on this, they established an intra‐city epidemic model to predict the epidemic curve of COVID‐19 in Shenzhen, China. This model quantified the difference in the epidemic transmission caused by intra‐city human mobility. However, it was still an extension of the traditional SEIR model, which lacks the movement and interaction of individuals in the city. Zhou, Zhou, Zheng, and Lu (2021) used an agent‐based model to establish an interaction system for intra‐city individuals. They combined it with the SEIR model to simulate the spatio‐temporal process of COVID‐19 transmission in Guangzhou. This model described the interaction behaviors (e.g., exposure, movement, immunity) of multi‐type individuals. Still, it used small‐group agents to represent a large population and did not consider the transmission differences caused by different urban facilities. Koo et al. (2020) constructed multi‐type human agents based on Singapore demographics (e.g., age, marital status, religion, ethnicity) and then assigned agents to different urban facilities, simulating the spread of the epidemic in Singapore under various interventions. This model proposed a relatively complete framework of the intra‐city spatial epidemic model, but it used census data rather than precise mobility data, and the model is not validated by the actual data. The existing intra‐city epidemic models have intrinsic limitations that restrict their application to intra‐city spatial epidemic spread simulation. Further, their focus was the first pandemic wave. The sources and cause of the intra‐urban outbreaks are unknown, and detailed epidemic data are lacking, making the models inapplicable to the actual case.

In this study, with the help of mobile phone signaling data and point‐of‐interest (POI) data, we used a multi‐agent system technique to build a dynamic module of infected case movement within the city. We also incorporated the epidemic transmission module to establish an urban spatially epidemic simulation model (USESM). The aims of this study were to: (1) develop a USESM; (2) analyze and quantify Beijing’s intervention policies; (3) simulate the outbreak and verify the results; and (4) assess the impact of intervention policies in different scenarios on the epidemic.

METHODS

Model design

Most epidemic transmission models are deterministic (Jia et al., 2020; Kraemer et al., 2020; Tian et al., 2020), assuming that if all the parameters and initial values in the model are known precisely, the spread of epidemics can be accurately predicted at any time. This assumption is inconsistent with reality. The spread of epidemics is affected by random events, and small changes in the transmission process lead to entirely different results. We used a multi‐agent system technique to replace human social interaction, decision probabilities to replace real‐world coincidences, and roulette algorithms to translate decision probabilities into human actions to reflect this impact.

We adopted the definition of the compartment in the SEIR model (Kermack & McKendrick, 1991; Li, Graef, Wang, & Karsai, 1999), classifying the population into four agents—susceptible, exposed, infected, and removed in the multi‐agent system. The study area was divided into several sub‐regions, and mobile phone signaling data and POI data were used to calculate the decision probability of every agent moving, staying, and being infected in the sub‐regions. The roulette algorithm translated these decision probabilities into specific agent actions. Decisions with a higher probability have a greater chance of being executed, while those with a smaller probability still have a slight chance. We set the sub‐regions of the model as Beijing sub‐districts (a sub‐district is one of the smaller administrative divisions of China, it is a form of township‐level division which is typically part of a larger urban area), comprising 128 units. Data from the outbreak were used to validate the model’s accuracy and usefulness.

Model framework

The USESM framework is summarized in Figure 1. This model integrates the initialization module as well as the infection movement module and the epidemic transmission module. In the initialization module, the study area is divided into several sub‐regions and generates four types of agents (S: susceptible, E: exposed, I: infected, R: removed) in the sub‐regions (Step 1). The initial module generates a SIER dataset to store all model data, including current time, infection period, incubation period, agent status, and quantity. In the infection movement module, the movement possibility of each infection is calculated, and a binary roulette is built based on this probability to determine whether the infection can move (Step 2). The model calculates the probability of moving to each sub‐region for infections requiring movement and uses a multivariate roulette to determine the destination (Step 3). Moreover, the model calculates the probability that the infection will stay at the destination. If this decision probability is successfully converted into agent action, the infection will stay at the destination until removed. Otherwise, the infection will return to the origin sub‐region on the same day (Step 4). In the epidemic transmission module, POI and mobile phone signaling data are used to correct the basic reproduction number (R0) in different sub‐regions (Step 5) and then calculate the number of the next generation exposed through a random draw from a Poisson distribution for each infection (Step 6) (Lai et al., 2020). After this step, the temporarily moved infections will return to the origin sub‐region to avoid interference with subsequent epidemic transmission simulation. At the end of the epidemic transmission module, Bernoulli trials are used for each exposed person and infections to acquire the next generation infections and removed persons (Steps 7 and 8) (Lai et al., 2020). Finally, the model checks whether the SEIR data reached the termination condition: if not, the SEIR dataset will proceed to the next iteration (Time = time + 1); if so, the model will output the SEIR dataset for the entire simulation period.

FIGURE 1

Overall framework of urban spatial epidemic simulation model

Infected agent movement within the city

The movement of infections within the city differs from the movement across cities and countries, so we could not use the approach that allocated infections based on the ratio of the population going to the destination to the local population. Even if the infected persons have mobile behavior in the city, their destinations and stay times differ. We defined infected agent movement in terms of three perspectives: movement willingness, movement direction, and stay time.

Movement willingness

Before quarantine, the mobile behavior of the infected agent is no different from that of healthy people. We used movement probabilities of urban dwellers to quantify the movement willingness of the infected agents. The movement probability is related to the region population and number of people traveling, then the movement probability of each infection, as follows: where is the movement probability of each infected agent in sub‐region i at time t; is the number of the population moving from sub‐region i at time t; and is the number of the local population in sub‐region i at time t.

Movement direction

For those infections requiring movement, the direction is related to the number of the population in their sub‐region going to the destination sub‐region. To explore the movement direction of each infected agent, we calculated the probability of all candidate sub‐regions before the multivariate roulette algorithm determines the destination. The probability of each infected agent moving to the candidate sub‐region is as follows: where is the probability of each infected agent in sub‐region i moving to sub‐region n at time t; is the number of the population moving from sub‐region i to sub‐region n at time t; and is the number of the local population in sub‐region i at time t. At the current time step, if there is no infection in sub‐region i or no one has left, then all movement direction probabilities of sub‐region i are 0.

Stay time

The difference between intra‐city and inter‐city movements is that most intra‐city travelers will return to their starting point on the same day (i.e., daily commute). Therefore, most of the moving infected agents in the model will return to the starting sub‐region on the same day. But some exceptions stay at the destination for a long time (e.g., relocating to a different house, changing jobs). To calculate the probability of special cases, we propose the probability of the infection staying at the destination, as follows: where is the probability of each infected agent moving from sub‐region i to sub‐region n at time t and staying at sub‐region n until isolated; is the number of the population who moved from sub‐region i to sub‐region n at time t and stayed in sub‐region n longer than the infection period; and is the number of the population moving from sub‐region i to sub‐region n at time t.

All parameters in Equations (1)–(3) include local population [], moving population [,], and staying population [], calculated from mobile phone signaling data. The probabilities calculated from these parameters are conditional probabilities for each infected agent. The roulette algorithm determines whether the model can execute the conditional probability. Additionally, according to the probability order of movement willingness, movement direction, and stay, only when the model successfully executes the previous conditional probability can the following one enter the roulette algorithm.

R0 correction based on human mobility and POIs

COVID‐19 transmission efficiency is affected by human mobility and transmission places (e.g., transportation, supermarkets, farm product market, community, and workplaces) (Hu, Nigmatulina, & Eckhoff, 2013). We proposed an R0 correction method based on mobile phone signaling data and POI data to quantify these impacts in a city simulation.

Mobile phone signaling data

We extracted each trajectory's starting point, ending point, and length from the mobile phone signaling data. We determined the home location of the mobile phone users based on their ending point at night. The sum of home locations is the daily local population in each sub‐region, divided by the area of the sub‐region to obtain the population density (Equation 4). According to the starting and ending points of the mobile phone users' trajectories, we can determine the sub‐regions of their activities, then sum the number of active mobile phone users in each sub‐region as the daily active population. The ratio of daily active population to the local population can reflect the pressure of the population flow on the local population (Equation 5). After determining the sub‐region where the trajectories' starting points are located, the sum of all trajectories' lengths starting from each sub‐region can be calculated (Equation 6). Finally, we used population density (Equation 4), human mobility pressure (Equation 5), and travel lengths (Equation 6) to construct a dimensionless index (Kraemer et al., 2015). LogNormalize (Equation 7) was used to scale the index to the interval [0,1] to reflect the impact of human mobility on COVID‐19 transmission efficiency.

POI data

We summed the number of some POI types in each sub‐region, then divided by the area of the corresponding sub‐region to obtain the density of different POI types. Then a series of weight coefficients were constructed to combine different POI densities to reflect the impact of POI on COVID‐19 transmission efficiency (Equation 8).

Finally, we grouped the impact factors of human mobility and POIs to correct R0 in different sub‐regions (Equation 9). The maximum and minimum values of R0 in Equation (9) are theoretical values obtained from similar studies (Koo et al., 2020; Zhou et al., 2020). where is the number of the local population in sub‐region i at time t; is the population density of sub‐region i at time t; is the area of sub‐region i; is the population flow pressure of sub‐region i at time t; is the number of the active population in sub‐region i at time t; is the sum of all travel lengths starting from sub‐region i; is the sum of all travel lengths from sub‐region i to sub‐region n; is the impact index of human mobility in sub‐region i at time t; is the impact index of POI in sub‐region i; is the POI density of type j in sub‐region i; is the weight coefficient of POI type j, ; is the corrected R0 of sub‐region i at time t; and are the theoretical minimum and maximum values of R0; and and are the impact factors of POI and population movement, .

Epidemic transmission module

After infection movement and R0 correction, we can obtain the number of infections and R0 for each sub‐region on a given day. The model calculated the newly exposed people due to infectious people through a random draw from a Poisson distribution (mathematical expectation is R0) and multiplied by the proportion of susceptible people in the city (Equation 10) (Lai et al., 2020). Then, the next generation exposed for each sub‐region is the sum of new and existing exposed people (Equation 11). The model converted exposed people to infections by Bernoulli trials (the success probability is 1/incubation period) for exposed people (Lai et al., 2020), and the sum of new and existing infectious people is the next generation of infectious people (Equation 12). Similarly, the success probability of Bernoulli trials for quarantining infectious people is 1/infection period (Equation 13). where is the newly exposed people of sub‐region i at time t + 1; Random.poisson is a function of random sampling from the Poisson distribution, where λ is the mathematical expectation of the Poisson distribution; n is the number of random samples, and the result of this function is the sum of multiple samples; and are the exposed people of sub‐region i at time t and t + 1, respectively; and are the susceptible people and total population of sub‐region i at time t; and are the infectious people of sub‐region i at time t and t + 1, respectively; and are the removed people of sub‐region i at time t and t + 1, respectively; Random.binomial is a function of Bernoulli trials, where p is the probability of successful experiment (onset; quarantined); n is the number of experiments; and the result of this function is the sum of multiple experiments.

Data profile and processing

Mobile phone signaling data were acquired from China Mobile Ltd (Beijing, China). The data include phone number users (locals and tourists) active in Beijing from May 5 to June 30, 2020. Whenever a mobile phone user makes a mobile phone call, accesses the internet, sends or receives a text message, or moves across base stations, it is recorded by the nearest base station. This record includes user identification, time, and location of the base station. There are over 110,000 base stations in Beijing (Figure 2); we collected more than 120 million records to approximate the population and human mobility.

FIGURE 2

Map showing the geographic location, number of base stations in different sub‐districts of Beijing, and distribution of cumulative infection cases in Beijing (July 5, 2020)

Daily population

The local population size is the crucial parameter of the epidemic simulation model. Previous models were limited by data and only treated it as a fixed value, negatively impacting results (Kraemer et al., 2020). The local population of a city changes with mobility. To count the daily local population of every sub‐region, we defined the base station where a mobile phone user appeared most frequently between 10 p.m. and 6 a.m. as the user's home location, then counted the number of home locations in each sub‐region as the daily local population.

Human mobility

Human mobility is the number of people moving within or between different regions (e.g., cities, districts, sub‐districts) each day. To quantify human mobility, we performed trajectory segmentation (time threshold over 30 min, distance threshold over 200 m) on each trajectory in the mobile phone signaling data. We then extracted the midpoint of each segmentation as the stay location (Xu, Belyi, Bojic, & Ratti, 2018). We used these stay locations as the starting and ending point of each movement of the mobile phone user. Then, we counted the number of trajectories with the starting and ending point in the same region to measure the human mobility within the region. Similarly, we calculated the number of trajectories with starting and ending points in different regions to measure human mobility between different regions.

Point‐of‐interest

The POI data were crawled from the Gaode Map Application Programming Interface (https://lbs.amap.com/), including markets (agricultural, seafood, vegetable, and supermarkets), subway and bus stations, workplaces, and health facilities (hospitals, health centers, clinics, and pharmacies). Gaode Map does not explicitly state the time information of the POI data. However, as the Gaode Map is a navigation software, its POI data are updated in real time. Therefore, the time information of the POI data crawled from the Gaode Map should be consistent with the collection time (i.e., August 2020).

Basic parameters and scenario establishment

Basic parameters

Distribution of the infected person onset

The number of daily new confirmed infection cases and the infection period of 201 cases were obtained from the government website of Beijing (Liu, Gayle, Wilder‐Smith, & Rocklöv, 2020). We excluded some outliers (infection period data after June 14) and calculated the average infection period (time from onset to confirmation) for all cases to be 3 days. Then, we estimated the complete distribution of infected person onset by fitting a log‐normal distribution (the mean value was the average infection period, i.e., 3 days) to data on confirmed histories and onset dates in a subset of cases.

Simulation period

Our model aims to simulate the number of infections (i.e., the number of individuals who can spread infectious diseases). In this study, we hypothesized that only post‐onset individuals could spread infectious diseases. Therefore, we need to infer the earliest onset date and the latest in this outbreak. We estimated 1,000 groups for the distribution of daily onset cases based on publicly available epidemiological data (information about the onset date in publicly available epidemiological data is incomplete). In these distributions, the earliest onset date was May 28, and the latest was June 28. Therefore, our simulation period was May 28 to June 28, 2020.

Incubation time

The incubation time was 5.2 days (95% CI: 4.1–7.0 days) (Lai et al., 2020).

R0

The outbreak in Beijing was not severe, and there were significant differences in the transmission intensity of different sub‐districts. Therefore, we set to be 1.5 and to be 2.5 (Zhou et al., 2020). Using the R0 correction method, R0 was close to 2.5 in areas of high transmission risk and close to 1.5 in areas of low transmission risk. For most areas, R0 was close to 2, representing a state of moderate and likely outbreak (Koo et al., 2020).

Infection period

In the model, we define the infection period as the time during which an infected person participates in the transmission process. We then break this time into three parts (Equation 14) for quantification. The reciprocal of the infection period is the removal rate (Equation 15), which represents the probability of quarantining infected people from the epidemic transmission process (Liu et al., 2020). where is the time from onset to participation in nucleic acid testing; is the time from participation in nucleic acid testing to the publication of testing results; and is the time from publication of testing results to isolation.

We can acquire the infection period from the epidemiological data. Beijing quarantined people linked to Xinfadi market after June 14. Most of the new confirmed infection cases after June 14 came from isolated people who were unable to participate in the epidemic transmission process. Thus, before June 14, we set the infection period as the average infection period of daily new confirmed infection cases. After June 14, the infection period is set to 1 (removal rate = 1).

Nucleic acid testing intensity

Under the definition of infection period, we use the days required for nucleic acid testing (the time from participation in nucleic acid testing to the publication of testing results) to quantify nucleic acid testing intensity (“testing intensity”) in the transmission process [Equation (14), ]. The shorter the days required for nucleic acid testing, the higher the testing intensity, the shorter the infection period, and the higher the removal rate (Equation 15).

We used the average of the time (1.6 days) interval between the onset and the first visit of the infected person in the epidemiological data as the value of . is related to many unquantifiable factors (government organization capabilities, etc.), so we set this period to 0. According to Equation (14), we subtracted (1.6 days) from the infection period to obtain the daily testing intensity. In addition, we can adjust the testing intensity [Equation (14), ] to control the removal rate and build other scenarios.

Quantification of human mobility restrictions

Mobile phone signaling data from June 13 to June 28, 2020 were affected by human mobility restrictions. To obtain normal human mobility data not affected by interventions, we used mobile phone signaling data from May 9 to May 24, 2020 instead. Data for these two periods were collected from two Saturdays (May 9 and June 13, 2020) to two Sundays, 2 weeks later (May 24 and June 28, 2020), so that replacement did not affect the periodicity of the mobility data. Finally, we calculated the mobility reduction rate between the two data periods for each sub‐district in Beijing to quantify human mobility restrictions from June 13 to June 28, 2020 (Figure 3).

FIGURE 3

Quantification of human mobility restrictions: (a) human mobility reduction rate entering and leaving the sub‐districts; (b) intra‐sub‐district mobility reduction rate

Scenario establishment

Simulation scenario

Simulate the spatio‐temporal process of the outbreak in Xinfadi market, using basic parameters. In this scenario, the results of epidemic curve and spatio‐temporal distribution of infected persons were compared with the actual data to verify the model’s accuracy.

Non‐intervention scenario

Delete all intervention policies based on simulation scenarios. For the human mobility restrictions, we used mobile phone signaling data from May 9 to May 24, 2020. For the nucleic acid testing, we set testing intensity as 12.6 (Zhou et al., 2020). This scenario is used to analyze the spatio‐temporal process of epidemic transmission without intervention policies. By comparing with the simulation scenario, this scenario explains the effect of the intervention policies in Xinfadi market.

Human mobility and testing intensity adjustment scenario

Adjust the values of human mobility reduction rate and testing intensity in the non‐intervention scenario to quantify the impact of different intensity intervention policy combinations on epidemic control. Human mobility restrictions of different intensities (0, 20, 40, 60, and 80% reduction rates) were combined with nucleic acid testing (1, 2, 3, 4, and 5 days) to establish 25 policy combinations. The number and spatial distribution of infected persons under these combinations were analyzed by adjusting the reduction rate of human mobility and testing intensity according to the combinations (note: the other parameters are consistent with the baseline, i.e., the non‐intervention scenario). To avoid the impact of the scenario’s intervention policies on the model results, we chose the non‐intervention scenario as the baseline, which can reflect the slight policy adjustment.

Failed epidemiology survey scenario

Move the quantified human mobility restrictions (Figure 3) to Xicheng District, and testing intensity is the same as the simulation scenario. This scenario supposes the epidemiological survey mistakenly identified the outbreak source as Xicheng District and the government wrongly concentrated its intervention policies on Xicheng District. Compare the model results under this assumption with actual data to emphasize the importance of the correct epidemiological survey.

The fundamental difference between scenarios lies in the spatial scope of intervention policies. The spatial scope of the intervention policies for the human mobility and testing intensity adjustment scenario is citywide. In contrast, the simulation scenario is concentrated in Xinfadi market, the non‐intervention scenario has no scope, and the failed epidemiology survey scenario is concentrated in Xicheng District.

Model implementation

All analyses were done with Python, version 3.7 and ArcGIS Pro, version 2.3 (Esri, Redlands, CA).

RESULTS

Measures for non‐pharmaceutical interventions

Based on the COVID‐19 situation, Beijing’s government did not strictly reduce human mobility like the policies adopted by other Chinese municipal governments in January–March (Xi et al., 2020). Instead, it implemented non‐pharmaceutical interventions centered on community control, imposing different intensities of mobility restrictions and nucleic acid testing of 2019‐nCov on different sub‐districts according to the severity of the spread.

We used mobile phone signaling data to analyze human mobility from four perspectives (city, district, sub‐districts, and outbreak source). Beijing implemented human mobility restrictions on June 13, 2020. At the city level, the human mobility curve in June showed a downward trend after June 13, with a decrease of 10.36% (Figure 4a). At the district level after June 13, human mobility decreased by 19.76 and 14.99% in Fengtai and Daxing districts (Figure 4b), respectively, where the outbreak was severe (Figure 2), while human mobility decreased by only 6.58 and 4.75% in Huairou and Yanqing districts (Figure 4), respectively, where there were no confirmed infection cases (Figure 2). This indicates that Beijing’s human mobility restrictions were targeted. At the sub‐district level, mobility in Huaxiang, Huangcun, Qingyuan, Xincun, and Xihongmen sub‐districts decreased by 31.74, 20.55, 28.55, 28.52, and 29.03%, respectively (Figure 4c). These sub‐districts had the largest number of confirmed cases in Beijing (Figure 2); thus, their human mobility restriction intensities were higher than in other areas. In the area around Xinfadi market, human mobility decreased by 60.56% after June 13 (Figure 4d); even at the source of the outbreak, the area was not completely locked down. We selected the above case regions for comparative analysis based on the cumulative number of infections in different regions of Beijing.

FIGURE 4

Human mobility, daily new confirmed infection cases, and infection period in the outbreak of COVID‐19 linked to Beijing’s Xinfadi market: (a) city‐level daily population movement; (b) district‐level daily population movement; (c) sub‐district‐level daily population movement; (d) daily population movement in Xinfadi market; (e) daily new confirmed infection cases; and (f) average infection period of daily new confirmed infection cases. Human mobility refers to the number of people moving within the region each day, the definition and calculation method are given in Section 2.6. Infection period data of June 16, June 17, June 30, and July 1 in (f) are missing

Nucleic acid testing is also an essential part of the non‐pharmaceutical interventions for epidemic control because it can quickly identify the source of infections, shortening the infection period. Beijing isolated and tested all individuals linked to Xinfadi market on June 14, 2020 (Wu, Wang, et al., 2020). Based on the curve of daily new confirmed infection cases, the spread of the epidemic was controlled after June 14 (Figure 4e). The average infection period of daily new confirmed infection cases after June 14 was significantly lower (Figure 4f). (Data source: publicly available epidemiological data; Wu, Wang, et al., 2020. Infection period data missing on June 16, June 17, June 30, and July 1.) Considering that most of the new confirmed infection cases after June 14 came from populations already in isolation (Wu, Wang, et al., 2020), the actual infection period after June 14 should be much lower than this.

We may summarize the non‐pharmaceutical intervention of Beijing. It implemented limited human mobility restrictions based on the epidemic's severity in different areas. The nucleic acid test isolated high‐risk populations (people with abnormal body temperature or related to Xinfadi market) and stopped the epidemic with limited medical equipment. These policies reduce the risk of the epidemic spreading and do not significantly impact the daily lives of people in Beijing.

Simulation scenario

In the simulation scenario (Figure 5a), the coefficient of determination (R 2) between the simulated daily new infection curve and the actual daily new infection curve was 0.734, indicating that our simulation results fit the actual data well. In the simulation, the Beijing epidemic was predicted to peak on June 13, 2020, with a total of 329 (95% CI: 317–340) cases on June 28, 2020. In the simulation, the first infection case appeared at Xinfadi market on May 28 (Figure 5b). On June 13, the epidemic had begun to spread massively (Figure 5f). After June 13, Beijing implemented non‐pharmaceutical interventions (Figure 5g); subsequent new infections were caused by the onset of cases in the incubation period (Figure 5h). At the end of the simulation period (June 28), our model generated a total of 329 infection cases distributed across 35 sub‐districts (Figure 5i), whereas 335 confirmed infection cases were distributed across 45 sub‐districts of Beijing in the report on July 5, 2020 (Figure 2). The spatial distribution of the simulation results and the actual data are similar.

FIGURE 5

Simulation scenario: (a) simulated and realistic daily new infection curves; (b–i) distribution of cumulative infections in Beijing at different periods under the simulation scenario

Non‐intervention scenario

In the non‐intervention scenario simulation (Figure 6a), there were a total of 127,988 (95% CI: 108,115–147,862) cases on June 28, 2020, and the epidemic had not peaked at the end of the simulation period (June 28). In the spatial result, before June 7, the epidemic was only spreading through the sub‐districts around Xinfadi market (Figures 6b,c). By June 18, the epidemic had spread to most sub‐districts in Fengtai, Xicheng, and Dongcheng districts, and the number of infections already doubled simulated scenarios (Figure 6f). By June 23, the epidemic had spread to most sub‐districts within the Fifth Ring of Beijing (Figure 6f). On June 28, the epidemic had spread to 290 sub‐districts of Beijing (Figure 6h).

FIGURE 6

Non‐intervention scenario: (a) daily new infection curves of the non‐intervention scenario and reality; (b–h) distribution of cumulative infections in Beijing at different periods under the non‐intervention scenario

Human mobility and testing intensity adjustment scenario

We used human mobility reduction rate and testing intensity to quantify the non‐pharmaceutical interventions and explore their impact on controlling the outbreak. We adjusted the values of these two parameters based on the non‐intervention scenario and counted the cumulative infections of each result.

In the trend planes composed of cumulative infections (Figure 7a), three trend planes represented the upper limit, average value, and lower limit of the 95% CI of each simulation result. The maximum cumulative number of infections in the trend plane was 61,800 (95% CI: 50,200–73,300; 0%, 5 days); the minimum was 225 (95% CI: 183–267; 80%, 1 day). We adjusted the parameters in equal proportion (trend plane diagonal). For every 20% reduction in human mobility (20, 40, 60, and 80%) and 1‐day reduction in testing intensity, the cumulative number of infections fell by 61.33% (95% CI: 59.75–53.75%), 89.09% (95% CI: 88.88–89.09%), 97.85% (95% CI: 97.65–98.05%), and 99.64% (95% CI: 99.63–99.66%), respectively. The trend plane diagonal represents the general trend of the result; hence, we can conclude that the cumulative number of infections will decrease rapidly as the human mobility and the days required for the nucleic acid testing decrease.

FIGURE 7

Human mobility and testing intensity adjustment scenario: (a) trend plane composed of cumulative infections in each result; (b) section view for mobility reduction; and (c) section view for testing intensity

To quantify the impact of the two factors separately, we made section views based on the trend plane. For mobility reduction (Figure 7b), we set the testing intensity to the maximum value (5 days); for every 20% reduction in human mobility (20, 40, 60, and 80%), the cumulative number of infections fell by 31.23% (95% CI: 28.69–32.88%), 50.32% (95% CI: 49.93–50.80%), 59.55% (95% CI: 56.75–63.15%), and 67.96% (95% CI: 67.26–68.92%), respectively. For testing intensity (Figure 7c), the human mobility reduction rate was 0%, and with every 1‐day reduction in testing intensity, the number of infections fell by 49.19% (95% CI: 47.61–51.20%), 79.45% (95% CI: 78.85–80.32%), 95.57% (95% CI: 95.50–95.66%), and 99.36% (95% CI: 99.35–99.37%), respectively. There was an obvious tendency for the cumulative number of infections to decrease as the days required for the nucleic acid testing decreased when one parameter is set to the maximum (Figure 7c). Furthermore, compared to the other curves in the section views, fewer days required for nucleic acid testing led to a more moderate downward trend in the cumulative number of infections as human mobility decreased (Figure 7b). Therefore, the enhanced testing intensity was deemed more important than human mobility reduction for this outbreak.

The effect of these two parameters on spatial spread cannot be discovered from the line graphs alone. Therefore, we selected some combinations of parameters and present their spatial results.

When testing intensity is 5 days and the human mobility reduction rate is 0%, the epidemic spreads to most sub‐districts (Figure 8b). When testing intensity is unchanged, and the human mobility reduction rate is set to 80%, the epidemic still spreads to the central area of Beijing, but the spread range is significantly smaller (Figure 8c). When the testing intensity is 3 days and the human mobility reduction rate is 0%, the spatial extent of the epidemic is reduced to the Fengtai, Daxing, and Fangshan districts (Figure 8d). If the testing intensity is unchanged and the human mobility reduction rate is set to 80%, the epidemic spread is basically controlled in Fengtai District (Figure 8f). Mobility reduction reduces the scale of the epidemic to a greater degree than increasing testing intensity.

FIGURE 8

Spatial results of human mobility and testing intensity adjustment scenario: (a, e) daily new infection curves; and (b–d, f–h) the final distribution of cumulative infections in Beijing under different human mobility reduction rate and testing intensity

When testing intensity is set to 1 day, and the human mobility reduction rate is set to 0%, the number of sub‐districts to which the epidemic spreads and the number of people infected is very small (Figure 8e), but its spatial distribution is relatively scattered (Figure 8g). Suppose testing intensity is unchanged, and the human mobility reduction rate is set to 80%. In that case, the number of sub‐districts where the epidemic spreads and the number of people infected is still very small (Figure 8e), but its spatial distribution is relatively concentrated (Figure 8h), conducive to controlling the epidemic. However, considering the allocation of medical resources, it is practically impossible to achieve a 1‐day testing intensity. Human mobility reduction reduces the spatial spread, while enhanced testing intensity reduces the cumulative number of infections. Neither intervention alone can completely control the epidemic; a combination of both is required.

Failed epidemiology survey scenario

To verify the importance of accurately identifying the outbreak location, we established the failed epidemiology survey scenario. Xicheng District, where the first case occurred, was considered the outbreak site. We moved the mobility reduction policy (Figure 4a) to this district; testing intensity was the same as the simulated scenario. This error tripled the cumulative number of infections (961 [95% CI: 850–1,072]; Figure 9a). Xicheng District is located in the geographic center of Beijing; restricting mobility there would still allow spread to other sub‐districts within the Fifth Ring Road. Eventually, the outbreak would spread beyond the central city (Figures 9b–h).

FIGURE 9

Failed epidemiology survey scenario: (a) daily new infection curves of the failed epidemiology survey scenario and reality; and (b–h) distribution of cumulative infections in Beijing at different periods under the failed epidemiology survey scenario

DISCUSSION AND CONCLUSIONS

Our model was able to generate multiple human agents based on mobile phone signaling data to simulate the transmission process of the epidemic within Beijing. Although public health interventions vary among countries, most interventions are aimed at reducing crowd contact and infection period. We quantified policies associated with those two key factors: the human mobility reduction rate and testing intensity. Our model not only shows the epidemic curve, but also presents the spatial distribution of epidemic spread on a map; hence, it may be valuable for policymakers trying to assess the impact of various public health interventions on epidemic control, especially when seeking the optimal combination of interventions (e.g., intervention intensity in sub‐regions, key lockdown areas, policy implementation timing, etc.) and the magnitude of interventions (e.g., mobility restriction from 0 to 80% and testing intensity from 1 to 5 days). Our model can also weigh the potential benefits and negative effects of different policies. For example, increasing the testing intensity and reducing some mobility restrictions may ease the pressure of interventions on the regional economy.

Beijing’s interventions reduced the number of infections by 99.74% (99.69–99.77%). In the human mobility and testing intensity adjustment scenario, reducing human mobility decreased the cumulative number of infections by 49.59% (31.23–67.96%), while increasing testing intensity reduced the cumulative number of infections by 80.49% (61.33–99.64%). Human mobility reduction tends to reduce the spatial extent of the epidemic spread. Enhanced testing intensity significantly reduces the cumulative number of infections, but not the spatial extent of the epidemic spread. Therefore, we can infer from this case that a combination of the two measures can effectively control an outbreak like Beijing’s. Furthermore, erroneous investigation results significantly increase the cumulative number of infections.

The outbreaks in Beijing's Xinfadi market and Wuhan's Huanan seafood wholesale market were similar; both were accidental and unpredictable COVID‐19 outbreaks within the city. However, due to the different interventions of the Beijing and Wuhan municipal governments, the results of the two outbreaks were completely different. The Wuhan government failed to notice the high contagiousness of COVID‐19, resulting in an overwhelming number of infected person activities within the city, greatly increasing the scale of the epidemic. Despite the subsequent implementation of human mobility restrictions, there was a lack of equipment for nucleic acid testing, and the medical system collapsed due to insufficient medical resources and the large population in need of medical treatment. In contrast, Beijing discovered the outbreak site through an epidemiology survey after confirming the first infected case. Through targeted large‐scale testing and precise human mobility restrictions, Beijing's Xinfadi market COVID‐19 outbreak was stopped at an early stage. Therefore, when governments deal with similar outbreaks, they should identify the outbreak site early and use nucleic acid tests based on human mobility restrictions to reduce the number of infected persons active in the transmission process.

Because testing involves various factors such as medical personnel, medical equipment, and government organizational capacity, a shortfall in any of these factors will weaken the testing intensity and increase the days required for nucleic acid testing. However, it is difficult to perform large‐scale nucleic acid testing in all regions for cities where the current epidemic is severe due to resource constraints. We recommend that these countries use an epidemiology survey to differentiate the epidemic transmission risk level in different regions, implement different intensities of human mobility restriction according to risk level, and then focus medical resources on mass nucleic acid testing and isolate infected cases in high‐risk areas.

The mobile phone signaling data and POI data used in the epidemic model highlight the spatial heterogeneity of urban epidemic transmission. This spatial heterogeneity is transmission difference due to different population characteristics (e.g., population distribution, travel distance, active population, age, education, etc.) and urban facilities (e.g., markets, subway, workplaces, etc.). In this research, we used coefficients to link population characteristics and urban facilities to R0. But their relationship with the transmission characteristics of COVID‐19 is more complex, which requires further research. The mobile phone signaling data can provide detailed information to characterize behaviors and characteristics at the population level. However, these data also have limitations. Mobile phone users do not cover specific subgroups of the population, especially children and the elderly. Therefore, our population characteristics may also be an incomplete description of the urban population. POI data are easy to access, but public POI data do not have much available information; only the location and names of urban facilities are available. Different POI densities used in our research do not fully represent the impact of urban facilities on epidemic transmission. Access to key parameters of urban facilities such as area, ventilation conditions, and daily flow would improve our research.

This study has several limitations. First, the model is a probabilistic model based on an agent network, which is susceptible to real‐world random events, thus requiring multiple simulations to obtain relatively stable results. However, the model architecture is complex; therefore, multiple runs will take a long time. Second, we selected human mobility restrictions and nucleic acid testing for analysis. Other interventions, such as wearing face masks, ventilation indoors, and handwashing, can also control the outbreak. However, these interventions are not easy to quantify. Third, we set R0 to [1.5, 2.5] for each infectious individual in the model, without accounting for differences in infectivity between infectious individuals (e.g., asymptomatic, highly infectious individuals), but there were no relevant data for reference. Fourth, our model requires too many parameters. The main parameters include the transmission characteristics of COVID‐19 such as R0, infection period, simulation period, and incubation time, and changes in these parameters can significantly impact the results. However, the precise transmission characteristics of COVID‐19 are still unclear and require further research. In this research, the central parameters come from the epidemiological survey report, and its quality may be affected by the capacity for diagnosis and surveillance, regions, and time. Other parameters include the weight coefficients of some input data, such as POI density and human mobility.

On the one hand, these parameters can help our model establish many scenarios; on the other hand, the tuning of the model is highly complicated. Human mobility and a population estimated from mobile phone signal data are essential for our research, and we purchased the data from the service provider (China Mobile Ltd). However, it is unclear whether China Mobile’s Beijing market share (more than 60%) and base station distribution will affect human mobility and population estimates. Additionally, our use of high‐quality mobile phone signaling data may limit model promotion.

Our future research will aim to improve the estimation method of the model parameters, reduce the sensitivity of the model to random events, and reduce the model run time.

CONFLICT OF INTEREST

The authors declare no competing interests.