Evaluating the impacts of non-pharmaceutical interventions on the transmission dynamics of COVID-19 in Canada based on mobile network

The COVID-19 outbreak has caused two waves and spread to more than 90% of Canada’s provinces since it was first reported more than a year ago. During the COVID-19 epidemic, Canadian provinces have implemented many Non-Pharmaceutical Interventions (NPIs). However, the spread of the COVID-19 epidemic continues due to the complex dynamics of human mobility. We develop a meta-population network model to study the transmission dynamics of COVID-19. The model takes into account the heterogeneity of mitigation strategies in different provinces of Canada, such as the timing of implementing NPIs, the human mobility in retail and recreation, grocery and pharmacy, parks, transit stations, workplaces, and residences due to work and recreation. To determine which activity is most closely related to the dynamics of COVID-19, we use the cross-correlation analysis to find that the positive correlation is the highest between the mobility data of parks and the weekly number of confirmed COVID-19 from February 15 to December 13, 2020. The average effective reproduction numbers in nine Canadian provinces are all greater than one during the time period, and NPIs have little impact on the dynamics of COVID-19 epidemics in Ontario and Saskatchewan. After November 20, 2020, the average infection probability in Alberta became the highest since the start of the COVID-19 epidemic in Canada. We also observe that human activities around residences do not contribute much to the spread of the COVID-19 epidemic. The simulation results indicate that social distancing and constricting human mobility is effective in mitigating COVID-19 transmission in Canada. Our findings can provide guidance for public health authorities in projecting the effectiveness of future NPIs.


Introduction
Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) caused coronavirus disease 2019 . The COVID-19 has been declared as a public health emergency by the

Methods
In this section, we introduce the sources of data and present the network model that was developed to evaluate the impact of human mobility and NPIs on the spread of the COVID-19 epidemic in Canada.

Data collection and analysis
To capture the COVID-19 transmission dynamics in Alberta (AB), British Columbia (BC), Manitoba (MB), New Brunswick (NB), Newfoundland and Labrador (NL), Nova Scotia (NS), Ontario (ON), Quebec (QC), and Saskatchewan (SK), we fit publicly available incidence data of COVID-19 for Canada to the mathematical model incorporating human mobility. We obtain the number of COVID-19 cases in nine Canadian provinces from Johns Hopkins University Center for Systems Science and Engineering [12] (see Fig A.1 of S1 File). The mobility data is obtained from Google COVID-19 Community Mobility Reports (CMR) at https:// google.com/covid19/mobility [13] as is plotted in Fig A.2 of S1 File. The areas under study are divided into six different categories, which can be summarised as retail and recreation, grocery and pharmacy, parks, transit stations, workplaces, and residences. CMR show how visits and duration of stay at types of different places change around the baseline. The baseline is the median value for the corresponding day of the week during the five weeks from January 3 to February 6, 2020 [13]. We use the average of the previous seven days to approximate the missing values. Cross-correlation is a standard method of estimating the correlation between different sets of data (see Appendix B of S1 File). To determine which factor impacts the dynamics of COVID-19 most, we use the cross-correlation analysis [14][15][16][17][18] to derive the correlation between the Google Community Mobility Reports data and the weekly number of confirmed cases from February 15 to December 13, 2020. We also evaluate the impact of the time-lag of mobility data with respect to the incidence data on the COVID-19 outbreak.
The cross-correlation between the weekly number of confirmed cases and the Google Community Mobility Reports data for the nine Canadian provinces under study are as shown in Fig 1. The correlation between the weekly number of confirmed cases and the Google Community Mobility Reports data is statistically significant, and the cross-correlation between the weekly number of confirmed cases and human mobility data in residential areas are different from that between the weekly number of confirmed cases and human mobility data in other five classes of places (see Fig 1). Specifically, the human mobility in residential areas and the weekly number of confirmed cases are negatively correlated when the average time-lag in each province is four weeks, while the human mobility in the other five classes of places and the weekly number of confirmed cases are positively correlated when the average time-lag in each province was four weeks. The above results show that activities around the residences do not The variables x 1 (A), x 2 (B), x 3 (C), x 4 (D), x 5 (E), and x 6 (F) denote retail and recreation, grocery and pharmacy, parks, transit stations, workplaces, and residences. The gray area represents the 95% confidence intervals.
https://doi.org/10.1371/journal.pone.0261424.g001 accelerate the spread of the COVID-19 epidemic. The average time-lag between human mobility in parks (eight weeks) and the weekly number of confirmed cases is longer than that between human mobility in retail and recreation (four weeks), grocery and pharmacy (six weeks), transit stations (six weeks), and workplaces (one week) during the presence of positive correlation. The cross-correlation analysis between Google Community Mobility Reports data and the weekly number of confirmed death cases are shown in Fig A.3 of S1 File.

