Containing pandemics through targeted testing of households

Background While invasive social distancing measures have proven efficient to control the spread of pandemics failing wide-scale deployment of vaccines, they carry vast societal costs. The development of a diagnostic methodology for identifying COVID-19 infection through simple testing was a reality only a few weeks after the novel virus was officially announced. Thus, we were interested in exploring the ability of regular testing of non-symptomatic people to reduce cases and thereby offer a non-pharmaceutical tool for controlling the spread of a pandemic. Methods We developed a data-driven individual-based epidemiological network model in order to investigate epidemic countermeasures. This models is based on high-resolution demographic data for each municipality in Norway, and each person in the model is subject to Susceptible-Exposed-Infectious-Recovered (SEIR) dynamics. The model was calibrated against hospitalization data in Oslo, Norway, a city with a population of 700k which we have used as the simulations focus. Results Finding that large households function as hubs for the propagation of COVID-19, we assess the intervention efficiency of targeted pooled household testing (TPHT) repeatedly. For an outbreak with reproductive number R=1.4, we find that weekly TPHT of the 25% largest households brings R below unity. For the case of R=1.2, our results suggest that TPHT with the largest 25% of households every three days in an urban area is as effective as a lockdown in curbing the outbreak. Our investigations of different disease parameters suggest that these results are markedly improved for disease variants that more easily infect young people, and when compliance with self-isolation rules is less than perfect among suspected symptomatic cases. These results are quite robust to changes in the testing frequency, city size, and the household-size distribution. Our results are robust even with only 50% of households willing to participate in TPHT, provided the total number of tests stay unchanged. Conclusions Pooled and targeted household testing appears to be a powerful non-pharmaceutical alternative to more invasive social-distancing and lock-down measures as a localized early response to contain epidemic outbreaks. Supplementary Information The online version contains supplementary material available at (10.1186/s12879-021-06256-8).


