The ongoing SARS-CoV-2 pandemic, first detected at the end of December 2019 in Wuhan, China, has caused substantial morbidity and mortality globally. In the five months since the beginning of the outbreak, 188 countries/regions had reported the COVID-19 disease and there have been around 350,000 confirmed deaths and over 5.6 million confirmed cases globally [1, 2] .
In the absence of a vaccine and effective medical treatment, unprecedented national and community nonpharmaceutical, social distancing interventions (NPIs), including travel restrictions and border closures, limitations on gatherings and public events, and school, workplace and business closures, were introduced across the globe in addition to recommendations about personal hygiene and behavior (washing hands, covering cough, keeping 6-feet distance, wearing masks, self-isolation) [3] . The goal of NPIs is to reduce the transmission of the SARS-CoV-2 virus and thus slow the exponential growth of the epidemic [4] . Delaying the onset and reducing the magnitude of the infection peak helps to better utilize medical resources, to avoid the overburden of the health care system and to effectively "buy time" until a cure and/or vaccine is available. While reducing the impact of COVID-19 on population morbidity and mortality is the desired health outcome, the economic and social consequences of social distancing interventions are devastating. Therefore, it is critical to understand if social distancing interventions and policies introduced by governments and public health authorities are effective in achieving desirable health outcomes. Even more, the ideal scenario would be to establish a distinct quantitative relationship between a readily available metric of the population's response to policies in place and the growth rate of the epidemic.
One such metric is population mobility, which can be assessed and quantified either on an individual or a population level and compared before and after the implementation of various NPIs. The first studies about the impact of social distancing and mobility on COVID-19 epidemic came from China. A study by Kraemer et al. [5] used mobility data and fine-grained epidemiological data at the level of individual people to explain the spread and growth rate in early stages of the COVID-19 epidemic in China. The study found that variation in the COVID-19 growth rates outside of Wuhan could be almost entirely explained by population movement from Wuhan. The study also explained how drastic social distancing and the control of movement measures across the country resulted in negative growth rates. A related work that used fine-grained mobile phone data showed that movement of undocumented cases contributed to over 80% infections during the early stages of the epidemic in China [6] . Despite the useful insights, it is difficult to replicate these studies in other countries due to obstacles in obtaining fine-grained, individual-level linked mobility and case data.
Several publications report on the relationship between implementation of social distancing measures and population mobility. Klein et al. [7] used Cuebiq company aggregated and de-identified data collected from 17 million mobile devices between January, 1 and March 25, 2020, to conclude that within a week of implementation of social distancing, the major U.S. metropolitan areas experienced on average 50% reduction in typical commutes to/from work. A similar study in Italy [8] showed that during the early phase of the outbreak when only mild measures were in place (e.g., school closures) the traffic to/from the most affected provinces declined by about 50%. Recently, the public release of population mobility data by Google (as Google Mobility Index, GMI) [9] and Apple (as Mobility Trends Report) made it possible to understand how population mobility changed across the globe in response to social distancing interventions. While it is clear that most countries accomplished a significant reduction in population mobility, elucidating the precise relationship between population mobility and COVID-19 growth is still an open question.
There are in general two types of studies attempting to model and predict behavior of COVID-19 epidemic. The first type of studies relies on mathematical models that realistically reproduce disease dynamics. For example, Chinazzi et al. [10] used a stochastic disease transmission model to show that strict travel restrictions delay the disease spread by only a few days. Linka et al. [11] used a network simulation model to understand the impact of travel restrictions on COVID-19 spread in Europe; and Tuite et al. [12] used a mathematical model to study the epidemic spread in Canada. The mathematical models rely on first principles that specify fundamental dynamic of an epidemic typically require multiple parameters, some of which are selected based on the known information about the COVID-19 and others are fitted using historical data. The resulting dynamics for all but the simplest models cannot be described in a closed form, which reduces their interpretability. Despite calibrating their parameters to data and constraining their values based on domain knowledge, they can still overfit and produce incorrect behaviors if the number of parameters is large compared to the size of available data. Another type of studies [13, 14] develops statistical prediction models that are simple to describe but are not rooted in first principles. A representative of such an approach is a widely publicized Institute for Health Metrics and Evaluation (IHME) model [15] which assumes that the logarithm of positive cases follows a quadratic function with the global maximum and thus optimistically predicts a quick decline of an epidemic. The UTexas model [16] extends the IHME model by allowing the quadratic function to be governed by the population mobility. It is notable as one of the few studies that directly link population mobility and number of cases. A related approach [17] provides a more elaborate model with more realistic assumptions but a larger number of parameters that are fitted using Bayesian approaches. A common weakness of statistical models is that they are not derived from first principles such that, while they can result in accurate short-term forecasting, they underperform in reproducing the natural behavior of an epidemic and are unreliable for longer-term forecasting. Our objective is to combine a mathematical model of epidemic dynamics and a statistical prediction model, to elucidate a quantitative relationship between COVID-19 infection growth rate and mobility as a proxy of social distancing.
As illustrated in Figure 1 our modeling framework is based on 3 interconnected models -SIR model, Social Distancing Model and Statistical Model (see Materials and Methods for details). The augmented SIR (Susceptible-Infected-Recovered) epidemiological model establishes the underlying dynamics of an epidemic where individuals can be classified by their epidemiological status as susceptible to the infection (S), infected and therefore infectious (I), and recovered and no longer infectious (R). The central quantity is the infection growth rate ! ( ), which is the relative rate of increase in total number of infected individuals ( ) at time that can be calculated as ! ( ) = ( ( + 1) − ( ))/ ( ) . The infection growth rate depends on the transmission rate ( ) (how quickly an infected person infects somebody else), the recovery rate (how quickly an infected person recovers), and the susceptible ratio (what fraction of population is susceptible to becoming infected). In our framework, we assume we are observing an initial stage of a new epidemic, so the susceptible ratio can be assumed to be close to 1. Then, the infection growth rate is approximated as ! ( ) = ( ) − (see Materials and Methods). Under some additional assumptions, the growth rate " ( ) in newly infected individuals ( ) equals the growth rate ! ( ) in total infections, ! ( ) = " ( ). By augmenting the SIR model, it can be shown that the growth rate # ( + ∆) in new fatalities ( + ∆) at time + ∆ can be approximated by the growth rate " ( ) in new infections at time , # ( + ∆) = " ( ), where ∆ is the average delay (measured in number of days) between the time of infection and the time of death, if it occurs. Thus, it is also approximately correct that # ( + ∆) = ! ( ).
The social distancing model establishes that the transmission rate ( ) depends on the behavior ( ) of the impacted population, which is not easily observable. Population mobility ( ) is an observable, standardized, and publicly available quantity that depends on the population behavior ( ). The value ( ) = 1 corresponds to the baseline mobility immediately prior to the COVID-19 outbreak in that country, ( ) = 0 corresponds to the theoretical scenario when the population stops moving, and any number between those two extremes is a fraction of the baseline mobility. Because population mobility information is publicly available for multiple countries and regions (e.g. with Google Mobility Index) we use it as proxy for the behavior ( ). We hypothesize that ( ) can be used as a predictor of ( ).
In Figure 1 , the quantities in the boxes surrounded by single-line rectangles are either not readily observable (such as ( ), ( ), ) or are unreliable (such as ! ( ), " ( )) and can be treated as hidden variables. The quantities in the boxes surrounded by double-line rectangles ( ( ) and # ( )) are readily available and sufficiently reliable and can be treated as observed variables. The full-line arrows are causal relationships. The single-dashed line is an indirect causal relationship. By combining the social distancing model and the SIR model, the statistical model establishes a link (double-dashed line arrow) between the two observable variables: mobility ( ) at time and death growth rate # ( + ∆) at time + ∆ as # ( + ∆) = ( ( )), where the function is fitted from the observed data about mobility and daily death counts using a regression algorithm. Since the augmented SIR model established that # ( + ∆) = ! ( ), ( ( )) can be used to predict the current infection growth rate ! ( ) (orange dotted arrow). In addition, we can also directly estimate and ( ) as the intercept and the remaining terms of ( ( )), respectively.
Our framework assumes that during the period covered by our data (from mid-February to May 9 th , 2020) the 14 countries we analyzed were in the initial stages of COVID-19 epidemic, with large majority of population still being susceptible. Based on several population-based surveys that test for the presence of SARS-COV2 specific antibodies, the crude estimates of the proportion of individuals in the general population who developed antibodies and therefore have been infected is below 10% [18, 19] . Important exceptions are COVID-19 geographical hotspots like New York City where some reported prevalence of up to 24.7% or in very selected, small clusters of high-exposure risk populations like health-care workers and teachers, where estimates indicate that up to 50% or more individuals may have already been infected in some localities [18, 19] .
The relationship # ( ) = " ( ) = ! ( ) is valid under the assumption that population mobility ( ) changes slowly on a time scale of days. If ( ) changes rapidly, as was the case during the first days of social distancing in most countries, it is difficult to write the relationship in a closed form. Due to this observation, before fitting the regression function we remove data points corresponding to days when population mobility changes rapidly.
The proposed framework relies on data about COVID-19-related deaths to infer the relationship between population mobility and COVID-19 growth rate. We note that using data about infection cases rather than deaths would result in a simpler estimation problem because it would obviate the need to use a delay ∆ in the modeling. However, for several reasons, the publicly available numbers of confirmed COVID-19 cases provide a biased and unreliable view of the epidemic spread, due to the inability to accurately diagnose/identify all patients with the disease [20] . First, the majority of COVID-19 cases are only mildly symptomatic or even asymptomatic so they do not necessarily seek medical care; yet, these undocumented cases substantially contribute to the spread of the infection [6] . Second, laboratory testing for confirmation of the presence of the virus is the accepted method for verifying disease presence; however, it suffers from numerous sources of variability that are typically not reported in public data and are hard to control for. For example, testing capacity is a rate-limiting step for identifying new cases and it widely differs between countries. Similarly, different clinical criteria exist across countries that trigger lab testing so accessibility to testing varies across and within countries and across time. Further, numerous testing methods are used globally, all with different or unknown analytical accuracies for SARS-COV2 detection. Different test processing times (from swab to result) and delays in reporting of testing results (based on the day the test was administered or the day when the test results were received) add to the inconsistency and variability in reported number of daily new cases. Taken together, the inability to consistently determine the number of COVID-19 cases could not only underestimate true disease prevalence but also makes it difficult to compare across different countries and stages of the epidemic.
In contrast, although not perfect, the number of COVID-19-related deaths and death growth rate derived from it present a more reliable statistic for tracking the impact of COVID-19 across regions and times. Death reporting is generally standardized, and countries follow the "cause of death" classifications from the WHO's International Classification of Diseases guidelines [2, 21] . However, countries also have their own guidelines on how and when COVID-19 deaths should be recorded and the definitions of what is considered a "COVID-19 death" may have changed during the outbreak within a particular country [22, 23] . For example, as of April 14, 2020 US CDC expanded the definition of COVID-19 deaths to include both confirmed cases and "probable cases" [23] , adding clinical and epidemiologic evidence as being independently sufficient instead of COVID-19 laboratory confirmation only. Several countries, like UK, Italy, Belgium and France [22, 24] were initially recording only in-hospital deaths but as the epidemic progressed they included out-of hospital deaths (i.e., nursing homes), too. Similar to case reporting, there are also delays and seasonality effects in death reporting, especially over the weekends and public holidays, where underreporting occurs during weekends and holidays followed by overreporting in the days after [25] . Additionally, for some countries the daily death count represents recorded death count on a particular day, not the actual number of deaths that occurred on that day. Our data preprocessing procedure is (see Materials and Methods) specifically developed to alleviate some of the issues related to death reporting.
We downloaded daily counts of COVID-19 deaths for 14 countries collected up to May 9, 2020 (see Table S1 for a list of country names and important events and statistics) from Our World in Data [2] . We downloaded population mobility data at the country level from Google's COVID-19 Community Mobility Reports (https://www.google.com/covid19/mobility/) [26] , that are referred to as the Google Mobility Indices (GMIs). The data was released on May 9, 2020 and included mobility data from February 15, 2020 to May 3, 2020. For this study, we selected 14 countries if both population mobility and death counts data were available, and the death counts exceeded 100 for at least 3 days up to May 9, 2020.
As seen in Figure S1 .a, there are 6 categories of GMI, measuring population mobility with respect to retail and recreational venues ("retail"), groceries and pharmacies ("grocery"), parks ("parks"), transit centers ("transit"), workplaces ("work"), and places of residence ("home"). The changes in population mobility after social distancing differ across the 6 GMI categories, with retail, transit, and work categories being impacted the most overall. Because the GMI time series show strong weekly seasonality, we reduced the seasonality by smoothing over a 7-day window (see Materials and Methods), as seen in Figure S1 .b. Figure S1 and [27, 28] , the smoothed GMI transit mobility decreased below 25 on March 15 and stayed below 25 (we refer to it as Phase 2 later in the paper). Figure S2 shows COVID-19 death counts and death growth rates (in log-scale) for each country. As seen from the upper panels, the originally reported daily death counts are very noisy with outliers and weekly seasonality. The strongest weekly seasonality in the data is observed for Germany, Great Britain, Netherlands, and Sweden, which follow a pattern of evident dips in reported deaths during weekends, followed by subsequent jumps in reported deaths that compensate for the dips. In addition to seasonality, several countries such as France and United States have days with significant outliers. For example, during the period of April 3-6, 2020, the number of reported deaths in France is the following sequence: 471, 2,004, 1,053, and 518, indicating reporting inconsistencies. In addition, Sweden has several distinct underreporting periods, the largest surrounding the Easter holiday (April 12-14, 2020), that are not followed by compensatory upticks in counts. Because of the lack of consistent information about reporting issues in 14 countries, we developed a country-agnostic data preprocessing approach (see Materials and Methods). The preprocessed daily counts ( Figure S2 , middle panels) capture the overall trend and enable us to obtain a more reliable estimate of death growth rates ( Figure  S2 , bottom panels). It is evident from Figure S2 and Table S1 that all 14 countries experienced a sizeable decrease in death growth rates during the observed period, from the highs in range from 0.144 (Turkey) to 0.259 (Italy) at the start of the outbreak to the lows in range from -0.043 (Spain) to 0.058 (Brazil) by second half of April. For our analysis, we excluded all days when the corrected plus smoothed death counts were below 10 because they cannot reliably estimate death growth rate, and the last 3 days due to the effects of smoothing.
The objective was to use the collected data to fit a regression function that predicts the death growth rate from GMI mobility in all 14 countries. Using the preprocessed GMI and death count time series, we created a separate data set for each of the 14 individual countries. We created a different data set for each choice of daydelay ∆ ∈ [1,30] . Given ∆ , we formed data points ; ( ), # ( + ∆)< for each country. Mobility was calculated as ( ) = ( )/100, where ( ) is one of the 6 GMI categories. For each data point, we also recorded the variance of # ( ), needed to provide uncertainty weights for the regression, and daily change in population mobility, | ( + 1) − ( )|, needed to decide which data points to remove (see Materials and Methods).
In the first experiment, we assumed that the same delay ∆ is appropriate for all countries. Thus, for each choice of delay ∆ we joined data points from all 14 countries. Depending on ∆, the total number of created data points ranged from 622 to 705. To fit a regression function that predicts # ( + ∆) from ( ), we had to decide what is the appropriate delay ∆, which of the 6 GMI categories to use for population mobility ( ), which order k of polynomial is appropriate for the regression function , and what threshold to use to remove data points with | ( + 1) − ( )| > . Using a grid search, we determined that the best choices are ∆ = 20 days, and its accuracy is + = 0.784 (the adjusted coefficient of determination, which measures the fraction of the variance in # explained by ; + = 1 is the theoretical maximum). We refer to this model as CovTM2 (Covid-19 growth rate model with Transit GMI with order-2 polynomial without the linear term). In Figure 2 we show the CovTM2 model with superimposed additional country-specific information to help interpret the results. CovTM2 accurately predicts # ( + 20) from ( ) for 9 of the 14 countries, namely Italy, France, Spain, UK, Belgium, Turkey, Netherlands, Germany and USA. Moreover, the results for Italy reveal that CovTM2 accurately predicts both phases of social distancing in Italy. Of the remaining 5 countries, Sweden's growth rate is lower than predicted from transit GMI, while Brazil, Mexico, India, and Canada have growth rates higher than predicted by the model. A visual inspection of Figure 2 and the high accuracy of + = 0.784 indicate that GMI transit is a good proxy for measuring the effect of social distancing on the growth rate. It also shows that each country may have additional intricacies influencing the COVID-19 growth that are not fully captured by the CovTM2 mobility model. is ,-= 23.8 days; b) that the transmission rate at the baseline mobility level is * = + = 0.214 1/day (95% CI [0.186, 0.244] 1/day), indicating that the expected time needed for an infected person to infect another person is * ,-= 4.67 days, c) that the baseline reproduction number is * = * / = 5.06 (95% CI [4.38,7.19] ), which is the expected number of people infected by a single infected person, and d) that the critical mobility equals . = P− * / + = 0.444 (95% CI [0.361, 0.476]). The confidence intervals are bootstrap estimate (see Materials and Methods) based on 1,000 replicates (see Figure S3 ). The optimal delay ∆ = 20 days is relatively close to the recovery time ,-= 23.8 days, even though the regression model did not impose it. All those estimates are comparable to what is currently empirically known about the transmission dynamics and progression of the COVID-19 disease [29] [30] [31] [32] . These results indicate that population mobility needs to be reduced to less than half of its baseline value to be able to start reducing the growth of COVID-19 epidemic.
To demonstrate that the choice of parameters for CovTM2 (∆ , GMI category, polynomial order, and ) is stable we show the sensitivity analysis results. Figure S4 shows the + accuracy when the delay ∆ ranges from 1 to 30 days for each of 6 different categories of GMI, while keeping = 2 and = 0.05. The accuracy for all 6 categories changes smoothly over the range of delays. The work category of GMI achieves a comparable accuracy to CovTM2 of + = 0.773 for ∆ = 18 days, while the accuracy of other categories is substantially lower, with the park category being the least predictive (peaking at + = 0.326 when ∆ = 25 days). Figure S5 shows changes in accuracy as a function of the threshold for data removal. It could be observed that has a minor influence on the regression model. This is a strong indicator that the rate of change in population mobility does not have a pronounced negative effect on modeling, although our analysis (see Materials and Methods) implied it could be an important factor. We decided to use a removal threshold = 0.05 for CovTM2 because it results in a slightly higher accuracy than other choices. Table S2 shows the fitted regression models for 5 different polynomials and their accuracies when other parameters are the same as CovTM2. Figure S6 depicts the 5 learned regression models on a scatterplot. As can be seen, the first order polynomial is inferior, while 2 nd and 3 rd order polynomials have the same accuracy. The 2 nd order polynomial has a negative linear term and is therefore not an increasing function of ( ) in the whole range from 0 to 1. Thus, we also fitted the second order polynomial without the linear term, which is the best-fit quadratic polynomial that is constrained to be an increasing function of ( ) in the whole range. This constrained second-order polynomial is used in CovTM2.
In the second set of experiments, we explored a country-specific relationship between mobility and growth rate. In order to reduce the degrees of freedom in polynomial fitting, we set the constraint that each countryspecific order-2 regression function has an identical intercept * , by assuming that the recovery rate = − * is a general property of the epidemic and should not be country-dependent. We also noticed that the average post-social distancing growth rate in France and Spain is slightly lower than the intercept of CovTM2. Thus, we set * for all country-specific models to −0.05, slightly lower than that of CovTM2 (but, still within its 95% CI). In addition, we also allowed for the delay ∆ to be country-dependent, as long as it remained in the range from 15 to 25 days, to remain close to the overall optimum ∆ = 20 days. We selected the ∆ for each country that results in the highest predictive accuracy.
We show the results for Italy in Figure 3 and for all 14 countries in Figure S7 and Table S3 . As expected, for the 9 countries that were already a good fit for the CovTM2 model, their country-specific model is very similar to CovTM2. For Italy, in particular, we can observe that after adjusting the delay to ∆ = 15 days the fitting becomes highly accurate ( + =0.984) and closely follows both phases of population mobility reduction. For Brazil, Canada, and Mexico, their models are almost linear, while for India the resulting model has a concave shape. The reported country-specific + accuracies are very high for all 14 countries, ranging from 0.755 to 0.984 with the median of 0.903. In addition to the GMI transit category, we also trained country-specific models with other GMI categories. In Table S3 we list the highest obtained + with the corresponding GMI category and delay. Interestingly, for 11 out of 14 countries, GMI transit is the best or very close to the best category. The only exceptions are Belgium (GMI retail increases + by 0.064), Mexico (GMI work increases + by 0.074) and Turkey (GMI parks increases + by 0.049). By using the country-specific delays from Table S3 , we created a data set for each country using its best delay, and then joined those data sets to create a delay-optimized data set with 642 data points, using GMI transit mobility as ( . We refer to this model as CovTM2d (with country-specific delays) and this model is illustrated in Figure S8 .
Despite the widespread implementation of social distancing measures across the globe to contain the COVID-19 pandemic, its quantitative relationship with the epidemic growth rate has remained an open question. The results of this study indicate that there is a quantitative relationship between the population mobility, measured by the Google Mobility Index (GMI), and the growth rate of COVID-19 infections, namely: the growth rate is a quadratic function of the mobility index. This function was obtained via a regression model trained on publicly available data from 14 countries, and it provides a direct estimate of the important parameters of the epidemic such as the mobility-dependent transmission and mobility-independent recovery rates, as well as the critical mobility. An important finding of this study is that a single function that relates mobility and infection growth fitted on 14 countries achieves high accuracy. Given the significant differences in combination of social distancing measures introduced across these countries, as well as underlying social, demographic and cultural differences, the quality of model fit is very remarkable.
The regression model was derived from first principles encoded in an SIR compartmental epidemiological model, which is one of the simplest models that can still generate a dynamic that qualitatively resembles an epidemic such as COVID-19. In fact, thanks to the SIR simplicity and augmenting it by the mild assumptions needed to model the death growth rate, the resulting regression model is easy to interpret, easy to fit, and accurate, thus hitting a sweet spot between modeling simplicity and accuracy. There is a range of extensions of the basic SIR model that can generate more realistic behaviors, such introducing more compartments to model the effects of incubation or stages of an infection, or using heterogeneously mixing populations to account for age-dependent and spatial effects. However, more complex epidemiological models have a larger number of parameters that would necessitate more complex data and result in a less interpretable statistical model. Our study relied exclusively on population mobility as a proxy for changed behavior due to social distancing. Using additional population behavior (such as data about mask usage) and environmental (such as weather) variables would increase the expressive power of our statistical model. However, we are currently not aware of other high quality behavior variables with a global coverage that could have been used in our study. Due to the small size of our data, we reason that environmental time series (although they are publicly available and of good quality) could lead to overfitting of our model.
Our results imply that the COVID-19 transmission rate is a quadratic function of the GMI transit mobility. While understanding why is a topic for future research, a quadratic function would arise if we assumed that the population is homogeneous (as is the case in the SIR model) and that the sole impact of social distancing is in reduction of number of contacts between infected and susceptible people (infection probability per each contact does not change). In this case, the transmission rate would scale with the number of contacts, which would scale as the square of population mobility if the fraction of infected people is small.
The quadratic relationship implies that while an initial moderate mobility reduction would lead to a sizeable reduction in COVID-19 infection growth rate, any additional mobility reduction would have diminishing returns. Another undesirable implication of the quadratic relationship between population mobility and infection growth rate is that, starting from a given level of reduced population mobility, any increase in mobility will have a much larger negative effect (rapid increase in growth rate) than any decrease in mobility (mild reduction in growth rate). Reversing the growth of an epidemic requires a significant reduction in mobility (critical mobility below 50% of the baseline pre-COVID-19 mobility). To achieve a robust decline in infections and deaths, population mobility would need to go much lower and become comparable to the reductions achieved through strong social distancing measures in Italy, Spain, and France.
Overall, our results imply that any attempt to eradicate the epidemic before much of the population gets infected and in the absence of pharmaceutical interventions, would require disciplined, coordinated, and prolonged social distancing. Another option would be to allow a sizeable fraction of the population to get infected, in what case the susceptible ratio = / would be significantly below 1. Because in the SIR model the transmission rate is multiplied by the susceptible ratio when calculating the infection growth rate, this would lead to a reduction in the infection growth rate even if the population mobility is kept constant.
It is becoming evident that population mobility reduction due to social distancing is an overwhelmingly costly intervention. Thus, it is of utmost importance to implement less crude NPIs or to combine social distancing with other forms of NPIs. Our results on individual countries indicate that population mobility reduction is not the only effective control mechanism. First, it is evident that baseline COVID-19 growth rate prior to social distancing differed across countries (before social distancing the daily growth rate differed from 14.4% to 25.9%). Moreover, growth rates observed after introduction of social distancing show that mobility reduction did not have the same effect in all countries. Five countries, in particular, were outliers with respect to the learned regression function (Sweden, Canada, Brazil, Mexico, and India) and they might offer useful insights about COVID-19 control. Probably the most important research objective outside of finding a cure or vaccine will be related to precision social distancing: to understand how governments and public health authorities could control or modify select aspects of society and how individual people should fine-tune their behavior in a way that minimally impacts the economy and lifestyles while at the same time curbing or reversing the growth of the epidemic.
For a time-dependent function ( ), its rate of change is defined as time derivative / ( ) = ( )/ and its (relative) growth rate is defined as 0 ( ) = ′( )/ ( ). In this study we assume that the unit of time is one day. If the function changes slowly within a unit of time (a day), as is the case for the functions relevant for the COVID-19 epidemic, we can approximate the rate of change via finite differences of daily changes as ( + 1) − ( ), and the growth rate via ; ( + 1) − ( )</ ( ) or via log ( ( + 1)/ ( )).
Our framework for modeling the impact of social distancing on the COVID-19 epidemic growth, outlined in Figure 1 , depends on three models described below (epidemiological model, social distancing model, and statistical model).
An established way to capture the basic dynamics of infectious diseases such as COVID-19 is via compartmental epidemiological models. We start from the basic Susceptible-Infectious-Recovered (SIR) model,
which describes the temporal evolution of the numbers of susceptible, S(t), infected, I(t), and recovered, R(t), individuals at time t from a population of size N, such that S(t) + I(t) + R(t) = N. The parameter is the transmission rate and is the recovery rate. We assume that = ( ) is time-dependent due to its dependence on social distancing (see below), while γ remains constant (in the absence of a cure).
In this work, we are studying COVID-19 during the initial stages of epidemic, when a large majority of people are still susceptible. Thus, we make a simplifying assumption that ( )/ is sufficiently close to 1 and that
where the infection growth rate equals ! ( ) = ′( )/ ( ) = ( ) − . From this equation, we obtain the basic reproduction number * = / , which varies in time as changes. From (2) we also obtain the number of new infections per day happening at time , as ( ) = ( ) ( ).
A key shortcoming of ordinary differential equation models like (1) is that they do not explicitly capture delay and temporal non-locality effects. While this can be justified for the infection dynamics, for COVID-19 the time between infection and death is usually weeks, and thus the simple R(t) dynamics in the SIR model do not suffice. We therefore retain SIR as a simplified model of COVID-19 infections, but use a more realistic model for the dynamics of COVID-19 deaths. Let d(t) denote the number of new deaths per day happening at time t. We assume that the probability of death for a newly infected individual is # and that the time of death is distributed over some future time, which gives rise to the following law
Here the shift parameter ∆ is the average number of days between infection and death (if death occurs) and is a kernel (normalized, ∫ ( ) = 1, and centered, ∫ ( ) = 0) describing the distribution of times of deaths around the delay ∆.
Based on the infection law (2) and death law (3), we obtain for the death growth rate 
In the special case when ( ) is constant in time around time , i.e., ′( ) = 0, this relationship substantially simplifies to # ( + ∆) = ! ( ), which provides a direct expression for death growth rate as a time-delayed function of infection growth rate ! ( ). In addition, under the same assumption that ′( ) = 0, it also follows that " ( ) = ! ( ), such that we can write ! ( ) = " ( ) = # ( + ∆).
It must be stressed that this relationship (5) does not require any knowledge of the virus's lethality # .
In turn, when and are varying in time (i.e., immediately after the implementation of social distancing measures), there is no simple relationship between death growth rate and infection growth rate (the additional integral terms in (4) do not vanish). Thus, the data from periods during which experiences rapid change will be ignored in our statistical model, as will be described later in this section.
The transmission rate ( ) in the reduced SIR model (2) depends on the population behavior ( ) at time . However, it is not clear at the moment what population behavior variables would be sufficient to precisely predict ( ) and how those variables could be reliably measured over many countries. Instead of using population behavior, in this study we rely on the Google Mobility Index (GMI) [9] , which is a publicly available data set about population mobility in a large number of countries calculated using a standardized methodology. We hypothesize that GMI can be used as a proxy for ( ), because population mobility reflects some aspects of population behavior change in response to social distancing. Specifically, to obtain ( ) we use the publicly available Google Mobility Index (GMI) [9] . The ( ) variable is created by adding 100 to the original GMI values and can be interpreted as the percent of the baseline mobility. Thus, ( ) = 100 represents the baseline mobility and ( ) = 0 refers to the hypothetical case when the population stops moving completely. We define the mobility variable as ( ) = ( )/100, such that it measures the fraction of the baseline mobility.
Our main hypothesis is that the transmission rate can be statistically described as a function of population mobility, i.e. ( ) = ; ( )<. We denote * = (1) as the baseline transmission rate. We constrain to be a non-decreasing function of for 0 < < 1. Moreover, we require that (0) = 0, reflecting the assumption that the infection transmission halts in the theoretical limit when the population stops moving. Thus, we model ( ) as a polynomial with zero intercept. Hence, we obtain the infection growth rate ! ( ) = ; ( )<, where ( ) = ( ) − .
For the data fitting, we choose polynomials of order , i.e., ! ( ) = 2 ; ( )<, where is a hyperparameter that is determined experimentally. Moreover, we constrain the function 2 ( ) to be non-decreasing on the interval 0 < < 1, as there is no principled reason why an increase in mobility should ever result in a decrease in transmission rate.
By combining the social distancing model and the SIR epidemiological model, the resulting relationship between population mobility and death growth rate is
It states that death growth rate at time + ∆ is a function of mobility at time . We note that both ( ) and # ( ) are measurable and publicly available, unlike other quantities such as ( ), ( ), , that are not easy to measure, or ( ), ( ), and ( ), that are unreliable. Thus, it will be possible to fit the function 2 using observed data about population mobility and death growth rate.
Once a functional relationship between death growth rate and mobility is learned from data, from (5) we also know the infection growth rate as ! ( ) = 2 ( ( )). Moreover, the recovery rate and the transmission rate ( ) can also be inferred from 2 . The root of the function 2 ( ) determines the critical mobility . (where 2 ( . ) = 0) below which the infection decays.
The Google Mobility Index (GMI) exhibits a strong weekly seasonality, where mobility during weekdays is noticeably different from weekends. Because the kernel G spreads deaths around ∆ days after infection, weekly seasonality in mobility becomes a source of covariate noise in our statistical model. To reduce this noise we remove seasonality by smoothing GMI with moving averaging over 7 days, (344$5 ( ) = ∑ ( + 6 "7,6
). In our experiments, we use ( ) = (344$5 ( )/100. We note that it has been estimated that COVID-19 deaths occur in a range from 2 to 8 weeks after infection, but the actual time distribution is not known. Thus, the 7-day smoothing of mobility should not have an adverse effect on the growth rate prediction.
Daily death counts ( ) are very noisy (see Figure S2 , top subplot) and require data cleaning and processing. Our approach has three steps: (Step 1) outlier removal, (Step 2) undercount/overcount correction, (Step 3) smoothing. For outlier correction (Step 1), we manually corrected only 10 entries. We replaced a single zero entry occurring in Germany, Spain, Sweden, and U.S. by an entry from the previous day. We clipped two U.S. entries over 3,000 to 3,000 and four France entries over 1,000 to 1,000.
For undercount/overcount correction (Step 2), we allowed transferring some overcount deaths to previous days if they had an undercount. Because deaths are sampled with probability from the newly infected, we treat ( ) as the binomial random variable. When ( ) is sufficiently large, it could be approximated as a normal random variable with mean ( ) and variance ( ). By recognizing that the growth in deaths is exponential, we first calculate the local geometric mean for deaths, _ ( ), as log`_( )a = . The same correction procedure is repeated for every overcount day. Then, the entire process is repeated a several times until convergence (no more reassignments of deaths). After the convergence, we compare the difference ( ) − _ ( ) for each day, and, if it is beyond one standard deviation, P _ ( ), we set ( ) at one standard deviation from _ ( ).
After the correction (Step 2), the time series ( ) will still have seasonality effects. Thus, in Step 3, we calculate its local geometric mean _ ( ) over a moving window of size 7 (see Figure S2 , middle subplot) and
use it to calculate the death growth rate as # ( ) = log . Thus, # ( ) is a log-ratio of two binomial variables, which can be approximated [33] as a normal random variable with mean # ( ) and variance a, where can be neglected because it is small. Eventually, we removed # ( ) on days when _ ( ) < 10 and the last 3 days (because of the smoothing).
Given a time series of preprocessed daily mobility indices, { ( )}, death growth rates, { # ( )}, and assuming a second-order polynomial, we can fit the following regression function,
The parameters * , -, and + can be fitted using least squares and ∆ is a hyperparameter that can be obtained by line search. We impose a constraint that 2 is an increasing function of in range 0 < < 1. Assuming ∆ is known, the data set for fitting consists of matched pairs ( ( ), # ( + ∆)). Because equation (5) is correct only when ( ) does not vary too rapidly, we remove data points when | ( + 1) − ( )| > , where is the removal hyperparameter. Because # ( ) is treated as a random variable with variance that depends on ( ), to fit (7) we use Weighted Least Squares (WLS), where the weights in WLS are the reciprocals of the variance. To reduce overfitting we use feature selection to select the most appropriate GMI category as ( ) and discard the others.
Because our data set consists of time series from 14 countries, we use bootstrapping with two levels of sampling to provide confidence intervals for our regression model. For each bootstrap sample, we first select 14 countries randomly with resampling from the 14 countries. Then, for each sampled country we select data points randomly with resampling. We fit the regression models on 1,000 bootstrap samples and report 95% confidence intervals for regression parameters and critical mobility. Table S1 . COVID-19 death and mobility statistics for 14 selected countries: Countries and their abbreviations, maximum recorded number of daily deaths (Count D Max), date of the maximum recorded number of deaths (Date D Max), date when number of daily deaths reached 10 (Date D ≥ 10), date when population mobility measured by GMI transit value decreased below 85 (Date GMI < 85%), minimum smoothed GMI transit mobility (GMI Min), average estimated death growth rate before social distancing corresponding to the days when GMI transit was over 90 ( ! ( ) Before), and average death growth rate after social distancing corresponding to the days when GMI transit was within 10 of the GMI Min ( ! ( ) After). Smoothed Deaths Figure S3 . Plot of regression models fitted on 1,000 bootstrap replicates of the original data corresponding to the CovTM2 model. Figure S4 . R 2 accuracy of regression models (order 2, without linear term) as function of delay ∆ (measured in days) for 6 different categories of GMI. Figure S7 . Country-specific regression models (14 figures for 14 countries). The description is the same as for Figure 3 . 