key: cord-0890969-ddw2vojc
authors: Franco, C.; Ferreira, L. S.; Sudbrack, V.; Borges, M. E.; Poloni, S.; Prado, P. I.; White, L. J.; Aguas, R.; Kraenkel, R. A.; Coutinho, R. M.
title: Percolation across households in mechanistic models of non-pharmaceutical interventions in SARS-CoV-2 disease dynamics
date: 2021-06-09
journal: nan
DOI: 10.1101/2021.06.07.21258403
sha: 605c9f5e1ea395718426b54782ef8449b096dfdb
doc_id: 890969
cord_uid: ddw2vojc

Since the emergence of the novel coronavirus disease, mathematical modelling has become an important tool for planning strategies to combat the pandemic by supporting decision-making and public policies, as well as allowing an assessment of the effect of different intervention scenarios. A proliferation of compartmental models was observed in the mathematical modelling community, aiming to understand and make predictions regarding the spread of COVID-19. Such approach has its own advantages and challenges: while compartmental models are suitable to simulate large populations, the underlying well-mixed population assumption might be problematic when considering non-pharmaceutical interventions (NPIs) which strongly affect the connectivity between individuals in the population. Here we propose a correction to an extended age-structured SEIR framework with dynamic transmission modelled using contact matrices for different settings in Brazil. By assuming that the mitigation strategies for COVID-19 affect the connections between different households, network percolation theory predicts that the connectivity across all households decreases drastically above a certain threshold of removed connections. We incorporated this emergent effect at population level by modulating the home contact matrices through a percolation correction function, with the few remaining parameters fitted to to hospitalisation and mortality data from the city of Sao Paulo. We found significant support for the model with implemented percolation effect using the Akaike Information Criteria (AIC). Besides better agreement to data, this improvement also allows for a more reliable assessment of the impact of NPIs on the epidemiological dynamics.

Since the emergence of the novel coronavirus disease, mathematical modelling has become an important tool for planning strategies to combat the pandemic by supporting decision-making and public policies, as well as allowing an assessment of the effect of different intervention scenarios. A proliferation of compartmental models was observed in the mathematical modelling community, aiming to understand and make predictions regarding the spread of COVID-19. Such approach has its own advantages and challenges: while compartmental models are suitable to simulate large populations, the underlying well-mixed population assumption might be problematic when considering non-pharmaceutical interventions (NPIs) which strongly affect the connectivity between individuals in the population. Here we propose a correction to an extended age-structured SEIR framework with dynamic

Since the emergence of SARS-CoV-2, mathematical modelling has become an important tool for planning strategies to combat the pandemic by supporting decision-making and public policies [1, 2, 4, 22] . As a result, a proliferation of SEIR-like age-structured compartmental models assessing 5 the effect of different intervention scenarios have been proposed and discussed [12, 16, 17] . As simplified representations of reality, these models come with assumptions regarding the nature of the underlying network of interactions [25] . More explicitly, most SEIR compartmental models assume that the population is homogeneously-mixed, meaning that every individual in each 10 compartment has the same probability of contacting others, regardless of spatial distribution.

Though these underlying assumptions might be appropriate to make robust predictions for well-mixed populations, they may not be suitable in a context with significant contact network heterogeneity. Therefore, several 15 methods have been proposed to account for contact structure heterogeneities in compartmental models [20, 26, 28] . Nevertheless, they usually consider static underlying network structures, which are only good approximations for slow epidemic spreads, in which contact rates do not change significantly over time [30, 5] . This is not the case in a pandemic such as the one caused connectivity between individuals. The need for a more flexible and dynamic approach became evident, hence motivating the present study.