Description of individual-based modeling (IBM) software
When selecting and designing our computational modeling framework, we have explicitly taken into consideration the following aspects: (A) It should be possible to directly test a wide variety of interventions, such as partial workplace and school closings, social distancing, and testing with subsequent isolation of infected individuals. (B) Each infected person should have a disease-state development where the stages have realistic time delays. If a lag in response time is not included, we anticipate that all estimated response times to interventions will be markedly wrong. This is especially evident for the correct estimation of R, and thus, the capacity to determine how quickly R responds to interventions. The mentioned challenges are inherent to most metapopulation-based models, whereas an IBM approach is not hampered with such issues.
Guided by the above reasoning, we have developed a complex system modelling framework to describe the spread pattern of COVID-19 in mainland Norway. The model design and build features are chosen with the intent that it can be used to carefully assess a multitude of relevant intervention strategies. Briefly, the model is based on using an IBM based on complex network theory for each Norwegian municipality. The features of the model's layers and each individual all adhere to highresolution demographic data specific for the respective municipality. The disease-state of each individual follows an epidemiological SEIR-type model.
Our total model for all combined municipalities contains approximately 5.3 million people and simulates dynamics of spread through social interactions in layers such as households, schools and workplaces.
network are constrained such that known high-resolution demographic data for each municipality is matched: 1) The household layer consists of separate households with a size and age distributions that follow demographic data for that municipality.
2) The number of schools, type, and their student populations are based on demographic data.
3) The number of day care facilities is based on demographic data.
4) The number of nursing homes and their population sizes are based on demographic data. 5) If a household has multiple children in e.g. day-care age, these children are assigned the same day-care. Similar for primary and secondary schools. For high-school age, the assignment to school unit is random. 6) In the work layer, the number of companies and their sizes are based on demographic data.
Note that this layer is intended to represent spread between co-workers. For professions with large amounts of exposure to the general public, contacts with customers are represented through the generic contact network.
7) The generic contact layer is designed to capture heterogeneity in a person's daily contact patterns. While the other layers consist of fixed cliques, in the random layer, each node connects to a random selection of other nodes for each simulated day. For each connection between an infected and a susceptible node, there is a fixed probability of the susceptible node becoming infected. The daily number of connections for each node is redrawn each day, according to a uniform distribution ranging from 1 to a maximum number of daily contacts C D , which is a node attribute assigned at the initiation of the simulation. 8) We include two different modes: (A) For young (< 20) and elderly (80+) we assign a maximum number of daily contacts CD following a normal distribution, whereas the remaining age groups are assigned contacts following a combination of a normal and a power law distribution: Contact number is drawn from (i) normal distribution, (ii) power law before being added together. The latter generates larger contact heterogeneity. 9) Using data from Statistics Norway 1 , adults with different work and household municipalities are assigned work locations in the correct municipality. Supp. Table 2 lists the data tables   used. 10) When a person is committed to a hospital, they are removed from their domicile (household or nursing home). 11) When an infected person manifests a symptom, is confirmed COVID-19 positive, or an asymptomatic person is confirmed COVID-19 positive, the model assumes they will selfquarantine from activities in all layers except their domicile.
All individuals participate in the generic contact network, which is designed as a random timedependent scale-free network to capture heterogeneity in contact patterns. A new instance of this network (per municipality) is generated new every day. Except for the generic contact network, each grouping in each layer is represented as a k-clique, where every individual is connected to all others in the group, thus representing a well-mixed group. Supp. Figure 1 shows a schematic of the IBM layers inside which there are smaller groups (left panel) and a simplified example of a possible resulting contact network (right panel). The probability of infection depends on which layer the infected contact is located. The specific values used for the parameters in each layer are detailed in Supp. Table 1. Note that the infection probabilities apply between all nodes in a given layer, and as such represent a combined probability of contact and infection. While we may assume that some layers would see an actual close contact between all members of a given clique for each day (for instance, households), other layers, such as schools and workplaces, are not likely to actually see this form of all-to-all contact on any given day. This is reflected in a substantially lower combined contact and infection probability for these layers.

Population size scaling
We developed a schema for scaling up the population size of a municipality that preserves its demographic distributions. The scale parameter  determines the new size of a municipality:  =1 for the unscaled population,  = 0.5 results in a municipality with half of the population, and  =2 gives a municipality with twice the original population. As the statistical data, age data and household-, nursing home-and workplace-compositions obtained from Statistics Norway take the form of percentages, the distributions are independent of the size of the network and can be used as provided, simply by multiplying by  the number of desired individuals/households/workplaces (provided by Statistics Norway and passed directly as inputs to the generating code when  = ). On the other hand, schools and daycares are explicitly and individually defined, and their numbers are too few to be mapped to fine-grained size distribution. Therefore, in order to scale these systems to match the scaling of the whole municipality, we instead opt for an approach where we randomly select and duplicate (or erase, for   ) schools and daycares in the original municipality data until we reach the number specified by the scale factor.

