Prioritizing allocation of COVID-19 vaccines based on social contacts increases vaccination effectiveness

A.1 Generating synthetic populations and networks

A synthetic population of a region may be regarded as a digital twin of the real population of that region. Here we provide a compact summary of the model and the methodology behind constructing synthetic populations and their contact networks in the case of the US; see [47] for details. Our work builds on earlier techniques for a first principles approach for constructing synthetic populations [7, 23, 24]. These populations and networks are central to the EpiHiper simulation model.

To construct a population for a geographic region R (e.g., Virginia), we first choose a collection of person attributes from a set 𝒟 (e.g., age, gender, and employment status) and a set 𝒯_A of activity types (e.g., Home, Work, Shopping, Other, and School). The precise choices of 𝒟 and 𝒯_A are guided by the particular scenarios or analyses the population will serve. Described at a high level, we (i) construct people and places, (ii) assign activity sequences to people, (iii) map each activity for each person to a location (including the time of the visit), and (iv) from this, we derive a contact network using co-occupancy to infer edges. The construction is broken down in a sequence of steps outlined as follows.

Using iterative proportional fitting (IPF) [8, 20] the base population model constructs a set of individuals 𝒫 where each person has assigned demographic attributes from 𝒟. By design, this ensures that 𝒫 matches the actual distributions and Public Use Microdata Sample (PUMS) data from the US Census [64], which is the input data for the model. Additionally, this model partitions 𝒫 into a set ℋ of households, where the notion of household encompasses the traditional notion of “family”, but also any other subset of individuals residing in the same dwelling unit (e.g., dormitories, army barracks, or prisons).

After household assignment, each individual p ∈ 𝒫 is assigned a week-long activity sequence α(𝒫) = (a_i,p)_i where each activity a_i,p has a start time, a duration, and an activity type from 𝒜. Data sources used for this step include National Household Travel Survey (NHTS) [66], American Time Use Survey (ATUS) [65] and Multinational Time Use Study (MTUS) [61]; these sources are fused to form consistent, week-long activity sequences. We write for the mapping assigned to each person. For this construction, we use Fitted Values Matching (FVM) for adults [39], and Classification And Regression Tree (CART) for children (see, e.g., [11]).

The location model constructs a set of spatially embedded locations ℓ consisting of residence locations where households live, and activity locations where people conduct their non-Home activities. This construction is highly granular and is rooted in data such as the MS Building data [44], HERE/NAVTEQ data [33] for points-of-interest (POIs) and land-use classifications, National Center for Education Statistics (NCES) [49] data for public schools, as well as LandScan³, OpenStreetMap⁴, and Gridded Populations of the World (GPW) v4⁵. A point plot of the locations of Virginia is shown in Figure 12.

For each person p ∈ P, the location assignment model assigns a location ℓ = ℓ(a_i) to each of their activities a_i. We denote the sequence of locations visited by 𝒫 as λ_p = (λ_i)_p. The location assignment model uses American Community Survey (ACS) commute flow data [63] to assign the target county c for Work activities, and a particular location randomly within c work weights assigned to each location in c. School activity locations are assigned based on NCES data, with remaining activities anchored near home and work locations.