In this context, network theory provides a picture of a homogeneously 25 mixed population as a (highly connected) regular random network of individuals (vertices) connected through possible disease transmission contacts (edges) with long-range connectivity properties [5] . During an outbreak, the disease would spread across these links. Social distancing nonpharmaceutical interventions (NPIs) could, therefore, be thought of as mod-30 ulators of the strength and even persistence of such links. One of the main results of network theory is that, as contact networks become less connected, a critical transition occurs, and the network becomes fragmented into disconnected (or weakly connected) sub-networks. This phenomenon is known as percolation [7, 14] , and the existence of a critical percolation threshold 35 in the mean number of contacts is established for many kinds of networks. Across this threshold, even small changes in the number of contacts can lead to large changes in the connectivity of the network. In compartmental models, SARS-CoV-2 mitigation strategies are modelled by altering contact rates, thus changing the force of infection. This can 40 be done at different degrees of granularity, depending on the level of detail of contact matrices. For instance, many works [2, 21, 11] have used the categorisation employed by Prem et al. [24] , which divides contacts into four settings, namely home, work, school, and other. In this way, the effectiveness of different NPIs are reflected in reductions in the contact rates for each setting, 45 depending on the nature of the intervention -e.g. school closure reduces contact rates at school. Since the contact matrix is a linear combination of the contact matrices for each setting, with coefficients dependent on the NPIs' coverages and efficacies, each NPI contributes linearly to reduce the infection force. By itself, this change in contacts among individuals does not 50 affect the relationship between the mean number of contacts in the underlying network and the force of infection of the compartmental model. This is adequate if the structure of the network is not greatly affected. However, when social distancing NPIs are applied with high coverages and connections between individuals are continuously being removed, this approximation is 55 prone to break down, and the compartmental model may no longer provide a satisfactory portrayal of the epidemiological dynamics.

To account for this change in network structure and resultant reduction in the force of infection, here we propose to modify how the effect of NPIs are modelled in compartmental models. We integrate results from percolation 60 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) theory into an age-structured SEIR model by using a non-linear correction between the strength of NPIs and the resulting contact matrix, thus directly affecting the force of infection. In the following sections, we build on the previously implemented and widely used COVID-19 Modelling Consortium (CoMo) model [2], adapting its compartmentalisation to the Brazilian hospital system. Then, we present a nested model version that takes into account the loss of long-range connectivity (percolation) effect. Finally, we fit both model versions to hospitalisation and mortality data for SARS-CoV-2 during 2020 in São Paulo, Brazil, the country's most populous city, with over 12 million inhabitants, and the first to detect COVID-19 cases. Through model comparison using the Akaike Information Criteria (AIC) [6] we find that the data strongly supports the model incorporating percolation effect.

To model the epidemiological dynamics of COVID-19 in São Paulo, we ap-75 ply an age-structured SEIR model with infected compartments stratified by severity of symptoms, treatment requirement and accessibility to healthcare. The main framework for this SARS-CoV-2 epidemic model was developed in collaboration with the CoMo Consortium [2] with slightly different treatment seeking compartments, an adaptation applied to better represent the Brazil-80 ian organisation of hospital beds. The progression of individuals through the infection cycle, for this version of the CoMo model, is represented by the diagram in Figure 1 . Each model compartment depicted in Figure 1 is divided into 19 subcompartments, comprising all 19 age classifications from the Brazilian Insti-85 tute of Geography and Statistics (IBGE) [19] . Transmission between each pair of age classes is therefore evaluated based on the contact rates between them, which are the contact matrices estimated for different settings (school, work, home, and others) in Brazil by Prem et al. [24] .

We considered non-pharmaceutical intervention (NPI) scenarios including 90 social distancing, work from home, and school closure. Here, the effect of mitigation strategies enter as linear corrections to the contact matrices, since these interventions aim to reduce contact rates or the risk of contagious at each possible contact (mask use, for example). For instance, a certain coverage on the school closure NPI would decrease contacts between all age 95 4 . 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 preprint this version posted June 9, 2021. ; Figure 1 : A diagram of the baseline model structure for SARS-CoV-2 in Brazil, representing the unmitigated epidemic spread scenario. The variables in the compartments represent individuals S: susceptible, E: presymptomatic infectious, A: asymptomatic infectious, I: infectious with mild symptoms, X: infectious with mild symptoms and self-isolated, H C : infectious, requiring hospital treatment but denied, H: infectious and hospitalised, ICU: infectious, receiving intensive care, ICU C : infectious, requiring intensive care but denied, ICU H : infectious, requiring intensive care but being admitted in regular hospital bed, R: Recovered, D: Deceased. All compartments are further subdivided into 5-year age classes from 0 to 90+ years old.

. 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 preprint this version posted June 9, 2021. classes in the school setting, which is modelled as a proportional decrease in contact rates in school contact matrices (see SM).

Though for school, work, and other environments such linear relationship between coverage of mitigating strategies and decrease in contact rates might be an appropriate assumption, we expect that the response should be 100 different in the household setting. The reason for this difference is that compartmental models such as SIR-like models assume a homogeneously mixed population in all settings, which is not a good approximation for the dynamics of infection in a network of households under social distancing, as we will discuss in the next session. 

