key: cord-0260034-t6efqz5h
authors: Berke, A. A.; Doorley, R.; Alonso, L.; Pons, M.; Arroyo, V.; Larson, K.
title: Using mobile phone data to estimate dynamic population changes and improve the understanding of a pandemic: A case study in Andorra
date: 2021-11-11
journal: nan
DOI: 10.1101/2021.11.06.21265955
sha: 1155db1eb70534864a515cb277e9deefa1a4de97
doc_id: 260034
cord_uid: t6efqz5h

Compartmental models are often used to understand and predict the progression of an infectious disease such as COVID-19. The most basic of these models consider the total population of a region to be closed. Many incorporate human mobility into their transmission dynamics, usually based on static and aggregated data. However, mobility can change dramatically during a global pandemic as seen with COVID-19, making static data unsuitable. Recently, large mobility datasets derived from mobile devices have been used, along with COVID-19 infections data, to better understand the relationship between mobility and COVID-19. However, studies to date have relied on data that represent only a fraction of their target populations, and the data from mobile devices have been used for measuring mobility within the study region, without considering changes to the population as people enter and leave the region. This work presents a unique case study in Andorra, with comprehensive datasets that include telecoms data covering 100% of mobile subscribers in the country, and results from a serology testing program that more than 90% of the population voluntarily participated in. We use the telecoms data to both measure mobility within the country and to provide a real-time census of people entering, leaving and remaining in the country. We develop multiple SEIR (compartmental) models parameterized on these metrics and show how dynamic population metrics can improve the models. We find that total daily trips did not have predictive value in the SEIR models while country entrances did. As a secondary contribution of this work, we show how Andorra's serology testing program was likely impacted by people leaving the country. Overall, this case study suggests how using mobile phone data to measure dynamic population changes could improve studies that rely on more commonly used mobility metrics and the overall understanding of a pandemic.

At the start of the COVID-19 pandemic, nonpharmaceutical interventions (NPIs) were 2 widely deployed in an effort to stymie the rate of new infections. These interventions 3 included stay-at-home orders and restrictions on economic activity, which were used as 4 a means to reduce contact and hence transmission rates, effectively limiting mobility. 5 Country border restrictions were also put in place to reduce the chance of importing the 6 November 5, 2021 1/22 study. 58 We develop and test multiple (SEIR) models that differ in how they parameterize 59 transmission rates based on the trips and entrances metrics developed in this work. The 60 models are simple, where their purpose is to illustrate how different types of mobility 61 information can be better incorporated into SEIR models. 62 Finally, we use the best model to simulate a hypothetical counterfactual, 63 representing a scenario where economic and border restrictions had not been put in 64 place, and trips and entrances metrics had not drastically reduced. 65 Outline. Before presenting our methods and results, we provide background 66 information, with a timeline of events around the start of COVID-19 in Andorra, and 67 the features of the country that contribute to a unique case study. We also provide 68 background information about compartmental epidemic models to guide the reader in 69 the presentation of our models. The study region of this work is the small country of Andorra, which is located in the 73 Pyrenees mountains and shares borders with only France and Spain. The country has a 74 population of approximately 77,000 [25] , yet attracts more than 8 million visitors 75 annually, mostly for tourism associated with skiing and nature-related activities [26] . In 76 addition, a large number of cross-border temporary workers reside in the country, 77 mainly employed in the tourism industry. Andorra lacks an airport or train service so 78 the primary way to enter or exit the country is by crossing the French or Spanish border 79 by car. The country is divided into 7 municipalities, called parishes. 80 Partly because of the country's small size and limited border crossings, Andorra was 81 able to implement comprehensive policies at the start of the COVID-19 pandemic, as 82 well as implement a serology testing program which more than 90% of the population 83 participated in. Furthermore, there is one telecoms provider for the entire country, 84 which contributes a comprehensive view of all mobile subscribers who spend any time in 85 Andorra, whether they are Andorran nationals or have foreign SIM cards. The telecoms 86 data and serology data are used in this work and are described in the Data sources and 87 preprocessing section. 88 Timeline of COVID-19 cases and policies 89 The first COVID-19 case in Andorra was reportedly imported via Italy and confirmed 90 March 2, 2020 [27] . Reported cases then rose rapidly in March before falling again in 91 April (see Fig 1) . On March 13, government officials ordered the closure of public 92 establishments and a quarantine was requested of the entire population. A series of 93 COVID-19 related policies followed and neighboring country borders were restricted. In 94 accordance with these policies, mobility within the country dropped and border 95 crossings ceased. Other NPIs, such as masks and hand sanitizer, were also deployed. 96 The lockdown measures in Andorra were gradually lifted in April and May, and fully 97 lifted starting June 1. Borders also reopened in June and border crossings resumed.