Model formulation
In order to mimic the spread of the COVID-19 epidemic, we construct a disease transmission model with a population mobile network. The whole population is divided into nine classes, namely, S(t), E(t), A(t), I(t), Q(t), H(t), H i (t), R(t), and D(t) that represent the numbers of individuals who were: (i) susceptible; (ii) exposed; (iii) infectious (asymptomatic); (iv) infectious (symptomatic); (v) isolated symptomatic and infected (mild to moderate); (vi) hospitalized (severe) not requiring ICU care; (vii) hospitalized (severe) requiring ICU care; (viii) recovered; and (ix) who had died, respectively. The total population at time t is denoted by N(t). The population flow among those compartments is shown in Fig 2. The population mobile network consisted of nodes representing six categories of places, namely, retail and recreation, grocery and pharmacy, parks, transit stations, workplaces, and residences. Since the human mobility of individuals in each place is different during the day and night [19], we develop different models for the daytime and the night. During the daytime, the human mobility mainly occurs in retail and recreation, grocery and pharmacy, parks, transit stations, workplaces, and residences. During the night, people usually stay in residences. The model for the daytime is For Models (1) Table 1. Here, λ j (t) is the rate of infection in place j at time t, where N j (t)(j = 1, 2, 3, 4, 5, 6) represents the number of people in retail and recreation, grocery and pharmacy, parks, transit stations, workplaces, and residences at time t,  [20] θ j The factor for reduced transmission probability of asymptomatic infected individuals (dimensionless) 0.55 [21,22] δ j The factor for reduced transmission probability of exposed individuals (dimensionless) 0.55 [21,22] ρ respectively. Therefore, the rate of infection in places j at time t can be expressed as Thus, the rate of infection in places j at time t can be rewritten as During the COVID-19 epidemic, the NPIs taken by the Canadian public health authorities [6] are shown in Fig 3. The specific NPIs are as follows: T 1 (Mar 1, 2020): Entertainment/cultural sector closure. T 2 (Mar 13, 2020): School closure. T 3 (Apr 6, 2020): Mask wearing. T 4 (Apr 21, 2020): Imposed entertainment/cultural sector closure. T 5 (Jul 7, 2020): Imposed entertainment/cultural sector closure. T 6 (Sep 1, 2020): Imposed entertainment/cultural sector closure. T 7 (Nov 2, 2020): Imposed entertainment/cultural sector closure. Therefore, the transmission rate can be expressed as the piecewise function The effective reproduction number R ej (t) in places j at time t is The average effective reproduction number R e (t) at time t is where the three terms represent the numbers of new cases generated by exposed, asymptomatic, and symptomatic individuals, respectively. The probability of a susceptible individual being infected in places j at time t is defined as Similarly, the model for nighttime is

