The temporal spreading pattern
We display the difference between the overall spreading centroid and cumulative spreading centroid until ith period in Fig. 1a, c, and e for Wuhan, Beijing, and Urumqi, respectively. Interestingly, as the situation of COVID-19 evolves, cumulative spreading centroids at different periods in Wuhan, Beijing, and Urumqi are close to the overall spreading centroid: the mean absolute errors (MAE) between cumulative and overall spreading centroids in Wuhan, Beijing, and Urumqi are 0.3 Km, 0.4 Km, and 0.7 Km, respectively. It is clear to see that the temporal spreading centroid of COVID-19 has a feature of time invariance. That is, the spreading centroid is stable and nearly does not migrate during COVID-19 spreading.
a–f The cumulative spreading centroid and spreading radius in Wuhan, Beijing, and Urumqi, respectively. g Relation between the mean travel distance of people in each city and the corresponding COVID-19 growth rate.
We also conduct analysis for the other six cities with a relatively large number of confirmed cases in China: Xiaogan, Suizhou, Xiangyang, Huanggang, Guangzhou, and Wenzhou (Supplementary Fig. 6), where cases are mainly imported. Similar conclusions for spreading centroid and spreading radius can be made in these cities. We proceed to perform sensitivity analysis by varying the number of spreading period L (Supplementary Fig. 5) and find that L does not have much impact on the observed temporal pattern. Therefore, the temporal spreading pattern of COVID-19 in China features time invariance of spreading centroid and slow growth of spreading radius.
As can be seen in Fig. 1 and Supplementary Fig. 4, there are significant disparities in growth rate of spreading radius in different cities and different time periods. To find intrinsic mechanisms for these disparities, we first divide the spreading period T (T in Beijing for instance lasted from June 11 to July 10 with 30 days) of each city into two periods L1 and L2 (L1 in Beijing for instance lasted from June 11 to June 25 with 15 days). Then, two spreading radii (R1 and R2) can be calculated based on activity centroids of confirmed cases reported in the spreading periods L1 and L2, respectively. We define the growth rate of spreading radius as 2(R2−R1)/∣T∣, where ∣T∣ denotes number of days in spreading period T. Further, we randomly select 20,000 smartphone users in each city and leverage their trajectory data during the outbreak of COVID-19 to compute their mean travel distance. Since these smartphone users are randomly selected, we use this mean distance (over 20,000 users) to approximate that of all citizens in each city. Clearly, a large value of mean travel distance reflects a strong willingness of people for long-distance travelling. Considering that different control measures imposed in Wuhan and Urumqi since the outbreak of COVID-19 affected the corresponding mobility pattern and spreading of pandemic significantly (Fig. 2d, e), we divide the spreading period of these two cities into two sub-periods: before and after the implementation of travel restriction, and then calculate mean travel distance and growth rate of spreading radius in these two sub-periods, respectively. The correlation analysis results for all nine cities are illustrated in Fig. 1g. Interestingly, we observe a clear positive correlation between mean travel distance and growth rate of spreading radius, indicating that mobility pattern accelerates the urban spreading of COVID-19.
a–c A visualization of the number of confirmed cases in discretized grids in Wuhan, Beijing, and Urumqi, respectively. d–f The human mobility distribution (dots) as a function of traveling distance and the fitted models (lines) for Wuhan, Beijing, and Urumqi, respectively. g–i The spatial distributions (dots) as a function of distance from the overall spreading centroid and the fitted regression lines for these distributions.
The spatial spreading pattern
To characterize and visualize spatial spreading pattern, we divide the geographical area into grids of 1 Km × 1 Km26. The overall spreading centroid is set as the original point of grids. Confirmed cases of COVID-19 are then projected into grids according to their activity centroids. As illustrated in Fig. 2a–c, three-dimensional histograms are used to describe the spatial distributions of confirmed cases in Wuhan, Beijing, and Urumqi, respectively, where the height of each bar represents the case count in each grid.
To analyze the spatial distribution function F(d), driven by human mobility pattern, we apply the logarithm to the actual distribution of human mobility pattern as well as the number of confirmed cases and the distance from the overall spreading centroid in all cities. The human mobility distributions and spatial distributions F(d) of Beijing, Urumqi, Xiaogan, Suizhou, and Huanggang exhibit a prominent linear pattern (Fig. 2d–i, Supplementary Figs. 6 and 7). Although a significant change in human mobility patterns in Wuhan (Jan. 23) and Urumqi (Jul. 16) after control measures can be observed intuitively, corresponding human mobility distributions and spatial distributions in these cities are surprisingly demonstrated as power-law model. Therefore, we adopt F(d) = dα for linear regression. Specifically, we have α = −1.80 for Beijing (the Pearson correlation: −0.93) and α = −2.15 for Urumqi (the Pearson correlation: −0.93), respectively.
As shown in Fig. 2g, the spatial distribution of confirmed cases in Wuhan is power-law-like since it deviates slightly from power-law model when d is small. Due to the influence of human mobility patterns in Wuhan during the lockdown period (Fig. 2d), initial cases have a higher probability to infect susceptible individuals around the spreading centroid. As a result, the risk of infection around spreading centroid is much higher than that at distant locations (a significant 92% of cases are close to the spreading centroid). To have a more accurate quantification of the spatial distribution of Wuhan, we divide the area into two parts by their distance to the spreading centroid, and adopt two different models to fit the data. Specifically, spatial distribution of confirmed cases around the spreading centroid (d ≤ 18) is fitted by an exponential model F(d) = αd, and when d > 18 is fitted by a power-law model. As illustrated in Fig. 2g, F(d), d ≤ 18 is well fitted by an exponential model with a Pearson correlation of −0.99, and F(d), d > 18 is well fitted by a power-law model with a Pearson correlation of −0.96. This indicates that the spatial distribution of confirmed cases in Wuhan indeed has different characteristics, depending on the distance to the spreading centroid. We also observe a similar phenomenon in Xiangyang (Supplementary Fig. 7f). Interestingly, we find that this is mainly determined by human mobility pattern (to be elaborated on in the next Section).
We also notice that for cities (such as Guangzhou and Wenzhou) where imported cases are widely scattered, the spatial spreading pattern is less prominent. It is clear that there are multiple clusters of confirmed cases in these two cities (Supplementary Fig. 8), which impacts the power-law-like spatial spreading. Therefore, the observed spreading pattern does not apply to the case with multiple infection sources.
The underlying mechanism
The classic susceptible-infected-recovered model (SIR) and its variants have been widely adopted to understand the transmission characteristics of infectious diseases. The Kendall model27,28,29 introduces the spatial dimension to the SIR model and can be used to explain the spatial-temporal evolution of infectious diseases. Its differential equations can be expressed as equations (1)–(3) of supplementary materials. Note that confirmed cases in China get isolated for medical treatment once they are confirmed and would not cause further infection. Under such a condition, the confirmed cases can be regarded as recovered individuals in the Kendall model. Then, the differential equation for the proportion of recovered individuals can be written as
$$\frac{\partial R}{\partial t}=-\lambda R(x,t)+\lambda {I}_{0}(x)+\lambda \left[1-\exp \left(-\frac{1}{\lambda }\int\nolimits_{-\infty }^{\infty }R(y,t)K(x-y){{{\rm{d}}}}y\right)\right].$$
(1)
where R(x, t) denotes the proportion of recovered individuals at location x and time t, satisfying R(x, 0) = 0. Note that λ can be obtained by inverse of the basic regeneration number R0 in the model, that is, λ = 1/R0 = γ/βξ, where γ, β, ξ represents recovery rate, infection rate and the number of initial susceptible individuals, respectively. Besides, the kernel function K(x − y) > 0, satisfying \(\int\nolimits_{-\infty }^{\infty }K(y){{{\rm{d}}}}y=1\), quantifies the probability that an infected individual at location y visits x. Here, we use power-law distribution to describe the city-level movement behaviors30, which can be written as K(Δr) = Δrη. Hereby, K(Δr) represents the probability for the step size Δr and η, the power-law exponential, denotes the travel willingness, which has a strong correlation with the mean traveling distance.
To fit the model for recovered individuals, we first calculate the parameter η in the power-law distribution by utilizing the mobility data of anonymous smartphone users, with which to capture the inherent human movement behaviors for each city. Moreover, the diagnosed date for each confirmed case and corresponding activity centroid is also calculated as input of the model. Through fitting the model (i.e., equation (1)) based on these precalculated parameters and Least Squares algorithm, we can finally obtain a set of optimal parameters (λ and I0(x)). Note that we assume initial confirmed cases originate from the grid (0, 0). The parameter I0(x) fitted in the model could therefore be written as I0(0, 0). A detailed discussion about Kendall model and parameter fitting process are provided in supplementary materials.
Figure 3a–c illustrate the overall model fitting performance of recovered individuals R(x, t) for Wuhan, Beijing, and Urumqi, in which Root-Mean-Square-Error (RMSE) is also added to quantify the performance of model fitting. The values of RMSE for Wuhan, Beijing, and Urumqi during whole spreading period are 0.55, 0.18 and 0.57 respectively, indicating that the evolution of recovered individuals R(x, t) during spreading period can be well captured by the proposed Kendall model. More results on the performance of model fitting can be found in Supplementary Fig. 9.
a–c The overall model fitting performance. d Relation between travel willingness η and Simpson divergence. e, h Relation between travel willingness η and daily reported confirmed cases. g Relation between travel willingness η and the growth rate of spreading radius. f, i Relation between travel willingness η and spreading radius.
We proceed to study the impact of parameter η in the model on spatial dispersion of confirmed cases, the number of daily new confirmed cases, and growth rate of spreading radius during the whole spreading period. To characterize the spatial dispersion of confirmed cases, we introduce the concept of Simpson Divergence: \(div=1/{{{\Sigma }}}_{i}{p}_{i}^{2}\), where pi represents the proportion of confirmed cases distributed in grids whose distance to overall spreading centroid is within [i, i + 1]. Therefore, a small value of div reflects a high clustering of the confirmed cases, i.e., a large proportion of confirmed cases are distributed in a small number of grids. Fig. 3d illustrates the impact of η on spatial dispersion of confirmed cases. We consider two scenarios under different basic regeneration number R0: R0 < 1 and R0 > 1. We see that, with the decrease of η Simpson Divergence decreases to 1, at which all confirmed cases are distributed in one grid. The impact of η on growth rate of spreading radius is illustrated in Fig. 3g. Clearly, the growth rate of spreading radius decreases with the decrease of η, which is consistent with results in Fig. 3d. Figure 3e and h show the impact of η on the number of daily reported confirmed cases. Interestingly, a large η, which means long mean travelling distance in the human mobility model, results in quick spreading of COVID-19, which makes the peak of daily reported cases arrives early. This also shows that travel restriction policies will delay the peak arrival. Finally, we study the impact of η on the growth rate of spreading radius. As we can see in Fig. 3f, i, spreading radii under different η increase with time and converge to a fix value when R0 > 1. However, when R0 < 1, the spreading radius increases only when η is relatively large. This indicates that when R0 < 1 and the mean travelling distance is also low, the pandemic will not spread spatially. Therefore, η in the mobility model drives the temporal-spatial spreading process. In practice, we can optimize η by implementing a specific travel restriction policy to have a desired control and prevention performance.