98 Table A .1 in S1 Appendix shows a timeline of COVID-19 related events. population of a country and one of the largest of its kind [28] . Anyone over the age of 2 103 was invited to participate in the study, including the country's temporary workers. The 104 testing was conducted in two phases: May 4 -14, and May 18 -28, 2020. The objectives 105 of the second phase were (a) to track the progression of COVID-19 between the two 106 surveys and (b) to account for indeterminate or potential false negative results from the 107 first survey. More than 90% of the population participated voluntarily in at least one of 108 the two surveys. However, an issue with the testing program was that many 109 participants in the first phase did not participate in the second, limiting the data 110 collection and impact of the two-phase study. This issue is further explored and 111 addressed in the Results section.

SEIR models and COVID-19 113 SEIR models, and their variations, are compartmental models used in epidemiology.

They have been widely used in forecasting COVID-19 transmission and modeling the 115 outcomes of government policies [14, 29, 30] . The basic concept of these models is that 116 the population is partitioned into sequential compartments, and transitions through the 117 compartments over time. This framework was first developed by Kermack compartments, as well as the transmission rate, β, and the total population size, N. The 129 standard model considers N constant, and the following conservation holds for any time, 130 t:

Transitions between compartments are modeled by a set of ordinary differential 132 equations (ODEs).

133 data. However, these mobility metrics are often based on sources that report on a small 168 fraction of the population, and where the mobility metrics are aggregated statistics 169 based on the number of trips to points of interest (POIs), which may not be the most 170 important indicators of COVID-19 transmission. This is in contrast to the telecoms data 171 used in this work, which covers all mobile subscribers within the country of Andorra, 172 and is provided as a complete and unaggregated dataset, not limited to trips to POIs. 173 We note that any of the models referenced or presented in this work are 174 oversimplifications of the complex dynamics of disease spread. They also suffer from 175 unreliable case reports data, limited by the availability of tests, and reactive to changes 176 in testing protocols [1] .

This section describes the SEIR models used in this work, and how they are trained and 179 tested. It then describes data sources and preprocessing methods.

Code and data availability 181 All aggregated metrics and code used in this work are made available and documented 182 in a public repository. The code includes analysis notebooks as well as the preprocessing 183 scripts that produced the aggregated metrics. The data reporting on individuals, which 184 was used to compute aggregate metrics, is sensitive and kept private. The aim is to evaluate the relative impact of the trips and entrances data on model 192 performance; the aim is not to build a state-of-the-art, accurate predictive model. To 193 this end, the models are highly simplified.

Comparison models 195 In SEIR models, β(t) typically represents the average number of people an infected 196 person would expose per-unit time if everyone were susceptible. In particular, β(t) is 197 used to model the transition from the Susceptible to Exposed compartments. The use of 198 β(t) in our models is captured by the following equation from Eq (2) .

We develop multiple models that only differ in how they define β(t).

In the following descriptions, b0, ..., bn are parameters of β(t) and are estimated 202 during model training for each model in which they are included.

November 5, 2021 6/22

. CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint One model uses trips without entrances data (model ii). Another model uses both 204 trips and entrances data (model iii). A model that uses neither data source is used as a 205 baseline (model i).

Each of the models use the same framework, methods, and training and testing 207 periods, described further below. In this model, the average rate at which the susceptible population is exposed can be 215 impacted by the behavior of people within the country (e.g. mobility measured in trips) 216 as well as the import of new cases (entrances). virus. The term I(t) N reflects the assumption that the likelihood of new country entrants 223 being infectious tracks with the timeline of infection rates in Andorra. This assumption 224 is based on the fact that during the study period, the timeline of infections in Andorra 225 was highly correlated with the timeline of infections in Spain and France (with Pearson 226 correlation coefficients of 0.922 (p=0.000) and 0.932 (p=0.000), respectively), and the 227 primary way to enter Andorra is through the Spanish or French borders. Furthermore, 228 telecoms data showed that 86% of entrances by foreign SIMs were either Spanish or 229 French, and when accounting for entrances by Andorran SIMs, 68% of all entrances 230 were by Spanish or French SIMs. See section A.3 in S1 Appendix.