Parameter estimation
We fit Models (1) and (2) to the weekly number of confirmed cases to quantify the dynamics of COVID-19 in nine Canadian provinces. We simulate the epidemic for each province from February 15, 2020. Since daily numbers of new cases were not reported in Canada on February 15, 2020, we assume that the initially daily number of newly infected individuals on February 15, 2020 is zero in simulations. The size of total population and the age pyramid for the nine provinces are obtained from Statistics Canada [31]. In the simulations, we assume that the population is spatially uniformly distributed before February 15, 2020. After February 15, 2020, we use the Google Community Mobility Reports data to quantify the number of people in six types of places (retail and recreation, grocery and pharmacy, parks, transit stations, workplaces, and residences). Since the incubation period of COVID-19 is around 5.2 days [20], the incubation rate σ = 1/5.2. Regarding the survival rate and length of ICU stay duo to COVID-19 infection, we assume that the average length of ICU stay and the average time from hospitalization to deaths are both 21 days, and the average length of hospital stay rather than ICU is ten days, and the duration of quarantine for exposed cases is 14 days [27], that is, 1/γ Hi = 1/d H = 21, 1/γ H = 10, and 1/γ Q = 14. Moreover, we assume that the asymptomatic infectious period is eight days [25]. Thus, we let γ A = 1/8 per day. We assume that the average time to case detection is six days [24], that is, 1/� = 6. Around 30% − 60% of people infected with due to COVID-19 are asymptomatic or only have mild symptoms, and the transmissibility is lower, but still significant [23]. Thus, we assume that the probability of asymptomatic infected individuals is ρ = 60%. We set θ j = δ j = 0.55 [21,22] by assuming that exposed and asymptomatic infected individuals have lower transmissibility. We additionally assume that the probability of home isolation is 84.3% and the probability of ICU admission is 10% [24,[28][29][30], that is, 1 − α = 84.3% and 1 − P QH = 10%.
In order to simulate the number of new cases per week in nine Canadian provinces, the feasibility of the model is verified by the actual number of newly infected cases. Next, we use the MCMC method [32,33] for 50000 iterations with a burn-in of 45000 iterations to fit the Models (1) and (2). We estimate the unknown parameters and initial conditions for the Models (1) and (2) ; P HR Þ), using the weekly number of confirmed cases in Canada. Let Cðt;ŵÞ represents the cumulative number of confirmed cases. The dynamics of Cðt;ŵÞ is dCðt;ŵÞ dt ¼ �I: The weekly number of confirmed cases is where P C represents the weekly number of confirmed cases, and the time step is hours in the simulations, then C(0) = P C (0). We have C independent observations from the data, representing the number of confirmed cases on the i-th week, where i = 1, . . ., C. Let ε be the fitting error, and ε follows the additive independent Gaussian distribution with mean zero and unknown variance ξ 2 . Thus, the observations � Y C can be expressed as follows For simplicity, we assume that the unknown parametersŵ of system is an independent Gaussian prior specification. Therefore, where M is the number of unknown parameters. We further assume that the inverse of the error variance follows a gamma distribution as prior with the form where S 2 0 and n 0 are the prior mean and prior accuracy of variance ξ 2 , respectively.
The likelihood function pð � Y C jŵ; x 2 Þ for C independent identically distributed observations from Eq (4) with a Gaussian error model is where SSðŵÞ represents the sum of squares function The conditional distribution pðx À 2 j � Y C ;ŵÞ can be expressed as follows According to the conditional conjugacy property of the Gamma distribution [34,35], the conditional distribution pðx À 2 j � Y C ;ŵÞ is also a Gamma distribution with based on which sample and update ξ −2 for other parameters within each run of Metropolis Hastings simulation steps. Since we assume independent Gaussian prior specification for parametersŵ, the prior sum of squares for the given parametersŵ can be calculated as follows Hence, for a fixed value of variance ξ 2 , the posterior distribution of parametersŵ can be expressed as follows and the posterior ratio needed in the Metropolis-Hastings acceptance probability can be written as follows Accordingly, the new unknown parameter valueŵ 2 will be accepted with probability We find the local minimum of SSðŵÞ using nonlinear least square approach to obtain the estimated value of each sample. The minimum value of these estimatedŵ is used as the initial guess in MCMC simulation. Prior information of unknown parameters is given by E(0)2(0, 1000), a ji 2 (0, 10)(j = 1, 2, � � �, 6, i = 1, 2, � � �, 8), P HR 2 (0, 1), and the proposal density follows a multivariate normal distribution. By fitting Models (1) and (2)

Results
In this section, we calculated the effective reproduction numbers, transmission rates, infection probabilities, and quantified the impact of human mobility on COVID-19 in nine Canadian provinces.

The impact of mobility on COVID-19
Mobility in Canada dropped sharply in March 2020: e.g., in March and April, 2020, the average percentage of retail and recreation visits changed from the baseline by -37.5%, the average percentage of grocery and pharmacy visits changed from the baseline by -12.6%, the average percentage of park visits changed from the baseline by -2.3%, the average percentage of transit station visits changed from the baseline by -49.0%, the average percentage of workplace visits changed from the baseline by -39.8%, and the average percentage of human mobility in residential areas changed from the baseline by 15.8%. However, after April, 2020, the number of people going to the park increased significantly. From May to November, the average percentage of park visits increased by 100.1% from the baseline.
In this regard, we simulate the spread of COVID-19 when the average percentage of park visits changed from the baseline after April 25, 2020 (see Fig 6). We observe that when the average percentage of park visits changed from the baseline after April 25, 2020, the cumulative number of cases in Canada as of December 13, 2020 is reduced by 32.97%. More predictions for the impact of reducing park visits on the cumulative number of cases in nine Canadian provinces are summarized in Table 2. In addition, we simulate the scenario when the movements of humans are reduced by 20%, as of December 13, 2020, the cumulative number of cases in Canada is reduced by 56.72%. More predictions for the impact of reducing population movement on the cumulative number of cases in nine Canadian provinces are summarized in Table 2.