Finally, the contact network model uses the location assignment to derive the bipartite people location graph G_{P L} with vertex sets V₁ = 𝒫 and V₂ = ℓ and a labeled edge (p, ℓ) whenever 𝒫 visits ℓ where the label includes activity type, time for start of visit, and duration of visit. From this, we derive the list of visitors to each location and the co-location graph G_max with vertex set 𝒫 and edges all e = (p, p^’) for people p and p^’ that are simultaneously present at the same location. Merely being present at a location at the same time does not imply a contact, and sub-location contact modeling is applied at each location to determine which of the edges of G_max should be retained to form the contact network G which is also referred to as the person-person contact network and denoted as by G_{P P} (rather than simply G) to make this explicit. In this work, we use a random graph model referred to as the Min/Max/alpha model at each location to obtain G. Let ℓ be a location and let N = N_l denote the maximal number of simultaneous visits to ℓ. Define the function where Min < Max are non-negative numbers and α > 0. Given p = p_ℓ(N) one samples from this random graph model in the same manner as for the standard model G_n,p by independently at random applying to each edge e the probability p corresponding to the location ℓ where e ∈ G_max originate. Thus the parameters Min and Max bound the degree of each vertex locally at ℓ (in expectation) for each visit; note, however, that the degree of person p in the resulting graph G is the accumulation of degrees across their trajectory to locations visited while executing their activity sequence. Thus the choices of Min, Max and α will induce the degree of each vertex in a bottom-up manner, see [47] for full details. Finally, for the applications and scenarios of this paper, we project from G, the week-long contact network, to G_Wednesday, representing the contact network on a “typical day”.

We will use agent-based network models to study the epidemic process in this paper. Agent-based networked models (sometimes just called agent-based models) extend metapopulation models by explicitly capturing the interaction structure of the underlying populations. In this class of models, epidemic dynamics are modeled as a diffusive process on a specific undirected person-person contact network G(V, E) = G_{P P} (V, E) on a population V (see Section A.1 for precise definitions) – the existence of an edge e = (u, v) ∈ E implies that individuals (also referred to as nodes) u, v ∈ V come into contact⁶. Let N (v) denote the set of neighbors of v.

A.2 Interventions and vaccine allocation policies

Interventions are implemented to inhibit disease transmission. Interventions can be thought of as individual behavioral adaptations or policy mandated changes, such as closing certain facilities or reducing their capacity. Of course, policies also lead to further behavioral adaptations. Our agent-based models have a rich set of interventions implemented. The specific ones we use in the study are detailed in Section F. In our simulations, all policy changes and behavioral adaptations except vaccine uptake can be seen as processes that continually and adaptively change the social contact network and disease transmission parameters.

Here we focus on policy concerning vaccine allocation. Given a schedule S that specifies the amount of each vaccine available at each time, and the characteristics of the vaccines, a prioritization scheme is a policy that assigns at each time period the individuals (nodes) that are to be vaccinated. In other words, a vaccine prioritization scheme can be thought of as a Markov Decision process – at each time step we know the current state of the system and the available vaccines, and we need to decide who gets the vaccine. In this paper, we only consider non-adaptive policies – i.e. policies that do not change how vaccines are allocated at each time step. The amount depends, of course, on the schedule. We focus on four types of policies, each of which partitions the population into priority groups using different characteristics:

• 1. age group,

• 2. occupation,

• 3. number of social contacts (degree), and

• 4. total duration of social contacts (weighted degree).

• Within each type, the policies are distinguished by the fraction allocated to each priority group (further details can be found in Section 2.2).

Diameter

The diameter is the length of the longest shortest path between any two vertices of the network. The diameter of G is diam(G) = 12.

Degree distribution

The degree of a person u in the people-people contact network G is the number of different persons that u has contact with during a day. The average person degree in G is = 43.5. In the people-location network G_{P L}, the degree of a location ℓ is the number of distinct visitors to ℓ during a day, and the degree of a person u is the number of distinct locations u visits during a day. Figure 14 shows the degree distributions in G and G_{P L}. The weighted degree of a node u is the total time of contact with its neighbors per day. The weighted degree distribution of the network G is shown in Figure 15a.

• Download figure
• Open in new tab

Figure 14:

Degree distributions in people-people network G: (a) as a histogram in normal scale, and (b) in log-log scale showing a power-law tail with an exponent of about −7; and degree distribution in people-location network G_PL: (c) of locations in log-log scale showing a power-law tail with an exponent of about −1.7, and (d) of people as a histogram in normal scale.

• Download figure
• Open in new tab

Figure 15:

Distributions for connectivity properties of G here including (a) the weighted degree distribution, (b) the clustering coefficient distribution, and (c) the normalized core number distribution for G. of the people-people network.

Clustering coefficient

This is a measure of the degree to which nodes in a graph tend to locally cluster together. For a node v, its local clustering coefficient is the fraction of pairs of its neighbors which have a link between them. It quantifies how close the immediate neighbors are to being a complete graph. If the local clustering coefficient is 1, it means that the node and its neighbors induce a complete graph, and if it is 0, then they induce a star graph with v at the center. The average clustering coefficient is the average of all the local clustering coefficients [70]. The average clustering coefficient of G is = 0.092; the clustering coefficient distribution is shown in Figure 15b. k-core. A k-core of a graph G is a maximal connected subgraph of G in which all vertices have degree of at least k. Equivalently, it is one of the connected components of the subgraph of G formed by repeatedly deleting all vertices of degree less than k [71]. The core number of a graph is the maximum core it belongs to. The maximum core (max. core) of the graph is the set of all nodes with the maximum core number. It is another way of characterizing centrality of nodes. It is considered that for certain epidemiological models, the k-core number is a good predictor of the final outbreak size [50]. Also, the k-core decomposition is in many cases a good predictor of spreading efficiency. The k-core decomposition of G is shown in Figure 15c. The maximum k for which there is a k-core is 57. While there are only 1,700 nodes in the induced subgraph for the 57-core, there are approximately 35,000 nodes (around 0.5% of the total population) in the 45-core, and more than a quarter of the population belongs to the 30-core showing high global connectivity

Graph spectrum or eigenvalues

The spectrum of a graph is the set of eigenvalues of its adjacency matrix. Also popular is the Laplacian spectrum, which is the set of eigenvalues of the Laplacian matrix of the graph. There are several works that relate spectrum, particularly the first eigenvalue of the adjacency matrix, to disease spread in SEIR-like models [27, 52]. The common result that highlights the impact of the network structure on the dynamics is that epidemics die out “quickly” if λ₁(G) ≤ T, where λ₁(G) is the spectral radius (or the largest absolute value of an eigenvalue) of graph G, and T is a threshold that depends on the disease model. This relationship has motivated a number of works on epidemic control where the objective is to find an optimal set of nodes (or edges) to remove from the network that leads to maximum reduction in its spectral radius [57, 67, 74].

Activity-based structural analysis

We analyzed the constructed social contact network using the structural measures described above. We considered activity induced sub-networks, where a contact edge is retained only if both individuals corresponding to that edge have the target activity assigned to them. For example, in the case of School network, only School–School edges are retained. For School activity, we observed an average degree of 29.7, which is between the average degrees of 13.5 and 47.3 reported for the school networks of the SocioPatterns collaboration networks [14, 28]. Our average degree for Work activity is 16.3. In comparison to the SocioPatterns data, this seems to be on the higher side for a typical office environment (< 7) but comparable to their hospital network (14.0). This is expected, since our Work activity includes a wide variety of workspaces that include office environment, factories, restaurants, etc. The max. core for the School activity network is 32. In SocioPatterns, the two schools have max. cores of 24 and 47 respectively. For the Work network, it is 37, while in SocioPatterns it is 11, 25 and 23 respectively for the three office networks. Again, in our case the max. core is the maximum among all locations.

C.1 Epidemic curves with various strategies

Figures 16 to 18 show the cumulative numbers of infections, hospitalizations, and mortality in the first three months of 2021 under different vaccination schemes, assuming no relaxation of NPIs (Figure 16), slow relaxation (Figure 17), and fast relaxation (Figure 18). Note that for the age-based, essential worker-based, and degree-based schemes, we only show the moderate level of priority enforcement 80%. All curves show the uncertainty of one standard deviation above and below the mean.

• Download figure
• Open in new tab

Figure 16:

Cumulative counts under different vaccination strategies for (a) total infections, (b) total hospitalizations, and (c) total mortality, assuming current non-pharmaceutical interventions remain at the same level (as-is).

• Download figure
• Open in new tab

Figure 17:

Cumulative counts under different vaccination strategies for (a) total infections, (b) total hospitalizations, and (c) total mortality, assuming slow relaxation of non-pharmaceutical interventions.

• Download figure
• Open in new tab

Figure 18:

Cumulative counts under different vaccination strategies for (a) total infections, (b) total hospitalizations, and (c) total mortality, assuming fast relaxation of non-pharmaceutical interventions.

C.2 Sensitivity analysis

This section shows some of the variability in network measures for a statistical design over the network construction parameters Min, Max and α. Specifically, we consider the hyperplane given by Min = 5, α = 1000 and Max ∈ {35, 40, 45} with Max = 40 representing the base case used in all the simulations of this work. Accordingly, we have three graphs G-35, G-40, and G-45. As can be seen in Table 2, there is a slight increase in properties such as , , and d_max while the diameter drops by 1 for Max = 45. The slightly bigger change occurs in eigenvalues, but we note that the spectral gap (i.e., λ₁ − λ₂) remains nearly constant.

Why do degree-based heuristics work

The efficacy of degree-based heuristics has been discussed in several earlier papers, e.g., [3,10,22,51,73], as discussed in Section H. This phenomenon is also well understood in certain random graph models [10, 51]. Pastor-Satorras et al. [51] derive an immunization threshold in terms of the fraction of high degree nodes that are immunized through an analysis that relies on a degree-based mean field (DBMF) assumption. This is derived more rigorously by Bollobas et al. [10], who show that there is a threshold fraction of high degree nodes whose removal shatters the graph—this can be viewed as an immunization strategy in their random graph model for the case when the transmission probability 𝒫 = 1.

However, it is important to note that the theoretical results are not directly applicable in our context for the following reasons: (i) many of the theoretical results show the efficacy of these methods for power law networks – the networks we generate are similar to power law networks but with a very different exponent; and (ii) many of the results are shown when vaccines are applied at the start of the epidemic process, and the results do not say anything of what happens when the vaccine is applied temporally – this is important, because the temporal epidemic process infects individuals already and thus changes the network structure substantially, including the application of the non-pharmaceutical interventions.

As nodes are infected or vaccinated, they are either not susceptible or less susceptible (i.e. lower probability of getting infected). Consider the network structural changes due to such dynamics. For simplicity, we remove all edges incident on nodes who are or have been infected and those who are vaccinated, to form a sparser snapshot network at different time point, especially right before and right after each mass vaccination. In Figure 21, we compare the degree distributions of such snapshot networks from age-based and degree-based vaccinations, before and after each vaccination. The networks are extracted from a simulation run under as-is scenario with fast vaccine distribution and 100% prioritization. Before the first vaccination, the snapshot network from both strategies have the same degree distribution. After the vaccination, since more edges are removed with the degree-based strategy, the snapshot network has more low degree nodes and fewer high degree nodes, compared to those from the age-based strategy. Similar is observed between the snapshot network before and after the second vaccination. After that the degree distribution does not change much with vaccination, since the network is already very sparse.

• Download figure
• Open in new tab

Figure 21:

As nodes are infected or vaccinated, we remove edges incident on them from the contact network, to form sparser and sparser snapshot networks. We consider such snapshots under age-based and degree-based strategies, assuming no relaxation of NPIs and 100% distribution. After the first round of vaccination, the degree-based strategy reduces degrees more than the age-based strategy: we observe more low degree nodes and fewer high degree nodes in the snapshot networks from the degree-based strategy.