E (t) = β(t) S(t)I(t) N (t) + S(t)f (entrances(t)) − σE(t) 218 where 219 β(t) = b0 + b1 × trips(t) b2 220 f (entrances(t)) = I(t) N × b3 × entrances(t) b4

The above functions using entrances and trips can be combined into one equivalent 232 expression representing transmissibility. We do this to simplify modeling and maintain a 233 common expression for E (t).

The SEIR framework used in this work is illustrated in Fig 2 and . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; Fig 2. Schematic representing the SEIR model framework used in this work. The population is divided into compartments where individuals transition through the compartments: Susceptible, Exposed, Infected, Removed, Case reported, where the transitions are described by ODEs (Eq (3)). reported cases at the start of the study period. Initial values for E, I, are estimated by 252 model training, along with γ and parameters of β(t). The reporting rate, r, is set to 1 11 , 253 estimated from the serology and case reports data (Data sources and preprocessing 254 section). The latent rate, σ, is set to 1 5.2 , estimated by prior work [43] . The reporting 255 delay, d, is set to 7, consistent with related works [19, 44] . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint to August, 2020. The period of March 14 -May 31 is used for model training and the 272 following 10 weeks are used for testing.

Training: Parameters and initial values for E(t), I(t) at t = 0 were fit with maximum 274 likelihood estimation (MLE). Log-likelihood was computed by comparing time series 275 values of predicted cumulative reported cases (C) to the time series of actual cumulative 276 reported cases:

Where the sum is over all days in the training data, P k,λ (k, λ) is the Poisson 278 distributed probability mass function, k is actual reported cases, λ is predicted reported 279 cases.

Parameters were optimized by minimizing the negative log-likelihood using the 281 L-BFGS-B method [45] . See section A.5 in S1 Appendix for details. The value of C(t) is corrected to the true reported cases at time t and further integration 288 over the ODEs is used to continue the simulation over the test period. The resulting C 289 estimated over the test period is compared to actual reported cases via MAPE. The program was conducted for a previous research study, in which the methods and 301 results are detailed [28] . The study was approved by the Institutional Review Board of 302 the Servei Andorra Atencio Sanitaria (register number 0720). The data were also 303 provided to researchers in our lab as part of a research partnership. The dataset 304 includes a unique identifier for each participant and results from the 1st and 2nd round 305 of tests; test results were left empty when there was a lack of participation. The dataset 306 also includes demographic information for participants, including their home parish and 307 whether they are a temporary worker. As previously described, an issue with the 308 serology testing program was that many of the participants from the first phase of 309 testing did not participate in the second phase (see table A.3) .

From the serology data, Bayes Theorem [47] was used to estimate the portion of the 311 population infected up to May. With this number and the official reported cases data, 312 we estimated a case reporting rate of 1 11 . This reporting rate is used in the epidemiology 313 models described in this work.

November 5, 2021 9/22

. CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint Telecoms data and metrics 315 Andorra has one telecoms provider (Andorra Telecom), which provided the data for this 316 study. Since they are the sole provider, the dataset covers 100% of mobile subscribers in 317 the country, unlike most telecoms datasets where the market is fragmented. Each data 318 point includes a unique ID for the subscriber, a timestamp, the coordinates of the 319 device, and nationality for the subscriber's home network. The data have been further 320 described in [48] .

The stay-point extraction algorithm of Li et al. (2008) [49] was used to reduce the 322 series of data points for each subscriber into a series of stay-points of 10 minutes or more 323 within a radius of 200m or less. The stay-points represent a more concise and reliable 324 series of places the subscriber spent time; stay-points were used to infer presence in the 325 country, dynamic population changes, and compute the trips and entrances metrics. 326 There are gaps in the available telecoms data and the resulting trips and entrances 327 metrics during the period of study (data gaps were June 28-29, and July 21-27, 2020). 328 Missing values were imputed by taking the mean across the values from the 7 days 329 surrounding each missing period of data.

Dynamic population inference and metrics: On each day, a subscriber was 331 considered present in the country if they had a stay-point in the country within a 7-day 332 window. The window accounts for unobserved subscriber devices due to a combination 333 of inactivity, lack of reception in certain areas, or noisy data. The beginnings and 334 endings of periods of presence were counted as entrances to and departures from the 335 country, respectively.