The transmission rate in types of different places
In order to study the spread of the COVID-19 epidemic in Canada, we plot the transmission rates in types of different places of the nine Canadian provinces from February 15, 2020 to December 13, 2020, as shown in Fig 7. Table 3 shows that the transmission rate changes in different phases are mostly between -10% and 10% in Ontario and Saskatchewan, which shows that the effect of NPIs is not obvious in these two province. However, in Alberta, British Columbia, Manitoba, New Brunswick, Newfoundland and Labrador, Nova Scotia, and Quebec, the transmission rates have changed dramatically at different phases.
We investigate where the transmission rate (β j (t)) is the most important. We calculate the transmission rate (β j (t)) of each place in the simulations from February 15 to December 13, 2020 (see Fig 7), and we find that the average transmission rates over the entire simulation period are relatively large in the grocery and pharmacy, and transit stations, which are 0.2584 and 0.2587, respectively, due to the high population density in the grocery and pharmacy, and transit stations. The average transmission rates over the entire simulation period are relatively small in the workplaces and parks, which are 0.2370 and 0.2388, respectively, because most people stay home during the epidemic, thereby reducing the population density of the workplaces. In addition, the average transmission rates of the retail and recreation, and residences during the entire simulation period are 0.2464 and 0.2520, respectively.

The infection probability in types of different places
Based on the estimated parameters in Table C.1 of S1 File, and the value of variables E, A, I and N obtained by the model, we calculate the probability of susceptible individuals being infected in types of different places in nine Canadian provinces, as shown in Fig 8. We observe that the trend of the daily infection probability is consistent with the number of daily confirmed cases. As of December 13, 2020, in Alberta, British Columbia, Manitoba, New Brunswick, Newfoundland and Labrador, Nova Scotia, Ontario, Quebec, and Saskatchewan, the   The first measure represents the proportion of people who visited parks after April 25, 2020 as the baseline, and the second measure represents a 20% reduction in the number of population movements.
the average infection probability in Quebec is higher than those of the other eight provinces from mid-March 2020 until mid-October 2020. After mid-October 2020, the average infection probability in Manitoba exceeded the average infection probability in Quebec. After November 20, 2020, the average infection probability in Alberta has increased very rapidly and exceeded the average infection probability in Manitoba, and the average infection probability in Alberta has reached the highest value since the onset of the COVID-19 epidemic.

Uncertainty and sensitivity analysis
As there are lots of the parameters involved in the model, variations of the parameters may affect the trend of the epidemic greatly. To study the sensitivity of parameters on the model prediction results, we conduct Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCCs) analysis [36]. The goal is to identify most important parameters affecting the average effective reproduction number R e (t). The input parameters are the average transmission rate (i.e. � b), the factor for reduced transmission probability of asymptomatic infected individuals (i.e. θ), the factor for reduced transmission probability of exposed individuals (i.e. δ), the mean length of incubation period (i.e. 1/σ), the asymptomatic infectious period (i.e. 1/ γ A ), the average time to case detection (i.e. 1/�) and the proportion of asymptomatic infected individuals (i.e. ρ), the output variables is the average effective reproduction number R e (t). We choose that all input parameters are normally distributed. We consider absolute values of PRCC greater than 0.4 as indicating the high correlation between input parameters and output variables, values between 0.2 and 0.4 as moderate correlations, and values between 0 and 0.2 are not significantly different from zero [37,38]. The results of the sensitivity analysis of parameters are shown in Fig 10.  Fig 10 illustrates the PRCCs for the dependence of the average effective reproduction number in nine Canadian provinces on each parameter form February 15, 2020 to December 13, 2020. We found that the parameters � b, θ and δ are highly positively correlated with the average effective reproduction number, and the correlation coefficient between � b and the average effective reproduction number is the highest, which indicates that reducing parameter � b has a greater impact on the COVID-19 epidemic than parameters θ and δ. The parameters σ, γ A and � are highly negatively correlated with the average effective reproduction number, and the correlation coefficient between γ A and the average effective reproduction number is the lowest, which indicates that the early detection of asymptomatic infectious cases has a greater impact on the COVID-19 epidemic than parameters σ and �.