The model described in the previous session introduced a phenomenological implementation to model the effect through a proportional reduction of contact rates in the matrices. However, this approach fails to substantially decrease the transmission even when considering high coverages of NPIs. For 110 instance, at hypothetical coverages and efficacies of 100% for school closure, work from home, and social distancing, the infection rate would still not go down immediately. Consequently, we found that data fitting was challenging when trying to simulate NPI coverages as closely as possible to the patterns seen in real life.

In a non-NPI scenario, we could picture the whole population as a collection of households assigned as vertices with links between households (created due to interactions at work, schools, or other places) represented by edges. This highly connected network would form one giant connected component, which could therefore be approximated as a homogeneously mixed population 130 [5] . However, by introducing social distancing measures and increasing their coverages, we break connections in the network until long-range connectivity 6 . 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 preprint this version posted June 9, 2021. ; is lost. The critical value of number of connections (or edge density) where this transition happens is the so-called percolation threshold [14] . Above the percolation threshold of contact loss, i.e., in a high coverage social distancing 135 scenario, we would expect to see the emergence of small household clusters that, though well connected within themselves, would be poorly connected between them.

Here we model the percolation effect on home contact matrices assuming that, while contacts outside households are kept above the percolation 140 threshold, the effect of mitigation strategies is less apparent in home settings. On the other hand, the probability of infection of a susceptible individual drops drastically when its mean number of contacts goes below a threshold [23] . Mathematically, we correct the home contact matrix by a factor dependent on all NPIs' coverages and efficacies, which we write as a percolation 145 correction function, f perc , that has to satisfy several requirements:

• 0 ≤ f perc ≤ 1;

• f perc → 0 as NPIs are completely lifted; for low adherence to NPI, no connectivity loss should be noticed since different households would still be strongly connected;

• 1 − f perc 1 for high effectiveness of NPIs; as connections between different households are widely severed, so f perc approaches its maximum value.

A hyperbolic functional form, as follows, is proposed to model this effect:

where 0 < h ef f < 1 is the maximum reduction in contacts, a limit introduced due to the percolation effect; 0 < T perc < 1 is the percolation threshold, i.e., the fraction of connections we need to remove so that the network no longer percolates; and h steep is the steepness of the percolation correction function determining how fast the network percolates: a large h steep implies f perc changing abruptly from 0 to h ef f when W N P I ≈ T perc (see Figure 2 ). Finally, W N P I is defined as the weighted non-pharmaceutical interventions, combining the efficacy and coverage of each type of intervention, as well as the age and contact structure of the population they're applied to, namely W N P I = p work a work (t) + p school a school (t) + p other a dist (t)

7 . 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 preprint this version posted June 9, 2021. where, for instance, work adherence (a work ) can be written as the product of the respective effectiveness (work ef f ), coverage (work cov ) and duration of home-office strategies (translated as the step function θ work (t)):

The same holds for school and dist, forming the set {a j (t)}, j = {work, school, dist}, which throughout the text will be referenced as adherences. The weights p j , j = {work, school, other}, are calculated as

where P(t = 0) is the initial age distribution of the population and c j the contact matrices in each setting. Here we used T perc = 0.6 for the percolation threshold, which is indicated as the curve's inflexion point. Note that changing its value simply dislocates the curve along the x-axis. For low W N P I , f perc → 0, meaning no connectivity is lost at low adherence levels. Near W N P I = T perc , f perc grows rapidly. At the high W N P I regime, where more connections are lost, f perc → h ef f ≈ 1.

. 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 preprint this version posted June 9, 2021. ; https://doi.org/10.1101/2021.06.07.21258403 doi: medRxiv preprint For each NPI (j = {work, school, dist}, i.e., {work from home, school closure and social distancing}), we parametrize the respective adherence, a j , with efficacy and coverage parameters, here modelled as a combination of step functions allowing for changing values over time [18, 29] . Now, clearly 0 < W N P I < 1. Note also that W N P I has this specific form 160 to consider population structure, accounting for how much each of these NPIs effectively reduce contacts in each age-class. Also note that W N P I can vary according to the current implemented interventions, adding flexibility to the model. Finally, considering all effects due to NPIs and percolation, the effective 165 contact matrix,ĉ, is written aŝ

In the simulations considered here, the populations and contact matrices are obtained from publicly available data [24] ; all a j are estimated based on the mitigation strategies adopted by the city of São Paulo; whereas h steep and T perc are obtained through fitting to epidemiological data. 

The data used were time series of hospitalisations and deaths of Severe Acute Respiratory Illness (SARI) in São Paulo, Brazil, reported in the 2020 SIVEP-Gripe database [9] . To be sure records were already consolidated, we restricted the data to weekly aggregated data points comprising the first 23 175 weeks of the pandemic in the city, from 15 March to 31 August 2020.