Trips metrics: Daily trips for subscribers were counted as their daily number of stay 337 points minus 1, since a new stay point is recorded when a subscriber moves beyond a 338 200m radius. Daily trips by subscribers were summed as a total daily trips metric.

Home inference: The home parish of each subscriber was inferred from the telecoms 340 data, to come up with a population count for each of the 7 parishes of Andorra. This 341 was done by first assigning each stay-point to the parish in which it was contained.

Each subscriber's home parish was then determined to be the parish in which they 343 spent the most cumulative time during night-time hours (12:00am to 6:00am). Related 344 studies of human mobility that use cellular data have employed similar methods [50] [51] [52] [53] . 345 These inferred parish-level populations were compared to the published 2020 346 population statistics [25] . There is a Pearson correlation coefficient of 0. Before presenting our main findings, we first present the start of the pandemic in 363 Andorra through a series of plots, and compare this period to the same period in 2019, 364 when Andorra experienced a normal economy with tourism. There were also already fewer total daily trips being taken at the start of March, 2020, 369 compared to 2019. This is largely due to fewer people in the country making the trips. 370 This metric also substantially dropped at the start of the lockdown. This drop was 371 partly due to even fewer people in the country making trips, and due to the government 372 imposing restrictions on movement. The number of trips gradually rose again before the 373 border restrictions were lifted in June, indicating that the population increased internal 374 mobility. The number of daily entrances to (and departures from) Andorra also The time series of reported COVID-19 cases is shown with the time series of the trips 381 and entrances metrics in Fig 1. Other studies have implied that changes in case growth 382 often lag changes in behavior and mobility metrics by 14 or more days [16, 19, 23] . while newly reported cases remained low. Case growth did not increase again until daily 385 entrances increased again when the border restrictions were lifted in June. This suggests 386 that the entrances metric is more related to case growth than the trips metric in this 387 case study. The relative predictive power of these metrics is further shown by the model 388 results (Models results section).

Serology tests and country departures 390 Andorra's nationwide serology testing program conducted in May, 2020 involved two 391 phases of testing (see the Andorra and COVID-19 section). An issue with this program 392 was that many of the participants from the first phase of testing did not participate in 393 the second phase, limiting the impact of the study. An important question for a country 394 conducting such a program might be why this happened.

This drop in participation might be particularly concerning, as we found the drop in 396 participation was more than 3 times higher among temporary workers versus the 397 general population, and results from the testing program showed that temporary 398 workers had higher seroprevalence (infection rates) versus the general population. See 399 table A.4 in S1 Appendix. This might imply that a more infected demographic group 400 was then less monitored.

By combining the serology test data with information inferred from the telecoms 402 data, we find that test participants likely left the country after their first test. 403 We counted the number of mobile subscribers, by inferred home parish, who were in 404 the country during the first and second phases of testing (May 4-14 and May 18-28, 405 2020). Subscribers were counted as present during a testing period if they had at least 406 one "stay" within the period. We estimated how many subscribers left the country after 407 November 5, 2021 11/22

. CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint A.4 in S1 Appendix.

Simple models based on the SEIR framework, were developed to compare the impact of 430 trips and entrances data on transmission rates and predicted infections.

The baseline, dummy model (i) assumes a constant transmission rate. For model (ii) 432 transmission is a function of mobility measured by trips data, and for model (iii) Models were evaluated by their prediction performance over the weeks that followed 440 the training period. This was done using MAPE, based on the framework used by Fig. 3 and Fig. 5 in [46] . Note their evaluation used 446 cumulative deaths data whereas this work uses cumulative cases data.) 447 The model (iii) using both trips and entrances data outperformed the other models 448 in all but excluding the first week that followed the training period. More importantly, 449 the model (ii) that used trips data to model transmission rates (without entrances data) 450 had results similar to, and slightly worse than, the baseline model (i) which assumed a 451 constant transmission rate. This is not surprising, as the data indicated trips were able 452 to increase without impacting transmission rates (Fig 1) . 453 This is also shown in that the best fit for model (ii) had parameters that flattened 454 the impact of the trips data, resulting in a nearly flat reproduction number, . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint (similar to model (i)). This is in contrast to the model (iii) that used both trips and 459 entrances data, and where predictions for new reported cases closely tracked with actual 460 predictions. See Fig 4. Overall, these estimated R 0 values are reasonable and within the 461 range of values estimated by previous works [57] . Table 1 . MAPE results for the 3 models. The models differ in whether they incorporate trips and entrances data to model transmissibility. Model (i) is a baseline, dummy model where transmissibility is constant, model (ii) uses trips data, and model (iii) uses trips and entrances data. All models used the same framework and methods.