Discussion
The government of Canada has initiated a series of NPIs to contain the spread of the COVID-19 epidemic [6]. Public heath policies vary in space and time, resulting in varied patterns of movement and behavioural changes throughout the country. Due to the changes in movement patterns and behavioural changes of humans, the severity of the COVID-19 outbreak in Canada varies from place to place.
To understand how human mobility, the transmission rate, and infection probability affect the spread of the COVID-19 epidemic, we developed a meta-population model for the transmission dynamics of COVID-19, where the human mobility is characterized by a dynamic network. The model takes into account the heterogeneity between different provinces in Canada, such as the timing of implementing NPIs, the human mobility in retail and recreation, grocery and pharmacy, parks, transit stations, workplaces, and residences due to work, and recreation.
We only analyzed the epidemic in nine Canadian provinces from February 15, 2020 to December 13, 2020, including Alberta (AB), British Columbia (BC), Manitoba (MB), New Brunswick (NB), Newfoundland and Labrador (NL), Nova Scotia (NS), Ontario (ON), Quebec (QC), and Saskatchewan (SK), since the COVID-19 epidemic in other Canadian provinces is not serious due to low population densities.
As is known, the effect of human mobility on the spatial and temporal spread of COVID-19 is delayed a few weeks due to the incubation period of COVID-19 infection and reporting delay [39][40][41]. To determine which activity is most closely related to the dynamics of COVID-19, we used the cross-correlation analysis [14] to identify the places where the human mobility data are statistically significant and correlated with the weekly number of confirmed cases from February 15 to December 13, 2020. We found that activities in locations around the residences can reduce the severity of COVID-19, and when the time-lag is about 16 weeks, the positive correlation is the highest between the mobility data of parks and the weekly number of confirmed COVID-19 cases. Our model has incorporated many key factors responsible for the temporal and spatial spread of the COVID-19 epidemic, such as the timing of implementing NPIs and the human mobility. Since intra-province travel contributes more to the increase of cumulative number of cases than inter-province travel [9,11] due to substantially reduced long-distance travels, our model only considered intra-province human mobility. We found that the average effective reproduction numbers in nine Canadian provinces are all greater than one from February 15, 2020 to December 13, 2020. NPIs have little impact on the epidemics in Ontario and Saskatchewan. After November 20, 2020, the average infection probability of Alberta has reached the highest since the COVID-19 epidemic in Canada. We carried out the sensitivity analysis of the average effective reproduction number with respect to the model parameters. The result identified six key parameters: the average transmission rate (i.e. � b), the factor for reduced transmission probability of asymptomatic infected individuals (i.e. θ), the factor for reduced transmission probability of exposed individuals (i.e. δ), the mean length of incubation period (i.e. 1/σ), the asymptomatic infectious period (i.e. 1/γ A ), the average time to case detection (i.e. 1/�), play an important role in the size of the average effective reproduction number. This suggested that the transmission rate, the mean length of incubation period, the asymptomatic infectious period, and the average time to case detection play an important role during the COVID-19 epidemic. Our findings provide insights that help support the planning of emergency management and decision-making for public health authorities. However, our study still has several limitations. First, since the numbers of asymptomatic infected individuals and recovered individuals are not publicly available yet, our simulations only used incidence data and death cases. Second, behavior and personal risk management such as masks, no indoor gatherings, social distancing, compliance with quarantine measures are intertwined with movement and hard to disentangle from the movement data. In fact, we are likely seeing the results not of movement alone but of a combination of movement, NPI, and individual behavior change as a result of public health orders but also perceived risk and information from surrounding areas. Third, we did not consider the effectiveness of contact tracing and the speed of tracing, and the human mobility between provinces and the imported cases from overseas and the role of environmental factors, which will be studied in future work when such data become available.
Supporting information S1 File. Data collection and analysis, the details of the cross-correlation function, the mean values and standard deviation of parameters and initial for Models (1) and (2). (ZIP)