Based on the current literature and proxy data for mobility and coverage of NPIs [18, 29, 27], we set values for all parameters with reasonable data sources or proxies (see SM for all parameter values and references). Thus only a few free parameters could not be inferred from literature and, hence, 180 needed to be fitted to the epidemiological data. These fitted parameters were the probability of infection given a contact (p), start date of community transmission (startdate), and both h steep and T perc from Equation 1 for the model that takes percolation into account.

. 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 preprint this version posted June 9, 2021. ; https://doi.org/10.1101/2021.06.07.21258403 doi: medRxiv preprint

We used a Levenberg-Marquardt nonlinear least square regression fitting algorithm [13] to minimise the squared residuals, which is equivalent to maximising the likelihood under the assumption that the data points errors are normally distributed (details in SM). With that, we obtained the most probable set of values for p, startdate, h steep and T perc . This algorithm was applied • Standard model + 30% NPIs: the model and parameters are the 195 same as before, but we consider a much higher (30%) adherence to all NPIs. This was done to ensure that a subestimation of the effect of NPIs would not by itself favour the model with percolation -given that we already know that it always reduces the contact rates. However, notice that this leads to unrealistically high values of adherences.

• Model with percolation: derived from standard model, including the percolation effect, where T perc , h steep , p and startdate were fitted to the data.

Since the model with percolation introduces two extra parameters fitted to data into the standard model, we can consider them nested models and, 205 hence, compare the quality of the fittings using the AIC.

We can see the effect of implementing the percolation correction into our model by calculating the resultant basic reproduction number (R 0 ) considering how it changes in an environment subject to NPIs. R 0 is calculated 210 using the Next Generation Matrix (NGM) method [3, Chapter 6] . In Figure 3 , we see that the calculated R 0 diverges among the two model versions as the combined NPI adherence approaches the fitted percolation threshold (indicated as T perc ). Percolation implies a significantly lower R 0 values for higher NPI adherence, which is consistent with the steep decrease in NPI 215 effect near the percolation threshold modelled by the percolation correction function (Eq. 1).

. 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 preprint this version posted June 9, 2021. ; After implementing both model versions with equivalent fixed parameter sets, and also with increased NPIs in one scenario (see complete table of parameters in the SM), we fitted the free parameters (among T perc , h steep , p 220 and startdate) in each model to weekly new cases and new deaths from 2020 SARI data for São Paulo. Figure 4 shows the resulting curves for all fitted models, with confidence intervals assuming that fitted parameters follow a multivariate normal distribution. These confidence intervals are estimated by bootstrapping the co-225 variance matrix from the Levenberg-Marquardt algorithm. The comparison between the two standard model curves shows that increasing intervention coverages (even to the point of assuming unrealistic values) can modulate the epidemic curve to a certain extent, but still does not result in a good fit to data.

The resulting parameters and the computed AIC for each model version are shown in Table 1 . Even though in Figure 4 the curves for the standard model + 30% NPI and the model with percolation might look somewhat close, the AIC clearly shows a much higher level of empirical support for the model version with percolation. 235 

. 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. 

Modelling the ongoing COVID-19 pandemic represented an unprecedented challenge for modellers, as we tried to make urgent recommendations and predictions with scarce information about the virus's underlying biology, its main spread mechanisms, and severity and fatality rates [10, 8] . Furthermore, 240 we lacked individual level data on mobility patterns to accurately understand the dynamic effect of social distancing NPIs on individual behaviour and contact patterns.

Although previous works considered heterogeneous contact structures in epidemiological models [20, 26, 28] , and many others accounted for NPIs' 245 effects on transmissibility [12, 16, 17, 1, 2, 4, 22] , few to none made the 12 . 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 preprint this version posted June 9, 2021. ;  connection between both phenomena. Strictly speaking, NPIs can affect the underlying network structure over time due to large scale social distancing interventions. In this context, we identified a demand for more flexible and dynamic modelling approaches taking into account the non-linear phenomena 250 emerging from social distancing NPIs.

In the case of the COVID-19 pandemic, governments enforced wide movement and contact restrictions, essentially removing connections between individuals. To take into account the structural heterogeneities introduced by these social distancing measures on a population's contact network, we pro-255 posed a phenomenological correction to compartmental models that only depends on quantities already computed when modelling NPIs, without the addition of an excessive number of parameters which could lead to model overfitting. Using the Akaike Information Criterion (AIC) we found ∆AIC = 191 between the model with percolation and the standard model (Table 1) , which 260 is orders of magnitude higher than the conventional threshold of ∆AIC ≥ 2 to consider significant support for one model [15] , that is, the model with percolation correction had very strong empirical support with the data analysed.