Variation of household-size distribution
In order to evaluate the impact of the household size distribution on spreading potential, we needed to establish a systematic way to modify the distribution, preferably defined by a single scalable parameter α. We decided on either splitting or merging existing a fraction of households in the original distribution, yielding smaller or larger households respectively. In the first case, we identify households eligible for splitting, using the criteria that it must contain at least two adults. This is necessary in order to ensure that we do not create households containing only children. We then define α as the fraction of eligible households to split. The split itself is conducted by assigning the two first adults in the selected household to their own household, with each of the remaining individuals assigned to either household by coin flip (Bernoulli, 0.5).
We represent the merging of households by using negative values of α: As merging households with at least one adult each will naturally lead to a new household with at least one adult, all households are eligible for merging. We therefore randomly assign the households into pairs, with a chance |α| that determines the fractions of pairs which are combined to form a single household consisting of the members of the two original households, reducing the total number of households in the system by one for each merged pair. The households in the remaining pairs stay unchanged.
As an example, consider a fictitious small municipality containing 2000 individuals, labeled i1, i2, i3, … , distributed across 1000 households H1, H2, H3, … such that H1 contains individuals i1, i2 (using . We repeat this process for X2, adding H3 and H4 to N either separately (with probability 1-|α|) or as one combined household [i4, i5, i6, i7, i8, i9] (probability |α|), and similarly for each pair in X. If the number of households is odd, one household is copied directly into N without possibility of pairing. The remaining households are assigned to pairs as described.
Note that N will have a lower total number of households than H, but the same total number of individuals. For instance, setting α = 0.1 in a network starting with 1000 households, the average result is N consisting of 950 households (900 households in 450 pairs remain unchanged, with the remaining 100 households merged pairwise into 50 new ones). For α = 0.5 in a network starting with 1000 households, we would get (on average) 750 households, with 500 households unchanged and the other 500 being combined pairwise in order to form 250 households.  Pre-symptomatic and asymptomatic individuals are infectious and active in all layers they are associated with (for most, this will be the generic contact layer, a household, and either a workplace or educational/daycare institution). Symptomatic individuals are assumed to be infectious in the household and nursing home layers only. Hospital and ICU protective measures are assumed to be sufficient to prevent spread from patients to staff. The direct transition from Is to D captures the disease trajectory for some people in nursing homes.

SEIR-type epidemiological dynamics
The transition from S to E is governed by disease transmission probabilities in the IBM contact network. After an individual is infected from a neighbor in the contact network, the SEIR-type dynamics of that individual's disease state (onward from E) are governed by stochastic processes with appropriate waiting times according to empirical data for COVID-19. We continually update our parameter estimations given the daily update in disease state numbers for Norway by region.   Table 1. Once the next state is determined, we generate the duration of the state according to a Poisson-distributed random variable (with an appropriate λ) plus 1 (in order to avoid two state changes in one day) and set the date of next change of state accordingly.
Note that our SEIR-model uses age-stratified transition rates, where the age of an individual decides which of the age groups that individual is part of Supp. Table 1 shows the rates used for the different SEIR transitions of Supp. Figure 2.

COVID-19 testing of (sub-) population
The testing itself is represented by drawing individuals from the symptom-free population (susceptible, latent, asymptomatic, pre-symptomatic, or recovered), and returning a positive test if they are asymptomatic or pre-symptomatic. Pooled tests are represented similarly, but with multiple people (forming a pool) simultaneously tested, and all tests marked as positive if any of the people in the test pool are in the asymptomatic or pre-symptomatic states. Our IBM implementation allows for random selection of individuals for testing, or testing based on a pre-conceived schema for node attributes.

Base-line fitting of IBM to empirical data
The parameters in the stochastic IBM are treated as fixed, and they are either set from literature or estimated directly from publicly available empirical data, as indicated in Supp. Table 1. We determined model parameters by fitting the predicted hospitalization rate of our Oslo model to hospitalizations for Oslo in the period from March 1 st to April 20 th . We assume a sudden transition between two regimes, with infections following original "unrestricted" probabilities until March 13 th .
The "unrestricted" regime is initialized with 20 cases at an unspecified date and run until the system reaches Nlock= 400 active cases. This point is marked as March 13 th and serves as a reference point for the remainder of the simulation. The motivation for fixing dates at this point (rather than the initial point of 20 cases) is two-fold: First, it is important that the model follows a progression of nonpharmaceutical interventions (NPIs) that actually corresponds to those implemented by Norwegian authorities. Therefore, we decided to use the earliest NPIs introduced (March 13 th ) as an anchoring point, with the dates of subsequent NPIs defined in relation to the first. Second, due to the heterogeneity in the infectious potential of individuals (with a small fraction of individuals being substantially more infectious than others due to their particular contact network), the model exhibits highly stochastic behavior when the number of infections is low, as the number of highinfectivity individuals infected is more variable. As small fluctuations in the early evolution of the epidemic can have dramatic effects on long-term forecasts, anchoring the network around a reasonably large number of cases allows the model to make more consistent predictions. This approach also solves the issue of highly stochastic behavior at low infection levels, as 400 active cases is sufficient to reduce this notably.
After these steps, all schools and day care facilities close, as do 50% of workplaces. Random contacts are reduced by 83%. Infection probabilities in nursing homes and within-households remain the same in both regimes. Model simulations start in March and continue until August 30 th , 2020.
Model fitting to the data was then done by adjusting infection probabilities in order to satisfy the following criteria: (1) matching slopes (both increasing and decreasing) and (2) the location (date) and duration of the peak of the hospitalization data (Fig 2A, red curve) Fitting was done by adjusting infection probabilities alone until a satisfactory fit (by visual judgment) was reached. While the epidemic curve is also affected by state durations, these were kept fixed based on reliable direct clinical data which was readily available (as substantial deviations from these would in any case be hard to justify).
Since the lockdown in Norway occurred in the early stages of the epidemic, the fitting of the IBM results to data depends very much on the number of symptomatic infected on the day of lock down (March 13, 2020). Thus, for a given parameter set and municipality, we have used a large number of simulations to determine the optimal number of infected (Nlock) on this day for the model predictions to fit to time series data of number of hospitalized persons. Simulations are subsequently conducted by initiating the network with a small number of infections and running the code until Nlock is reached. As described above, the date is set to be March 13 at this point. We anticipate that several of the IBM parameters can be estimated from other data sources when they become available, e.g. large-scale systematic testing regimes for COVID-19.

Description of the TPHT process Selection of households
The central aspect of TPHT is the prioritizing households for testing according to size. As an example, consider a municipality with the following distribution of households: • 6 persons and over: 3,000 households • Total: 200,000 households, population of above 420,000 (depending on the exact distribution of households of size 5 and more) In this case, a TPHT test fraction of 2.5% corresponds to 5,000 households, meaning we select all households 6 and over (3,000), and a random choice of 2,000 of the 5-person households. A TPHT test fraction of 5% would yield 10,000 households, meaning all households 5 and above. A TPHT test fraction of 10% gives 20,000 households, meaning all 10,000 households five and above, as well as 10,000 randomly chosen 4-person households. For a given fraction and test intervals, the selected households are distributed evenly across the interval. For instance, using a test fraction of 3.5% on the municipality above with a test interval of 7 days, we would have to test 7,000 households a week; with 1,000 periodically tested every Monday, another 1,000 tested every Tuesday, and so on for the duration of the TPHT test regime.

Estimation of R for a choice of simulation parameters
Due to the individual-based nature of our model, each infection can be traced to the individual causing it. This also makes it possible to keep count of the number of secondary infections caused by each infected individual, which should make computing R fairly straightforward. However, there are some caveats. We want to determine R for a given point in time, but a closer consideration of the definition of the reproduction number reminds us that it is the number of secondary cases generated by an original case, and as such the strict definition of R applies to an infected population rather than a point in time.
We can, however, define a daily reproduction number by taking the average number of secondary cases caused by each individual infected on a given day. This number cannot be computed in real time, as any individual that does not recover on the specified day may very well cause additional infections before it recovers. Therefore, we opt to determine R through a "two-pass" process. We begin by simulating the epidemic over a given time We make three key observations: first, there is some stochasticity in the determined value of R between consecutive days; second, the transition from one regime to the next exhibits an intermediary phase of a few days until R begins to stabilize around a new average. We propose two primary explanations for this. Firstly, as we compute R by taking the average number of infections caused by all individuals at a chosen time, individuals that were already infected before measures were introduced will have been more infectious during that period, driving up the average. Second, different layers are unevenly affected by measures. Consider a transition from free spread to a regime where only workplace spread is possible, but school, random, nursing home and even household spread is entirely prevented. As restrictions enter into force, a few workplaces will have active cases. These will cause further infections for some time, either gradually petering off or saturating each environment, in which case R would abruptly drop to zero.
Lastly, we see a more abrupt drop-off in the estimated reproductive number as we reach the last few days of the simulation period. This is because many of the last nodes to be infected do not have time to recover, and with the requirement that a node must recover in order to have its reproductive numbe included in the average, there is a bias towards individuals with shorter illnesses, which pulls down the average R.
Meanwhile, estimating R over longer periods of time leads to substantial discrepancies: if R > 1, there is a build-up of immunity, leading to a gradual reduction in R; if R < 1, infection levels will gradually drop off, leading to situations with very few infected and correspondingly high noise in the relative number of new daily infections. Because of this, for purposes of determining R, our process is as follows: • First, initialize the model with 20 infected individuals • Unfettered growth until 100 symptomatic cases are reached (this is enough to remove most of the noise, while not enough to build up substantial immunity in the population). The purpose of this approach (rather than just beginning with a 100 infected individuals on the first day) is described in the next section.
• Once 100 (simultaneous) symptomatic cases are reached, transition into the regime for which we wish to determine R, and run this for 40 days. Due to the high infectivity up to this point, there is also a substantial number of pre-symptomatic cases already infected at this point.
• Take the daily average R from 17 to 22 days after the transition mentioned above (equivalent to the 27-day interval in Supp.

Simulation of outbreaks
When conducting the simulation of an outbreak consisting of Nout=1,000 symptomatic infected, the simple approach is to randomly select Nout nodes in the network and abruptly change their state from susceptible to exposed. However, this will cause artificial oscillations in the data since all of the exposed individuals start out in synchrony because they were infected at the same time. It takes several infection generations for these oscillations to subside. Instead, we dynamically generate Nout by starting from a small seed of 20 randomly (simultaneously) selected and exposed individuals.
Subsequently, the IBM is simulated until we reach Nout symptomatic individuals. This approach for achieving Nout does not cause infections to be focused on larger households, instead spreading them randomly and asynchronously in the network due to the large fraction of infections taking place in the random contact layer. The state with N out symptomatic individuals is used as the starting point (day 0) of the outbreak simulation, from which interventions may be implemented. We go through this initiation process independently for each separate simulation, thus ensuring stochastic variation in the day 0 state.

Software availability
The IBM software used for performing the simulations is available on two general access platforms: R as a function of test fraction with 50%, 70%, 90% or 100% of households willing to participate in TPHT. The testing frequency is set to once every 7 days for each of the four curves. The x-axis corresponds to the number of households actually tested (and therefore actual burden on testing capacity). For instance, the number of tests actually administered for a test fraction of 0.1 is the same for 50% and 100%; however, 50% compliance means that half of the 20% of largest households are tested. In contrast, 100% compliance means that all of the 10% largest households are tested weekly. Parameter values in black are used in computer simulations for Fig. 2 in the main article. To generate Supp. Figure 4, parameter values in red were used instead when present in the table. The infection probabilities given under the section "Individual-based network model" correspond to daily secondary attack rates in a clique with one symptomatic individual (except for the generic layer, which is not clique-based). For example, in a household with one symptomatic infected individual, each other individual has a 5% chance of being infected per day. The infectivity of pre-symptomatic and asymptomatic individuals for a given type of interaction are adjusted by a factor 3 and 0.5, respectively. Children are half as likely to contract the disease from a given interaction. In cliques with multiple infected individuals, the probability increases multiplicatively. For example, consider a household with two symptomatic, one asymptomatic and one pre-symptomatic case (all adult).
Here, the daily infection probability for remaining household members becomes: The specified infectiousness of generic contacts is the base probability of transmission from a single random contact.