As a robustness check, all models were trained and tested over an additional set of 463 training and testing periods that ended slightly earlier than those used for the main What if Andorra had not imposed a lockdown, which caused reduced mobility? What if 470 border restrictions had not been put in place, which caused a drop in entrances? 471 Overall, what if the population mobility, measured in total trips and entrances, had not 472 dropped in March?

In this section we explore such a counterfactual scenario by using the best fit model 474 (iii) from the Models results section, which uses the trips and entrances data.

The lockdown in Andorra began on March 13, 2020, and there was a large drop in 476 trips and entrances surrounding this date (see Fig 1) . We again take a simplified 477 approach to modeling, and create hypothetical trips and entrances data for a 478 counterfactual scenario where mobility and border restrictions were not put in place. 479 We do this by using the true metrics up to March 13 of 2020, and then keeping the These interventions and other NPIs showed to be effective in Andorra, as the country 493 brought case growth under control from March -May 2020, before the restrictions were 494 fully lifted. The counterfactual scenario modeled in this work shows a stark alternative 495 had the mobility changes observed during this period not occurred, with more than an 496 estimated 3x as many cases, likely overwhelming the hospital system.

Numerous other works have also used mobility data collected from mobile phones to 498 model the impacts of mobility restrictions on COVID-19 transmission. However, these 499 studies have relied on data about trips, and the data represented a small sample. Other 500 works using meta-population SEIR models where the modeled sub-populations are 501 dynamic have been based on static census data. In contrast, this work leverages data 502 collected from mobile phones that represent 100% of subscribers in a country. 503 We showed how these data could be used to build on previous works by computing 504 daily trips metrics as well as estimating a dynamic, real-time population census. We 505 then showed how these data can be used to improve upon the understanding of a 506 pandemic in two main ways.

First, these data were used in order to better understand why participation in the 508 nationwide serology testing program dropped between the first and second phases of . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint testing. The drop in participation may have been concerning as the second phase of 510 testing was intended to help better detect and track the virus. This decreased ability to 511 track the virus might have been particularly concerning because the test results showed 512 that the temporary worker population had the highest infection rates and this 513 population also had the largest drop in test participation. However, the analysis, which 514 leveraged the telecoms data to estimate dynamic population changes, suggested that the 515 decline in participation was likely due to test participants leaving the country after their 516 first test.

Second, we showed how the dynamic population data could be used to improve 518 epidemiological (SEIR) models that otherwise rely on mobility measured by trips. In 519 our contribution, we developed simple SEIR models that differed in how they used the 520 trips and entrances metrics developed through this work. These models performed well 521 compared to the 7 global COVID-19 models evaluated by Friedman et al.