One alternative explanation could be that it might still be possible to 265 adjust other model parameters to obtain a better fit to data in the standard model by adjusting NPIs' coverages. However, with further exploration of the parameter space, even considering less realistic parameter sets (e.g. tampering with NPIs' coverage parameters), we could not obtain a good fit to data. For instance, a 30% increase on NPI coverage on the standard model 270 (red curve in Figure 4 ) resulted in ∆AIC = 154. Though it represented a better fitting compared to the standard model version, it was not as good as the one obtained with the version including the percolation. These results highlight the importance when using compartmental models of implementing this non-linear response to social distancing NPIs in order 275 to obtain a better representation of its effect on a large population.

The correction to an age-structured SEIR compartmental model described here, though motivated by the need to better represent the COVID-19 spread dynamics, was inspired by network theory and general individual level con-280 siderations that affected a population's underlying network structure. We implemented the effect of individual behaviour changes on a population's 13 . 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)

The copyright holder for this preprint this version posted June 9, 2021. ; macroscopic dynamics without drastically increasing the number of fitted parameters in our compartmental model.

The usefulness of incorporating this fragmentation process into our SARS-

CoV-2 spread compartmental model was evident, but note that its flexibility also permits it to be applied when modelling the transmission dynamics of other communicable diseases under similar high NPIs' coverage regimes. Therefore, our framework may be applied to any homogeneously-mixing compartmental models trying to represent the dynamics of a population suf-290 fering drastic changes in its connectivity patterns during an epidemic. This result contributes towards a more accurate epidemic modelling, potentially implying in better control and prevention policy recommendations at public health level. 

. 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)

The copyright holder for this preprint this version posted June 9, 2021. ; https://doi.org/10.1101/2021.06.07.21258403 doi: medRxiv preprint

Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brazil (CAPES) (Finance Code 001 to LSF) and the Brazilian National Council for Scientific and Technological Development (CNPq) (grant number: 315854/2020-0 to MEB, 310 313055/2020-3 to PIP and 311832/2017-2 to RAK). RA is funded by the Bill and Melinda Gates Foundation (OPP1193472). LW is funded by the Li Ka Shing Foundation

Special report: The simulations driving the world's response to COVID-19

Modelling the COVID-19 pandemic in context: an international participatory approach

Mathematical epidemiology

Mathematical modeling and the transmission dynamics in predicting the Covid-19 -What next in combating the pandemic

When individual behaviour matters: homogeneous and network models in epidemiology

Model Selection and multi-model inference: A practical information-theoretic approach

Profile of a killer: the complex biology powering the coronavirus pandemic

Srag 2020 -banco de dados de síndrome respiratória aguda grave -incluindo dados da covid-19

Covid-19's known unknowns

Age-dependent effects in the transmission and control of COVID-19 epidemics

Effects of non-pharmaceutical interventions on covid-19 cases, deaths, and demand for hospital services in the UK: a modelling study

Algorithm Found in MINPACK, Plus Support for Bounds

Percolation theory

Appendix E: Model Selection Criterion: AIC and BIC

Impact of non-pharmaceutical interventions (NPIs) to reduce covid19 mortality and healthcare demand

Estimating the effects of non-pharmaceutical interventions on COVID-19 in europe

Covid-19 mobility reports

Dynamical behavior of epidemiological models with nonlinear incidence rates

Covid-19 scenarios: an interactive tool to explore the spread and associated morbidity and mortality of sars-cov-2. medRxiv

Mathematical modeling as a tool for policy decision making: Applications to the COVID-19 pandemic

Epidemic processes in complex networks

Projecting contact matrices in 177 geographical regions: an update and comparison with empirical data for the COVID-410 19 era. medRxiv

Heterogeneity and network structure in the dynamics of diffusion: Comparing agent-based and differential equation models

On representing network heterogeneities in the incidence rate of simple epidemic models

Semi-empirical power-law scaling of new infection rate to model epidemic dynamics with inhomogeneous mixing

Medidas de distanciamento social e evolução da covid-19 no brasil

We thank all members of Observatório COVID-19 BR and CoMo Consortium for the collaborative work. The authors also thank the research funding agencies: São Paulo Research Foundation (FAPESP) -Brazil (grant number: 305