ABSTRACT
Background As COVID-19 continues to spread around the world, understanding how patterns of human mobility and connectivity affect outbreak dynamics, especially before outbreaks establish locally, is critical for informing response efforts.
Methods Here, in collaboration with Facebook Data for Good, we built metapopulation models that incorporate human movement data with the goals of identifying the high risk areas of disease spread and assessing the potential effects of local travel restrictions in Taiwan. We compared the impact of intracity vs. intercity travel restrictions on both the total number of infections and the speed of outbreak spread and developed an interactive application that allows users to vary inputs and assumptions.
Findings We found that intracity travel reductions have a higher impact on overall infection numbers than intercity travel reductions, while intercity travel reductions can narrow the scope of the outbreak and help target resources. We also identified the most highly connected areas that may serve as sources of importation during an outbreak. The timing, duration, and level of travel reduction together determine the impact of travel reductions on the number of infections, and multiple combinations of these can result in similar impact.
Interpretation In Taiwan, most cases to date were imported or linked to imported cases. To prepare for the potential spread within Taiwan, we utilized Facebook’s aggregated and anonymized movement and colocation data to identify cities with higher risk of infection and regional importation. Both intracity and intercity movement affect outbreak dynamics, with the former having more of an impact on the total numbers of cases and the latter impacting geographic scope. These findings have important implications for guiding future policies for travel restrictions during outbreaks in Taiwan.
Funding Ministry of Science and Technology in Taiwan and National Institute of General Medical Sciences in USA
INTRODUCTION
The Coronavirus Disease 2019 (COVID-19) was first reported in Wuhan, China in December 2019 and has since caused a global pandemic, with over 500,000 confirmed cases and over 23,000 deaths reported by March 27, 2020.1 Scientific discoveries have advanced at an unprecedented pace, with numerous clinical trials of drugs underway,2 and a phase 1 vaccine trial beginning in Seattle.3 Despite these advances, no pharmaceutical interventions are yet available. In the meantime, public health officials must rely on other interventions, such as social distancing and travel restrictions, to slow the spread and reduce the peak of the outbreak, in order to prevent health systems from being overwhelmed.4,5
In January 2020, as the epidemic in Wuhan grew, many countries implemented travel bans, and airlines canceled flights to attempt to slow the spread.6 A number of studies have estimated the risk of importation globally, with some suggesting up to two-thirds of all imported cases went undetected.7,8 For Taiwan, there have been 267 reported cases as of March 27, 20209, with 227 imported (85%) and 40 local cases (15%). 37 local cases (92.5 %) were linked to imported or known cases, and 3 local cases (7.5 %) have unknown origin.
As the number of cases globally due to community transmission grows relative to the number of imported cases, attention has turned to more local measures to decrease spread, such as cancellations of mass gatherings, business closures, and local travel restrictions.10 Mobility data can provide critical information for responding to outbreaks and understanding the impact of travel restrictions.11 Recent studies have analyzed the effects of human mobility and travel restrictions on disease spread in China.12,13 Here, to prepare for COVID-19 and its impact, in collaboration with Facebook Data for Good, we describe the metapopulation models we’ve built that include human movement data to better understand the high risk areas of disease spread and assess the potential impact of local travel restrictions in Taiwan.
METHODS
Movement data
We incorporated two different sources of movement data from Facebook into our models. In the first model, we used Facebook’s newly developed colocation matrices (Facebook colocation data), which give the probability that people from two different geographic units will be in the same 600 m × 600 m location for five minutes using data over the course of a week. In the second model, we used Facebook’s regular movement data (Facebook movement data), which aggregates the number of trips Facebook users make between locations over time (Figure S1).14 Facebook movement data between January 26th and March 9th were used and were disaggregated by Lunar New Year (January 26th to January 29th), weekdays (Monday to Friday between January 30th and March 9th, except February 28th, a National Public Holiday in Taiwan), and weekends (Saturdays and Sundays between January 30th and March 9th plus February 28th).
Colocation covered five different weeks, with one including the Lunar New Year (January 10th to 16th) and the other four including only regular days (February 16th to 22th, February 26th to March 3rd, March 4th to March 10th, and March 11th to March 17th). As the colocation data from the four regular weeks were very similar, we used data from February 16th to 22th for our model. The geographic unit used in this study was at the centrally-governed level of “city” (here “city” indicates city, county or special municipality in Taiwan). Shape files were downloaded from Government open data platform (https://data.gov.tw/dataset/7442). We excluded three cities outside of the main island of Taiwan from the analysis due to their low connectivity with the main island, leaving 19 cities.
Models
We developed susceptible-exposed-infectious-recovered (SEIR) models of the spread of COVID-19 throughout Taiwan. Note that this model does not attempt to fit epidemiological dynamics to data on COVID-19 cases, and we do not account for the many complexities associated with asymptomatic transmission and different clinical outcomes. Rather, we aim to examine the spatial dynamics that would result from transmission of SARS-COV-2 across Taiwan under various mobility scenarios. As such, these dynamics may be generalizable to other pathogens with similar transmission routes.
We ran a deterministic SEIR model to analyze the overall infection dynamics, including time to peak and peak size, and we also ran the models stochastically to better understand the initial stages of disease spread. We ran the stochastic model until (1) it reached n accumulated infections or (2) the total number of infections became 0 to estimate the probability of having more than n infections (denoted by Pn,k, where k represents the number of initial infections), the time it takes to reach n infections (denoted by Tn,k), and the standard deviation of infection numbers at Tn,k (denoted by Vn,k). To assess the initial stages of the outbreak, we used n=1000 and k=3 as our baseline values.
Let Si, Ei, Ii, Ri be the number of susceptible, exposed, infectious, and recovered individuals in location i, respectively, and Ni be the total population in location i. Let R0 (=2.5)15 be the average number of secondary infections from one infected individual, DE (=5)16 be the incubation period (i.e. time from exposure to symptom onset), and DI (=7)17 be the duration of infectiousness.
We incorporated the Facebook data into these models in two ways. In the “contact model”, we used the Facebook colocation data (Cij, the probability that a person from location i collocates with a person from location j) and assumed that the contact rate between a random pair of people was proportional to Cij. We scaled R0 by Cij*Nj (for j not equal to i) or Cii*(Ni−1), standardized to Cii*(Ni−1) in Taipei, which was assumed to be the same as Wuhan.
We defined R0ii as intracity R0 and Σj ≠ i R0ij as intercity R0 for location i. The sum of intracity R0 and intercity R0 reflects total risk of infection and was standardized to the highest sum of intracity R0 and intercity R0 as follows:
Facebook colocation data from regular days were used to calculate risk of infection.
In the “residence model”, we first estimated the proportion of time people living in location i spend in location j (Pij) based on Facebook movement data (see details in Supplementary Methods), and modeled the transmission dynamics by considering both that (1) non-travelers get infected by infectious visitors to their home location (the first part in the following equation) and that (2) naïve travelers get infected when they travel (the second part in the following equation).
The remaining equations are the same across the two models.
We defined 
as risk of regional importation (i.e. importation from other cities within Taiwan) for location i, where qj represents the average number of subscribers in location j and mji represents the average number of people moving from location j to location i per unit of time in Facebook movement data. Source of importation was defined as the number of travelers from each location i and standardized to the highest value.
Weekday movement data were used to calculate risk of regional importation and source of importation. Risk of regional importation, source of importation, and risk of infection defined above were calculated from movement data and represent different connectivity measures relevant for disease transmission.
One of the major differences between the contact model and the residence model is that the transmission rate within each city (R0/DI) varies with colocation matrices in the contact model, while it remains constant in the residence model. In this sense, the contact model is similar to the traditional density dependent model, where contact rates (and therefore transmission rates) vary with population density, and the residence model is similar to the frequency dependent model.18 As it is unclear which is most appropriate for COVID-19, we used both and compared the results.
Modeling travel reduction
To assess the potential effect of travel restrictions at multiple levels, we modeled either intra-city travel reductions, inter-city travel reductions, or a combination of both travel reductions (“overall reduction” in texts and figures) for 1, 2, 3, or 6 months or for the whole period of time. Travel reductions started from the beginning of the simulations or when there were 10, 20, 30, 50, or 100 accumulated infections. The proportion of reduction is denoted by G. In the contact model, intracity reduction was modeled by R0ii*(1–G) for all i, and intercity reduction was modeled by R0ij*(1–G) for all i not equal to j. In the residence model, intracity reduction was modeled by R0*(1–G) and intercity reduction was modeled by Pij*(1–G) for all i not equal to j and Pii +(1– Pii)*G for all i.
RESULTS
Initial stage of transmission
At the end of March 2020, most cases in Taiwan were imported or linked to imported cases. Therefore, we explored a variety of initial conditions when estimating the probability of having more than 1000 infections (denoted by P1000,k, where k represents the number of initial infections), varying both the number of initial infections and their locations. As expected, we found that, if we assumed that the transmission rates varied among cities (contact model), the probability of having more than 1000 infections also varied, with the cities with larger intracity and intercity R0 (and larger risk of infection) showing larger P1000 (Figure 1, Figure S2 and Table S1). The top 5 cities with high risk of infection due to larger contact rates included Taipei City, New Taipei City, Kaohsiung City, Keelung City, Hsinchu City (Figure 2 and Table S2).
P1000,3 from simulations with initial infections in Taipei City (representing big cities) or Changhua County (representing small cities) using both contact and residence models are shown. The difference between big and small cities was more significant in the contact model than in the residence model. Intracity and intercity travel reduction reduced P1000,3, while the impact of intercity travel reduction was minor. Increased human mobility during Lunar New Year increased P1000,3.
Three kinds of connectivity measures relevant to disease spread are shown. The values for bigger cities were larger. (A) Risk of infection. (B) Risk of importation. (C) Source of importation.
In simulations where 1000 infections were reached, the time it took to reach 1000 infections (denoted by T1000,k) was also shorter for cities with larger intracity and intercity R0 (Figure S3). Due to higher mobility surrounding the holiday season, P1000,3 was higher (Figure 3) and T1000,3 was lower during the Lunar New Year than regular days. Because of the high variation in transmission rates in the contact model, the majority of infections at T1000,3 were in cities with larger intracity R0 no matter where the initial infections were (Figure S4).
P1000,3 was significantly higher if human mobility during Lunar New Year was used (Wilcoxon signed rank test, p-value= 3.82×10−6). Each point represents a city in Taiwan, and the red line is the x=y line.
When assuming that the transmission rates in different cities were the same (residence model), the probability of having more than 1000 infections and the time to reach 1000 infections did not vary much among cities (Figure 1, Figure S3 and Table S3). P1000,3 during the Lunar New Year were not significantly different from regular weekends or weekdays either (Wilcoxon signed rank test, p-value= 0.57 [weekends] and 0.72 [weekdays]). The majority of infections at T1000,3 were in cities with initial infections and the cities well-connected to them (Figure S4; and more details at https://roachchang.shinyapps.io/TW_CoV_Dynamics/). The effect of intercity connectivity, however, was reflected in the variation in infection numbers across cities at T1000 (denoted by V1000; Figure S5). The variation in infection numbers was lower in more connected cities as the chance of spreading the virus to other cities was higher.
In both models, well connected cities played more important roles, as they spread the virus to other cities more quickly and more widely. The top 5 cities with the highest risk of regional importation included Taipei City, Hsinchu City, Chiayi City, New Taipei City, and Hsinchu County, and the top 5 cities with the highest source values of importation were New Taipei City, Taipei City, Taoyuan City, Hsinchu City and Hsinchu County (Figure 2 and Table S2). In the contact model, these places also had higher intracity and intercity R0 values (i.e. bigger cities are more densely populated and more well connected) (Table S2), further increasing the speed of the spread. Similarly, higher travel volume during the Lunar New Year increased the spread of the virus at the initial stage compared to regular days in both models.
The effect of travel restrictions at initial stage
We analyzed Facebook’s movement data and found that intercity movement did not decrease between January 30th and March 9th in Taiwan among Facebook users (Figure S6A), suggesting that without travel restrictions imposed by the government, the public in Taiwan may not have reduced their frequency of travel voluntarily due to the awareness of SARS-CoV-2. We then examined the level to which travel restrictions could reduce the spread of SARS-CoV-2 in Taiwan at the initial stage of an outbreak. In both the contact model and the residence model, decreasing intracity movement had a much larger impact on P1000,3 (Figure 1) and T1000,3 (Figure S3) than decreasing intercity movement. The impact of reducing intercity travel was most evident in influencing how widespread the virus was: the infections were located in only a few cities at T1000,3 if intercity travel was reduced (Figure S5).
We also examined the impact of different lengths of travel reductions (Figure 4). The probability of having more than 1000 infections decreased when the duration of intracity travel reduction increased in both models, while it did not change with the duration of intercity travel reduction. The results also suggest that higher levels of reduction and longer periods of reduction for intracity travel can have similar impacts. For example, a 60% intracity travel reduction for 2 months had similar outcomes as a 70% reduction for 1 month. While P1000,3 did not change with the length of intercity travel reduction, longer intercity travel reduction led to slower progression of the outbreak (higher T1000,3) in the contact model and more clustered infections (higher V1000,3) in both models (Figure S7). Furthermore, we found that the best timing for the travel reduction depended on both the duration of travel reduction and the level of reduction (Figure S8).
P1000,3 from the contact model (A) and the residence model (B) with initial infections in Taipei City. The color represents the level of reduction in P1000,3 (white to red represents smaller to larger reduction).). As the duration of intracity travel reduction increased, P1000,3 decreased in both models. P1000,3 did not change with the duration of intercity travel reduction.
Overall dynamics
To assess the impact of travel on the overall dynamics, we performed deterministic simulations. Although the location of initial infections influenced the overall dynamics, with well-connected cities like Taipei City leading to higher total infection numbers and shorter time to peak than less-connected cities like Changhua County (Figure S9), the difference was minimal and therefore we only showed the results from Taipei (Figure 5 and Figure S10). Intracity travel reduction reduced the total number of infections and delayed time to peak in both the contact and residence models (Figure 5 and Figure S10), while intercity travel reduction had only minimal impact on the overall dynamics in both models. Similarly, increased intercity and intracity colocation probabilities during the Lunar New Year (Figure S6B) led to increased total number of infections and decreased time to peak in the contact model, but increasing only intercity movements during Lunar New Year in the residence model only shortened time to the peak minimally (Figure 5 and Figure S10), suggesting reducing large gatherings during holiday times is more important than reducing intercity travel.
Total infection ratio represents the ratio of total infections with travel reduction to total infections without travel reduction. The ratios shown here are from deterministic simulations with initial infections in Taipei City. Intracity and overall travel reduction effectively decreased the total number of infections. The impact of intercity travel reduction was small. For some conditions, number of infections decreased from the beginning of simulations and no bar was shown.
DISCUSSION
By utilizing aggregated human movement data from Facebook, we built metapopulation models to understand the potential spread of SARS-CoV-2 in Taiwan given different initial conditions and to assess the potential impact of travel restrictions. We identified the top cities with the highest risk of infection as well as the top cities with the highest importation risk from other cities based on Facebook data and population sizes. We made a web-based interface showing the geographic distribution of infections at different time points (T100, T500 and T1000) in the initial stages of the outbreak given different locations of initial infections. This will help the public avoid high-risk areas, help public health professionals identify surveillance targets, and inform decisions on travel restrictions.
Consistent with previous findings showing that international travel bans are less effective than social distancing within the country,6,19 we found that intracity travel reduction has a higher impact on disease dynamics than intercity travel reduction and increasing the length of intracity travel reduction increases the impact. Intercity travel reduction, however, influences the variation in infection numbers across cities and can reduce the number of cities that have infections at the initial stage of the outbreak. While intercity travel did not decrease the probability of having high numbers of infections, containing the infections to a few cities has important public health impacts, as this means surveillance system can focus on fewer cities and control efforts can be more targeted.
Intracity travel reduction in our model is effectively the same as any measure that can reduce contact rates between individuals, such as social distancing, or transmission probability given contact, such as hand washing or wearing facemasks. These measures have been shown to be effective in reducing the transmission of respiratory viral pathogens in both modeling and empirical studies,20–24 and should be encouraged. Our study found that similar overall numbers of infections can occur with various combinations of length of reductions and percent reduction in intracity travel (e.g. 60% reduction for 2 months had similar outcomes as 70% reduction for 1 month). The effectiveness of travel reductions varied depending on how early into the outbreak they started and their duration; for example, if restrictions lasted at least three months, starting them as early as possible had the greatest impact on P1000,3, but if limited to only one month, starting them later sometimes had more of an impact and sometimes did not, depending on the context. Health officials can therefore take into consideration feasibility of different interventions, impact on society, and the capacity of the healthcare system to determine the optimal interventions and their duration.5
We found that the volume of intercity travel around Lunar New Year, where both intercity movements and intracity colocation rates increased, can significantly increase the speed of virus spread. Another traditional holiday, Ching Ming festival, is coming soon (April 2 to 5), and the volume of human movement between cities during this period is expected to increase like Lunar New Year, but at a lower level. Our results suggest that it is important to avoid travel or reduce the impact of travel through measures such as limiting social interactions and wearing facemasks when taking public transportation to reduce the spread of the virus.
Although we did not observe a significant change in intercity movements from January 30th to March 8th, Facebook data have the potential to track how the volume and pattern of travel change through time as the outbreak progresses, and we can incorporate any change in human mobility into the metapopulation models in nearly real time to help fight COVID-19.11,25 Moreover, our model utilizing human mobility data from Facebook is not limited to intercity or intracity level, or Taiwan. Facebook movement data are also calculated at finer geographic scales (such as towns) and for other countries, and our model can be easily applied in these settings to understand disease dynamics of COVID-19.
Data Availability
Modeling output data are in the manuscript or can be accessed through the publicly available web-based interface (https://roachchang.shinyapps.io/TW_CoV_Dynamics/). Due to the Data Sharing Agreement with Facebook, readers can not access the original Facebook data used in this study.
SUPPLEMENTARY MATERIALS
SUPPLEMENTARY METHODS
Estimating Pij
We built a travel model to estimate the proportion of time people living in location i spend in location j (Pij) by fitting the model to the Facebook movement data. Xij represents the proportion of people living in location i currently in location j, and 
represents the equilibrium state of Xij, and its value under the fitted model is used as our estimate of Pij. People living in location i travel with probability Fi, and the probability that a traveler from location i travels to location j is denoted by Tij. Travelers go back to their home location at probability λi per unit of time. Mij,t,t+1 represents the number of people moving from location i to location j between time t and t+1.
At equilibrium,
For simplicity, we assumed that the majority of travel is work-related travel and on average travelers spend eight hours in the travel destination (λi =1 given the unit of time is 8 hours) and that Tij is proportional to Mij, leaving Fi the only parameters to be fitted. We used a gradient descent algorithm to find the local optimum solution for Fi, where the cost function is defined by the sum of the squared difference between normalized mij and the normalized value of Mij from the model. We calculated 
under fitted parameters to obtain estimates of Pij.
SUPPLEMENTARY FIGURES
(A) Regular movement data. (B) Colocation matrices.
Both intracity and intercity R0 were positively correlated with P1000,3 (Pearson’s correlation test, r= 0.93 [intracity], 0.83 [intercity], and 0.94 [risk of infection], p-value= 1.28×10−8, 1.34×10−5 [intercity], and 4.03×10−9 [risk of infection]). Intercity and intracity R0 were also positively correlated (Pearson’s correlation test, r= 0.80, p-value= 3.75×10−5, not shown).
If initial infections were in a big city, it took less time to reach 1000 infections in the contact model. The difference between big and small cities was not significant in the residence model. Intracity and overall travel reduction delayed the time to reach 1000 infections in both models, while intercity reduction did not. For some conditions, P1000,3 was 0 and no bar was shown.
The location of initial three infections is labeled with a white border. In the contact model, the first 1000 infections were in big cities; in the residence model, the first 1000 infections were located in cities with initial infections and the cities well-connected to them.
Intercity travel reduction increased the variation in infection numbers across cities in both models.
Lunar New Year ended on January 29th. The frequency of travel did not decrease after Jan 30th. Three travel routes are shown here, including a small city to big city (Changhua County to Taipei City), a big city to big city (Taipei City to New Taipei City) and a big city to small city (Taipei City to Yilan County).
T1000,3 (upper panel) and V1000,3 (lower panel).
P1000,3 when travel reduction started from the beginning of the simulations (denoted by 0), or when there were 10, 20, 30, 50, and 100 infections in both contact (left) and residence (right) models. Two different lengths of travel reduction duration were shown: (A) 1 month (B) 3 months. Only intracity travel reduction was shown here because intercity travel reduction only had minimal impact on P1000,3 and the results from overall reduction and intracity reduction were qualitatively similar. The timings associated with the lowest P1000,3 differed for one-month and three-month travel reductions. It was best to reduce travel at the beginning if the duration was for three months, but the pattern for one-month travel reduction was more complicated.
Total infection ratio and the ratio of time to peak are the infection numbers and the time to peak standardized to the highest values, respectively. The total infection numbers were similar between the outbreaks seeded by infections from big and small cities. The time to peak was slightly longer if the outbreak was initiated with infections in a small city.
The time to peak was standardized to regular days in the contact model and to weekdays in the residence model. The time to peak increased significantly with intracity and overall travel reduction, but only increased minimally with intercity reduction.
SUPPLEMENTARY TABLES
Colocation matrices in regular days were used.
Movement data on weekdays were used.
Acknowledgments
We thank Brian Karrer, Alex Pompe, Laura McGorman for their contributions.
Footnotes
↵* shared first authorship