The connection between spectral properties and epidemic thresholds provides another insight into the efficacy of a degree-based strategy. As mentioned in Section B, it has been shown for SIS/SIR models, via different approaches (including different mean field approximations), that if the first eigenvalue (referred to as the spectral radius) λ₁(G) of the network is below a certain threshold, the epidemic is not very “large” [27, 45, 52]. Figure 19a shows the first eigenvalue λ₁(G) of the network at different points in time as interventions are applied. It is very significant to note that λ₁(G) shows a substantial drop when the degree intervention is applied. Thus, the intuition behind the efficacy of such methods is clear – high degree nodes protect not only themselves, but also confer a higher level of indirect protection on their neighbors as they interact with many individuals. The fact that real-world social networks have sufficient nodes of high degree ensure that such heuristics are effective. It is important to note, however, that just the presence of high degree nodes themselves does not ensure that degree-based heuristics work.

A dynamical property: Vulnerability and its relationship to degree

Given a generalized SEIR system G(V, E) along with a random initial configuration and a node v ∈ V; and an integer t ≥ 1, we define t-Vul(v) as the probability that v gets infected at time t. Vulnerability of v by time t (t-TotVul(v)) is the probability that v is gets infected by time t. Informally speaking, t-Vul(v)and t-TotVul(v)estimate the chance of a node getting infected when a random set of nodes are infected to begin with. Vulnerability in essence captures various ways the infection can reach a given node. Degree-based interventions effectively reduce the vulnerability of several nodes at once. Interestingly, nodes of high degree are also highly vulnerable, and this can be empirically observed by noting the plots in Figures 22 and 23. Vulnerability computations are computationally challenging (#P-hard) and our simulations provide good estimates of the same.

• Download figure
• Open in new tab

Figure 22:

Degree distribution among nodes that are prioritized differently. Essential workers and older people seem to have lower degrees than high degree people and highly vulnerable people.

• Download figure
• Open in new tab

Figure 23:

Vulnerability distribution of nodes with different priorities. People of highest degrees have high vulnerability, close to that of the most vulnerable people.

How can we estimate degrees

A central question is: can individual node degrees be estimated? We believe that this is possible, and more so as a result of digital devices and apps that have recently been deployed for contact tracing. The current set of apps built for contact tracing keep track of the social interactions of an individual in a privacy-preserving way. These apps can easily be modified to count the number of such interactions. Furthermore, our results show that duplicate counts that might result are okay, as the degrees are likely to be skewed some. These apps can be modified to not only store the keys, but also the encrypted form of the time a person spends with the person that generates such a key.

Furthermore, a person does not need to reveal their degree to anyone but the app can simply provide this information for prioritization. The alerting and scheduling of vaccination appointments can be conducted through the app itself. Additionally, demographic information such as age, or other risk-factor information such as additional comorbidities or race and/or ethnicity, may also be used to further prioritize individuals within their degree-based allocation groups.

How accurate can such measurements be

Let V be partitioned into groups V₁, . . . , V_r, and suppose p_i is the probability that a node in V_i has the app. Think of V_i as an age group. Survey results show that younger individuals have a higher propensity to use digital device as well as download the app. We assume each node v ∈ V_i decides independently with probability p_i to install the app; if v chooses to install the app, we say that it is “sampled”. Let H denote the sampled graph, which is induced by the nodes (node induced subgraph) which have the app, and let V (H) denote the set of nodes in H. For a node v, let N_G(v) denote the set of its neighbors in G. Similarly, let N_H (v) = {u: u ∈ N_G(v) and u ∈ V (H)} denote the set of neighbors of v in H. Let d_G(v) = |N_G(v)| denote its degree in G, and let d_G(v, i) denote the number of neighbors it has in group . If v is sampled, let d_H (v) = |N_H (v)| denote its degree in H. Let X(u) ∈ {0, 1} be a random variable which is 1 with probability p_i, if u ∈ V_i, and let . Observe that . By applying a Chernoff bound to Y (v), we have Therefore, with probability at least 1 − 2/n, every node . Let Q ⊆ V denote the subset of nodes in the top quartile with respect to d_G(v). Then, with high probability, for every node v ∈ V_h ∩ Q, we have d_H (v) ∼ f (v). Therefore, if the sets V_h and Q have high overlap, high degree in H corresponds to high degree in G. The expected number of nodes in V_h ∩ Q which get sampled is

Where is the tradeoff?

The phased approach is a balance between exposure/infection risk (policy equivalent of high social contact) and disease risk (age-based). Our proposal is a stronger version of the former, but seems to be significantly better across all outcomes. This is one of the counter-intuitive aspects of the proposed approach, since it does not take into account disease risk explicitly, and purely relies on structural measures of the network. While individualized policies aimed at minimizing mortality do exist (e.g., comorbidities) and are part of the phased approach, lack of technology thus far had made it difficult to ‘individualize’ policies targeted on minimizing transmission. Aggregate policies targeting essential workers and other high exposure categories (e.g., school children for influenza) while in a similar vein to our approach, may still have high variance in reducing transmission potential. The advantage of our approach is in leveraging the mechanistic and network based understanding of disease spread, and creating priority categories that cut across age-, risk- and other demographic characteristics. It is possible to construct pathological instances of networks in which vaccinating high degree nodes is not optimal. But as is well known in the field of network science, real world social networks are often characterized by power-law degree distributions and exhibit a hub-spoke structure. In such networks, targeting the hub nodes at individual level, ‘shatters’ the transmission network faster than aggregate policies. Especially such a strategy might be necessary in regions with very limited vaccine supplies and at higher risk of variant induced surges. We believe our work helps advance the case for better integration of social network and digital technologies in swifter public health response.

Disease Progression

covers the health state transitions within an individual 𝒫 that are independent of other people. For the EpiHiper model, a disease progression diagram describes all the possible health state transitions that take place within a person. The nodes of the diagram are the health states X = {X_i} and directed edges of the form e = (X_i, X_j) with an assigned probability p_e = prob(X_i, X_j) and a dwell time distribution D_e. The sum of all probabilities of transitioning out of a given state X_i must be either 1 or 0. Zero indicates a terminal state.

The System State at any point in time is given by the attributes of the individuals (nodes) and contacts (edges) of the network, the simulation time, and the value of user-defined variables.

The values nodeTrait and edgeTrait are user-defined attributes which may be used to govern interventions (described below). They do not influence the disease transmission or progression directly.

Interventions are external modifications of the state of the simulation where “external” means not governed by the disease model nor the contact network. An intervention comprises of a trigger, a target, and an action ensemble. The action ensemble, which is a set of instructions to apply to the dynamic state of the target (a collection of nodes and/or edges), is only applied when the trigger evaluates to 1 (or true). The trigger is a function of the system state and thus may depend on any of the above-mentioned attributes, including the Person Trait Database. The operations of an action ensemble are partitioned into the following three categories: (i) those triggered once per intervention (e.g., to update global variables), and which are thus independent of the intervention target, (ii) those that are applied to each element of the intervention target, and (iii) those applied to a sampled subset of the intervention target and optionally to the complement of the sampled subset. Moreover, sampling may be nested, thus permitting chained sub-sampling. Each individual action of the action ensemble is given a delay specifying the time (or offset) at which the action should execute, which permits flexible, fine-tuned scheduling within the simulation.

Before starting the simulation, the state of all individuals P must be initialized. In EpiHiper, initialization is simply a special case of an intervention where the trigger is omitted (and by convention is true).

Digital apps to estimate network properties

Digital contact tracing apps have recently been deployed in several countries [2, 5, 17, 25, 35] with mixed success. Reasons for the range of outcomes include: (i) low penetration levels, (ii) compliance, and (iii) accuracy of the apps in discovering neighbors accurately. Our allocation method is based on exploiting simple network properties that can be estimated using digital devices. We use two measures here: (i) degree and (ii) weighted degree. Digital contact tracing apps can potentially measure both of these quantities quite accurately.