(2021) [46] [14, 40, 41, [58] [59] [60] , but their purpose was not to be highly accurate; the 523 purpose of these models was to illustrate the relative importance of trips mobility data 524 versus real-time population data, namely country entrances. In particular, for the case 525 of Andorra, we find that the population was able to regain internal mobility measured 526 in daily total trips with limited growth in cases, and that total trips per day did not 527 have predictive value in the SEIR models while country entrances did. In general, the models were limited by their simplifications. For example, there was 534 likely an interaction effect between the trips and entrances metrics that was not 535 captured in the models. The models also assumed that the case identification rate (and 536 hence removal rate) and reporting rate were constant, which related works have as well 537 (e.g. [19] ). However, these rates likely changed with Andorra's increased testing. Future 538 works can more accurately model the impacts of mobility and entrances, and the 539 interaction between these metrics. This might also include incorporating data on the 540 infection rates for other countries whose populations contribute to entrances. Future 541 work can also incorporate data on testing rates to better model changes in the removal 542 and reporting rates.

Overall, this case study suggests how using mobile phone data to measure dynamic 544 population changes could improve studies that rely on more commonly used mobility 545 metrics and the overall understanding of a pandemic. . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint S1 Appendix

March 2

First COVID-19 case confirmed in Andorra.

Ski stations closed.

Partial confinement: school closures and recommended confinement.

Total confinement: all non-essential activities ordered closed. April 7

Beginning of masks delivery and progressive use by population. April 17

Allowance to 1 hour of walk in 1km radius every two days. April 20 Phase 1 reopening: low risk economic activities resume; 1000 people return to work. May 4

Start of 1st round of population serology screening. May 4 Phase 2 reopening. Additional 4760 workers return to normal activity.

Increase to 2 hours of activity every day to walk or exercise.

Start of 2nd round of population serology screening.

Phase 3 reopening: additional 3300 workers returned to normal activity. June 1

Confinement restrictions completely lifted.

French and Spanish borders opened (with restrictions on Spanish side).

Spain ended state of emergency and further lifted border controls. July 1

Remaining Spanish border restrictions lifted. September 1 Massive testing for teachers and school kids began. . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint

The parish-level populations inferred from the telecoms data for May 2020 were compared to published 2020 population statistics [25] . There is a Pearson correlation coefficient of 0.959 (p<0.001), suggesting that the telecoms data are representative of the true population. . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint

The SEIR model (iii) in this work incorporates trips and entrances data to model transmission rates. It includes an assumption that the likelihood of new country entrants being infectious tracks with the timeline of infection rates in Andorra. In order to check this assumption, we compare the rates of COVID-19 between the country of Andorra and its bordering neighbors, Spain and France, during our period of study. This is done by comparing the timeline of reported deaths per 1 million residents, where data is smoothed over a 7 day rolling window. Death reports are used instead of case reports as a more stable comparison indicator in this study and others because death reports were considered to be less impacted than case reports by the dynamically changing testing procedures which varied by country [19, 61] . The Pearson correlation coefficients between Andorra and France and Andorra and Spain are 0.916 (p < 0.001) and 0.928 (p < 0.001), respectively. We note there was an error in the Spain data with negative deaths values in late May. We changed the negative values to 0 for this comparison. Without the change, the correlation was still statistically significant with Pearson correlation coefficient 0.918 (p < 0.001). The timeline of infections is shown in Fig A. 2. During this same period, telecoms data showed that 86% of entrances by foreign SIMs were either Spanish or French, and 68% of all entrances were by Spanish or French SIMs when accounting for entrances by Andorran SIMs. The timeline of entrances by SIM nationality is shown in Fig A. 3. . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint Two cross sectional serological surveys were conducted in Andorra from May 4-28, 2020, using a rapid serological test (nCOV IgG/IgM) [28] . For each participant, the test data include the dates of their participation in surveys 1 and/or 2, the positive versus negative results for IgG/IgM antibodies for each round of testing, the participant's parish of residence, whether the participant is a temporary worker, and other demographic information. Table A .3 shows the number of participants and seroprevalence from survey 1 as well how many participants from survey 1 also participated in survey 2. Data for temporary workers is highlighted. Seroprevalence data is from previously reported results and was calculated based on the number of individuals who had a positive result of IgG and/or IgM [28] . The testing was voluntary; an issue with the testing was that many people who participated in the first survey did not participate in the second survey. Seroprevalence was higher among temporary workers. At the same time, temporary workers who participated in survey 1 were less likely to participate in survey 2 versus the general population. See We counted the number of mobile subscribers, by inferred home parish, who were in the country during the first and second survey periods (May 4-14 and May 18-28, 2020). Subscribers were counted as present during a survey period if they had at least one "stay" within the period. We estimated the rate at which subscribers left the country after the first test by counting how many subscribers were present during only the first period versus both periods.

These numbers were compared to the parish-level serology test participant populations. Namely, the portion of serology survey participants who did survey 1 but not survey 2 was compared to the estimated portion of mobile subscribers who left the country between survey periods. There is a statistically significant Pearson correlation coefficient of 0.937 (p=0.0019).

To further validate that the decline in test participation was related to people leaving the country, we repeated these tests using 2019 telecoms data. In this case, there is a Pearson correlation coefficient of 0.4928 (p=0.2612), which is not statistically significant. If the May 2020 subscribers had left the country for reasons not related to the pandemic, we would expect the correlation to be similar for the 2019 and 2020 data.

. CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint

All modeling work is publicly accessible in a Python Jupyter notebook: https://github.com/CityScope/CSL_Andorra_COVID_Public/blob/main/ analysis/SEIR_models_trips_entrances.ipynb.

Optimal parameters were estimated by minimizing the negative log-likelihood function using the L-BFGS-B method with the Python SciPy library [45, 62] . This step searches for optimal parameters by taking initial parameters which are then modified towards improved values, with specified bounds. The (minimum, maximum) bounds for γ were set to (1/10, 1/2). The bounds used for each of the parameters related to β -b0, b1, b2, b3, b4 -were (0, 2), with the observation that the parameters for the best fit models did not come up against these bounds. The (minimum, maximum) bounds for both E(0) and I(0) were set to (40, 4000) , where t=0 corresponds to March 14, the first day of the training period. The values of (40, 4000) were conservatively set where the minimum was based on reported active cases, and the maximum based on reported cumulative cases. In addition, the training routine discarded models where E(0) and I(0) values differed by more than 3000.

Before describing how initial parameters were handled, first note that the optimizing function is not convex. To avoid the optimization function terminating at local minima, a grid search was used for the initial parameters. The same grid search method was used for each of the models.

In addition to the grid search routine, another step was taken to find optimal model fits: The models using the trips and entrances data were initially fit using spline approximations of these metrics, where the splines were estimated from the true metrics using knots spaced by 7 days. Fig A.4 shows the comparison of the true metrics versus their spline approximations. This step was taken to further smooth the data and ease the computational complexity of the model fitting routine. Without this improvement, the model training was slow and rarely resulted in successful outcomes, and the estimated parameters rarely varied from their initial values. After the models were fit using the spline approximations of the data, the models were finally fit again using the true data, where the parameters found via the fitting routine with the spline approximations of the data were used as the initial parameters in the L-BFGS-B method.

. CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint A.6 Fit model values and parameters Table A . 5 shows the values of the parameters for the best fit models that resulted from model training. Fig A.5 shows the corresponding time series values representing each of the best fit models. The values include R 0 , the compartment populations, and the predicted reported cases.

The best fit models were determined as those with the best log-likelihood score when fit over the training data, where the training data period was March 14 -May 31, 2020. See section A.8 for values from the robustness check.

Models were trained based on the standard SEIR model where:

and

is cumulative case reports and accounts for reporting delay, d, and reporting rate, r.

The 3 models in this work differed based on how they defined transmission rate and the transition from S to E. See the Modeling section for details.

Model i: baseline model with constant transmission rate Model ii: daily transmission a function of daily trips data Model iii: daily transmission a function of daily trips and entrances data . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint A.7 Time from exposure to reported case The average delay in time from exposure to reported case in our model (i.e. the full transition through compartments E,I,R,C) is due to the average latent period, σ −1 , plus average infectious period, γ −1 , plus average reporting delay, d. σ −1 is set to 5.2 based on previous research [43] , γ −1 is estimated via model training, and d is set to 7 based on related work [19, 44] . Reporting delays, d, can be due to the time it takes to seek a test, for the test to be processed, and then officially reported. Note that at the start of the pandemic in Andorra, tests were sent for processing to Spain, potentially adding extra time to reporting delays.

A study in Singapore from March 2020 estimated an average reporting delay of 6.4 days (95% CI 5.8, 6.9) [44] . A reporting delay of d = 7 was implicitly assumed by Arroyo-Marioli et al. [19] . They estimated time series values for the effective reproduction number for 124 countries across the world and validated their work by correlating their estimates of R t to mobility data from the "COVID-19 Community Mobility Reports" collected by Google, where the lag between R t and mobility was 14 days (2 weeks). They assumed an SIR model rather than an SEIR model, with time from exposure to removed, γ −1 , of 7 days, implying a reporting delay of 7 days (14 − 7). We note that Arroyo-Marioli et al. produced time series estimates for R t in Andorra. However their estimates are not comparable to the R 0 estimates in this work because (a) they did not correct for the reporting error that caused an influx of 78 additional cases June 1-10 [56] , as was done in this work, and (b) their estimates were for the effective reproduction number versus the basic reproduction number. We also note that Google's "COVID-19 Community Mobility Reports" are not available for Andorra.

The best fit for the model estimated γ −1 ∼ 5.5 days (γ = 0.18). Combined with σ −1 = 5.2, d = 7 results in a total estimated average delay from exposure to case report ∼17.7 days. This is consistent with previous work over a similar study period in Andorra that studied the correlation between mobility metrics and transmission rates over various lags and found the best correlations were with mobility metrics lagged by 18 days [23] .

. CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint

As a robustness check, we trained and tested all models over an additional set of training and testing periods that ended slightly earlier than those used for the main results. The training period for the robustness check was March 14 -May 14, 2020. The models were then tested on the period that directly followed this training period. . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint

The copyright holder for this this version posted November 11, 2021. ; https://doi.org/10.1101/2021.11.06.21265955 doi: medRxiv preprint

Estimating the historical time series of infections

Ranking the effectiveness of worldwide COVID-19 government interventions

Association between COVID-19 outcomes and mask mandates, adherence, and attitudes

Inferring the effectiveness of government interventions against COVID-19

Multiscale mobility networks and the spatial spreading of infectious diseases

Forecasting the spatial transmission of influenza in the United States

Baidu Qianxi platform

The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak

Google COVID-19 Community Mobility Reports

Facebook Data For Good: Our Work on COVID-19

Shelter in Place Index: The Impact of Coronavirus on Human Movement

Citymapper Mobility Index

Mobility changes in response to COVID-19. arXiv preprint arXiv

The impact of COVID-19 and strategies for mitigation and suppression in low-and middle-income countries

Estimating effects of physical distancing on the COVID-19 pandemic using an urban mobility index. medRxiv

Public mobility data enables covid-19 forecasting and management at local and global scales

Early Detection of COVID-19 Outbreaks Using Human Mobility Data

Tracking R of COVID-19: A new real-time estimation using the Kalman filter

What about bias in the SafeGraph dataset?

Early spread of COVID-19 in Romania: imported cases from Italy and human-to-human transmission networks

Preparedness and vulnerability of African countries against importations of COVID-19: a modelling study

Mobility and COVID-19 in Andorra: Country-scale analysis of high-resolution mobility patterns and infection spread

Differential effects of intervention timing on COVID-19 spread in the United States

Estimacions de població, gener 2020 and Estadística dels censos parroquials, gener 2020. Govern d'Andorra (Reports No A001 and A003)

The World Factbook: ANDORRA; 2021

A 20-year old man is Andorra's first coronavirus case. Reuters

Mass SARS-CoV-2 serological screening, a population-based study in the Principality of Andorra. The Lancet Regional Health

How epidemiological models of COVID-19 help us estimate the true number of infections; 2020. Our World in Data

Data-driven inference of the reproduction number for COVID-19 before and after interventions for 51 European countries

Containing papers of a mathematical and physical character

Modeling infectious diseases in humans and animals

Modeling infectious diseases in humans and animals

Complexity of the basic reproduction number (R0)

Infectious diseases of humans: dynamics and control

Temporal changes in Ebola transmission in Sierra Leone and implications for control requirements: a real-time modelling study

Improved inference of time-varying reproduction numbers during infectious disease outbreaks

COVID-19 basic reproduction number and assessment of initial suppression policies in Costa Rica

Forecasting COVID-19 and analyzing the effect of government interventions

COVID-19 projections using machine learning

The TVBG-SEIR spline model for analysis of COVID-19 spread, and a Tool for prediction scenarios

Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia

Real-time monitoring the transmission potential of COVID-19 in Singapore

A limited memory algorithm for bound constrained optimization

Predictive performance of international COVID-19 mortality forecasting models

Determining the value of diagnostic and screening tests

Travel demand and traffic prediction with cell phone data: Calibration by mathematical program with equilibrium constraints

Mining user similarity based on location history

Exploring universal patterns in human home-work commuting from mobile phone data

Analyzing cell phone location data for urban travel: current methods, limitations, and opportunities

COVID-19 outbreak response: a first assessment of mobility changes in Italy following national lockdown. medRxiv

Socio-geography of human mobility: A study using longitudinal mobile phone data

An interactive web-based dashboard to track COVID-19 in real time. The Lancet infectious diseases

Coronavirus Pandemic (COVID-19); 2021. Our World in Data

Els tests d'anticossos permeten diagnosticar 78 positius de la COVID-19, que podrien haver contagiat unes 360 persones; 2020. Govern d'Andorra

COVID-19 R0: Magic number or conundrum?

Los Alamos national Laboratory COVID-19 Team. Los Alamos Natinoal Laboratory COVID-19 Confirmed and Forecasted Case Data

MRC Centre for Global Infectious Disease Analysis. Imperial College COVID-19 LMIC Reports

Fast and Accurate Forecasting of COVID-19 Deaths Using the SIkJa Model