Real-time Updating of Dynamic Social Networks for COVID-19 Vaccination Strategies

2.1 SIR model

The analysis of the diffusion is conducted using a standard SIR model [4] with an additional state describing the number of vaccinated people, as shown in Fig. 1. For each individual, S, I, R denote the susceptible, the infected and the recovered (patients who are not infectious anymore). The SIR assumption has been widely adapted to simulate COVID-19 propagation [16, 61, 60] since reported COVID reinfection cases (e.g [59]) are still rare compared to the total number of reported cases thus far. The SIR model has also been broadly used in network-based disease simulations via random-walk-based simulations [35]. Each node symbolizes an individual in the social network, whose status can alter from susceptible to infected (S-I), or infected to recovered (I-R), according to the random walk through temporal edges [21]. The transition from susceptible (S) to vaccinated (L) only takes place when required according to chosen vaccination strategies. In contrast to classical disease modelling, since recent research[43] shows that current COVID vaccinations can be significantly less effective when facing new variants (e.g.,B. 1.617. 2 (Delta)), the L-S and L-I transitions can be activated as shown in Fig. 1. More details about the transition probabilities are given in Sect. 2.2. In view of the fact that until these days the infection probability after vaccination is still unclear, L-S and L-I transitions are not considered in this study. Nevertheless, the developed model can easily incorporate these types of transitions when required.

• Download figure
• Open in new tab

Figure 1:

Illustration of network-based SIR model with a vaccination state L.

2.2 Graph-based vaccination strategy

Both disease spread simulation and optimal vaccination modelling based on social networks have been receiving increasing interest for different types of infectious diseases [51]. We consider an undirected graph 𝒢 that is a pair of sets 𝒢 = (V, E), where V = {v₁, v₂…v_n} represents the set of individuals (graph nodes) and the set E contains the edges, each connecting a pair of individuals. Each graph edge e ∈ E is represented by a triple e = (v_i, v_j, w_i,j) where v_i, v_j are the two endpoints and w_i,j ∈ ℝ is the edge weight. For unweighted graphs w_i,j ∈ {0, 1}, while for weighted graphs w_i,j could represent the frequency or the intimacy of the contact. In epidemic spread modelling, the infectious probability p_i,j from the individual i to j (and vice versa) is often in function of w_i,j, p_i,j = ℐ 𝒫 (w_i,j). We also note that p_i,j may depend on individual-level characteristics of v_i and v_j, such as age or activities. The connecting graph can be fully represented by the associated adjacency matrix A = {A_i,j }_i,j=1,…,n. We use three Boolean vectors {I_t, L_t, R_t}∈{ { 0, 1} ⁿ}³ to indicate the status of each individual, either infected, vaccinated or recovered in the SIR model, at time t. The recovery period T_γ ∈ N is an uniform distributed random variable generated individually for each individual.

If we adopt the edge-wise function ℐ 𝒫(.) in the whole network, the infectious probability vector at time t in this SIR model reads where 1_n = [1, 1…1]^T and ⊙ denotes the vector-wise Hadamard product.

Following a uniform probability distribution, the vector of infections I_t is simulated using and I_t−1. The only controllable variable in Eq. 2 is the vaccination vector L_t.

Different graph-based vaccination strategies can be employed to enhance the immunization impact with a limited vaccination capacity. The state of the art approaches are usually determined by observed individual- or community-level social connections, often involving classical graph measures, for instance, graph degree, betweenness centrality [26] or community links[12]. Much efforts have also been made to use these strategies in practical settings where significant positive impacts have been observed [31]. Since the available graph data often include non-negligible uncertainties (missing vertices or edges), statistical models are commonly employed to provide an optimal estimation of these graph measures. Practical approaches involve, for example, fixed choice designs (FCD) [45] and the nomination strategy [23], both based on an estimation of the graph degree. Even with partially observed dynamic networks, the vaccination strategy could be significantly improved in terms of reducing the maximum infected number and delaying the disease propagation, compared to a random choice [63]. Nevertheless, precise knowledge of the network structure is crucial to determining an efficient vaccination strategy. It is essential to use community-based approaches (e.g [30], [12]), since graph clustering algorithms can be sensitive to noises. However, the data collection of dynamic social networks remains cumbersome, especially for large dimensional problems. In this paper, we conducted our analysis based on three classical strategies, considered less sensitive to data noise, compared to community-based approaches.

3.1 Variational assimilation and BLUE

DA algorithms aim to combine different sources of noisy information in order to provide a more reliable estimation of the current system. The state variables could be either a physical field or a sequence of parameters. The true state, denoted by x^true, stands for the theoretical value of the state at some given coordinates/time, often out of reach in real-world applications. The objective of the assimilation is to gain an optimal approximation x^a of the true state x^true, based on the prior information which are two parts: an initial state estimation x^b (so-called the background state) and an observation vector y. The former is often issued from prior numerical simulations/predictions while the latter can be obtained via physical measures of some control variables. Their tolerances, regarding theoretical values, are quantified by ϵ_b and ϵ_y, where the observation operator ℋ from the state space to the observable space is supposed to be known. The probability distributions of the prior error are supposed to be centred Gaussian, characterized respectively by the covariance matrices B and O.

The key idea in variational methods is to find a balance between the background and the observations using maximum a posteriori (MAP) method [13]. This leads to the loss function weighted by the inverse of B and O, The optimisation problem defined by the objective function of Eq. (7) is called three-dimensional variational method (3D-VAR), which can also be considered as the general equation of variational methods without considering the transition model error. The output of Eq. 7 is denoted as x^a, i.e. If ℋ can be approximated by some linear operator H, Eq. 8 can be solved via BLUE (Best Linearized Unbiased Estimator) formulation, where P_A = Cov(x^a − x_true) is the analyzed error covariance and the K matrix, given by is so called the Kalman gain matrix. In the rest of this paper, we denote H as the linearized transformation operator. The case when ℋ is non-linear is more challenging for finding the minimum of Eq. 6, especially for high-dimensional problems. The resolution often involves gradient descent algorithms (such as “L-BFGS-B” or adjoint-based numerical techniques).

3.2 Online assimilation with graph data

The essential idea is to perform real-time updating of the partially observed dynamic networks based on other available information, such as sub-graph structures or the current number of those infected. To this end, the prior observed network at time t is considered as the background state (i.e.,), while other information is embedded in the observation vector y_t.

Once the current contact network is updated based on Eq. 7, vaccination strategies can be implemented on the analyzed network (i.e.,step 1 → step 2 in Fig. 2) which is a more accurate approximation of the true state. The degree and the betweenness centrality of the assimilated network is given by

• Download figure
• Open in new tab

Figure 2:

Illustration of real-time DA updating for partially observed contact networks

where denotes the element (k, v) of the adjacency matrix . Similar expressions of and g^b(v) on the background state can be given using A^b and . The principle of real-time assimilation with graph data is illustrated in Fig. 2 where the virus propagation is simulated using the SIR model, as described in Sect. 2.2 between two vaccination steps. Compared to the overlapped graph, the advantage of working with temporal networks is that the temporal correlation could be considered. In fact, an individual can be active for a relatively short period of time only, as shown below in Sect. 4.1. Therefore, instead of using an overlapped graph (if available), analysing temporal networks can result in an efficient real-time vaccination strategy.

A major challenge of implementing DA algorithms with graph data is the computational cost since the adjacency matrix A_t, considered as the state variable, is a two-dimensional vector. We can rely on the assumption of graph sparsity and appropriate parameterization to reduce the computational burden. In this work, we propose two DA frameworks for dynamic networks updating, respectively introduced in Sect. 4 and 5. The former aims to reconstruct the full network with observations of sub-graphs, while the latter attempts to adjust the parameterized community-wise infectious probability, relying on multi-layer modelling. These two modellings, relatively at the local and global scale, also show the flexibility of this data assimilation framework.

4.1 Assumptions and preliminary analysis

This study is based on recently (before the COVID outbreak) collected face-to-face contact data from a French high school [29], which has been used to simulate a COVID outbreak [44]. The connection networks of 329 students (coverage of 86% of the students) in a high school in Lyon are available for 7374 time steps in a week. For the sake of simplicity, we condense the dynamic graph to 78 time steps by overlapping every 100 consecutive networks. Each time condensed time step symbolizes 30∼ 60 minutes. The temporal networks remain sparse since the average graph density (i.e. number of non-zero edges divided by the number of node pairs) is equal to 0.76%. All contact networks are assumed to be undirected, which means the associated adjacency matrices are all symmetric (i.e.,) and the virus could spread in both directions of an edge. According to [44], the infectious probability (of a 20-second contact) in this network can be estimated as p ≈ 0.1%∼ 1%. However, this estimated probability might be contested for the newly discovered SARS-CoV-2 variants [32]. In this paper, in order to adequately investigate the optimality of different vaccination strategies, we fix the infectious probability to p = 2%. The average recovery period in the SIR model is set to 60 time steps (around 4 to 5 days), following a uniform probability distribution, i.e. T_γ ∼ unif(55, 65).

We begin by performing some preliminary analysis of the network data in order to better understand the underlying graph structures. The overlapped network (i.e. ) of all the time steps is shown in Fig. 3(a) where a clear community structure can be observed. Identifying these communities is crucial to simulating the disease spread [36], especially for a highly infectious virus like SARS-CoV-2 [55], and to determining optimal vaccination strategies. Much effort has been given to developing community-detection algorithms in social networks [1, 53]. In this work, we make use of the Fluid community detection algorithm proposed by [53], which is advantageous for sparse graphs since the algorithm complexity is linear to the number of non-zero edges in the network, i.e. 𝒪(|E |).

• Download figure
• Open in new tab

Figure 3:

Preliminary analysis of the dynamical high-school connection network: (a): Overlapped contact network. (b): Performance rate p^r(𝒞) against assumed community number (c): Reordered adjacency matrix after clustering (d): Node degree distribution of the first and the last 50 time steps.

In real applications, specifying the number of communities is usually difficult. Here, we apply several times the community detection algorithms against different assumed community numbers k_c, before evaluating the performance rate p^r(𝒞) [25] of the obtained partition 𝒞. The latter is defined as where indicate the number of edges of intra- and inter-clusters respectively. The performance rate is commonly used as an indicator for finding the optimal community number, which is a standard approach for graph clustering problems. According to the result presented in Fig. 3(b), where we clearly observe a stationary performance rate starting from k_c = 4, we choose to proceed with the optimal number of clusters . The final clustering result is displayed in Fig. 3(a) where clusters/communities are shown in red, green and blue. The three detected communities are equivalently distributed, as shown by the reordered adjacency matrix (Fig. 3(c)), with 106, 110 and 111 nodes respectively. From a practical perspective, these communities could be considered as different grades or classes in the high school, with a similar structure to the graph data presented in [27].

4.2 DA modelling and numerical results

Since it is infeasible to collect contact networks via wireless equipment in all educational settings post lockdown, the objective of this study is to enhance the vaccination strategy when only partial/noisy information is available, for instance, via tracing applications. For this reason, the full contact networks are supposed to be out of reach. In terms of background states and observations, we suppose that the temporal network is only partially observable a priori where 50% to 70% of nodes are missing in the background estimation of the network . The missing nodes are selected randomly and kept invariant at all time steps. In reality, the missing nodes could refer to, for example, people who haven’t installed the tracing application on their smartphone. We also use an observation vector y_t, which contains the sub-networks for each of these three detected clusters. Thus, we suppose that the intra-community contacts of students in each class/grade are fully observable with y_t. The objective is to perform DA algorithms sequentially to correct the knowledge of the background network relying on the observed sub-networks. The transformation operator H is thus linear (sub-Identity matrix) and the DA problem is solved via BLUE, as shown in Eq. 7. and y_t are vectorized with Identity error covariances B and O, as demonstrated in Fig. 4.

• Download figure
• Open in new tab

Figure 4:

DA process at the current time where intra-connections are observable

After each vaccination, the SIR model is applied to simulate the virus propagation until the next time step, as summarized in Eq. 2. An essential advantage of BLUE-type formulation with invariant prior covariances is that the Kalman gain matrix can be computed offline a priori since it is invariant to the current x_b and y. The computational cost of DA can thus be considerably reduced. The vaccination capacity is fixed 2%(= 6 individuals of all students for all strategies (random, highest degree, highest centrality) presented in Sect. 2.2, based on prior or assimilated graphs.

The evolution of the number of infected |I_t|, according to different vaccination strategies, is displayed in Fig. 5, where the percentage of missing nodes in the background state is fixed as 50%, 60% and 70% respectively. To acquire robust numerical results, each type of simulation with or without vaccinations is repeated 10 times and the average values are drawn in solid or dashed curves in Fig. 5. Standard deviations of the simulations (except dashed lines) are also displayed in transparent shades to ensure the robustness of the comparison. The averaged maximum number of infected for each strategy is shown in Table 1. We note that vaccinations take place at every time step for 6 selected students (≈ 2% of the population) after the simulation of virus propagations with a infectious probability of 2% for each temporal edge. The initial infected I_t=0, commonly used for all simulations, is randomly simulated with a probability of P ((I_t=0)_k) = 10% for k = 1, …, 329.

• Download figure
• Open in new tab

Figure 5:

Evolution of infected against different prior error level (percentage of unobserved vertices): (a) 50%, (b) 60%, (c) 70%. Standard deviations are also displayed by transparent shades.

From Fig. 5, we observe that almost all averaged curves rise to a high point and peak around t = 50− 60 when all individuals are either infected or vaccinated. Since the vaccination process takes place in a relatively short period (a week), we suppose that the infected individuals are not detected in real-time. As a consequence, a student can be vaccinated after being infected by the virus, leading to vaccine failure. This fact emphasizes the importance of the vaccination strategy chosen. What can be clearly observed from Fig. 5 is the decreasing infected number according to the vaccination strategy in the order of free (no vaccination) → random →background→ assimilated (DA). This order is globally consistent regardless of time. First, all vaccination strategies manage to significantly reduce the number of infected and delay virus propagation compared to the free simulation (green curve). In terms of maximum infected number, for all three cases, the peak value is reduced on average by 26%, 34%, 34%, 40% and 37%, respectively for random, background with highest degree, background with highest centrality, assimilated with highest degree and assimilated with highest centrality. All other strategies are dominated by the assimilated curves, especially when proceeding with the highest degree strategy. The difference, in particular between background and assimilated curves, is more significant when working with large-scale networks. On the other hand, for background-network-based strategies, a growth of maximum infected number against prior error level is noticed in Table 1 while the results based on assimilated networks remain robust. This fact promotes the use of data assimilation on network data when prior error level can not be precisely specified. We note that the missing nodes at each time step are generated independently with no temporal correlation, explaining why reasonably good results can be obtained with 70% missing nodes. In summary, numerical results show that the DA-based real-time updating of networks considerably improves the impact of vaccination, resulting in reducing virus spread.

In these experiments, the use of node degree (solid curves) and betweenness centrality, for both background (red) and assimilated (blue) cases, exhibits a similar performance. Such fact suggests a high-level (non-negligible) inter-clusters connections where a contrary case can be found in Sect. 5.

5.1 Multi-layer modelling of scale-free networks

As stated in recent research [42], the infectious probability of COVID-19 can differ significantly for different populations, based on their age, gender, activities and so on. For example, both the transmissibility and the mortality rate is reported to be higher for aged people, necessitating appropriate strategies to protect this fraction of the population. SARS-CoV-2 variants may also vary geographically [5], leading to inhomogeneous transition probabilities. Since the outbreak of the COVID-19 pandemic, continuous effort has been made to understand the behaviour of the virus infection with respect to individual-level (e.g. aged people [48]) and community-level (e.g. healthcare workers [58]) characteristics. These phenomena have led to the idea of using multi-layer networks, where different types of connections exist between graph nodes (see Fig. 6(a)) to simulate the virus spread in social networks. In general, multi-layer (also known as “multiplex”) networks [18] are widely used to study graph diffusion problems [28] and define generalized versions of Pagerank [19]. Recently, multi-layer modelling has also been applied to COVID-19 spread simulation [57] where each layer refers to a potential contamination community, such as school, workplace or transport. Appropriate use of the information on these layers can optimise vaccination strategies as mentioned in [8], by prioritising the populations with high risk and high transmissibility.

• Download figure
• Open in new tab

Figure 6:

(a): Illustration of multi-layers network modelling where the infectious probability depends on the layer of the reception node. (b): Different layers in one temporal contact network where, for example, the yellow layer could represent the community of academic staff in the department of Computing at Imperial College London and the other layers stand for students of different grades (CI_a).

Since the collection of large-scale face-to-face contact multi-layer dynamic networks is extremely complicated, we rely on conceptual modelling in this work to further examine the performance of the novel approach. Dynamic contact networks of 1000 individuals and 5 layers (each of 200 nodes) are synthetically generated, where each layer suggests a specific group in the population, according to their age or activities (e.g. students, healthcare workers). Assuming all the edges in the temporal networks are fully observable, our objective is to calibrate the time-variant infection probabilities {p_i,t} _i=1,…,5 based on the observation of infected number in each of the layers {I_i,t}_i=1,…,5. The temporal variance of {p_i,t}_i=1,…,5 can be a consequence of SARS-CoV-2 mutations. More precisely, the values of {p_i,t}_i=1,…,5 update every 5 time steps, following a stochastic process, where δ_p,m∼ unif (−0.04%, 0.04%) and the observation vector consists of incremental infected numbers ΔI_i,t = I_i,t — I_i,t−1. For inter-layer connections, the infectious probability is determined by the layer of the receiving nodes, i.e. as shown in Fig. 6(a). It is worth mentioning that the associated adjacency matrix A_t is no longer symmetric under this assumption. Nevertheless, the network virus spread modelling in Sect. 2.2 remains valid.

As for the generation of temporal networks, we depend on the concept of scale-free networks [52] where the degree distribution follows a power law, where P_sf(k) stands for the probability of a node to have k connections while 2≤ γ≤ 3 is a chosen parameter. To simulate intra-connections in each layer, we use the Barabasi-Albert (BA) model [2], which is scale-free with γ = 3, incorporating two important concepts in graph theory: growth and preferential attachment [38], which exist widely in social networks. Therefore, the BA model is a reference tool to generate real-world-like networks, including web connections or citation networks. To generate a BA network, nodes are added to the network consecutively where the probability of the new node to be connected with the existing node v writes The denominator in Eq. 17 represents twice the current number of edges in the network. Individuals with a higher degree have a stronger ability to grab links added to the BA network, which is an adequate assumption for social networks. Moreover, the inter-layer connections are generated randomly with a density of 0.5%, much sparser than intra-layer edges. Eventually, an example of a complete temporal network is drawn in Fig. 6(b) where the five layers are shown in different colors.

Since temporal edges are supposed to be known in this modelling, we aim to estimate {p_i,t}_i=1..5 based on the evolution of the infected number in all five layers. In fact, we can predict Δ {I_i,t}_i=1..5 via a prior estimation of {p_i,t}, establishing a state-observation mapping H ∈ ℝ^5×5 for DA algorithms. The DA problem could be addressed as where The simulation/vaccination framework is similar to the one in Sect. 4 with a vaccination rate of ≈2% of the population at each time step. This means that all people will be vaccinated before t = 50. For all assimilations, the error covariances are set to be identity matrices, as in Sect. 4. Our goal is to determine an optimal vaccination order based on available noisy information. In order to cover more possible scenarios, we set various initial probabilities {p_i,0}, as shown in Table 2, denoted as CI_a, …, CI_f. For the sake of simplicity, {p_i,0} always follow a decreasing order from layer 1 to layer 5. Typically, the initial probabilities in CI_f are more homogeneous compared to CI_a or CI_e. To give an example, CI_a could be used to simulate, for instance, a scenario in the department of computing at Imperial College where nearly 800 students plus faculty members can be found. The layer with high infectious probability may consist of professors, (senior) researchers and HR officers, while the other four layers can represent graduate or undergraduate students of different grades. The former community has a much higher average age, in contrast to the latter. Furthermore, each community holds a dense intra-connections, coherent with our model assumption. The diversity of the initial conditions (CI_a, …, CI_f) ensures the robustness of the proposed approach.

The experiments set-up is similar to the one in Sect. 4. While computing the node degree and the betweenness centrality, the graph edges are weighted by either the background or the analyzed layer probabilities. Since the layer information is unattainable a priori, background networks are set to be homogeneous (i.e., ). The evolution of the infected number, issued from a Monte Carlo test of 10 simulations, is illustrated in Fig. 7. The stand deviation is represented by colored transparent zones. We also display the result of using exact {p_i,t} (instead of (red) or (green)) for vaccination in yellow. This curve is thus considered as the optimal target for the assimilation-based strategy. When vaccinating the nodes with the highest degree, a substantial advantage of the DA approach (solid green line) compared to the background one (solid red line), can be noticed in all 6 sub-figures of Fig. 7. In fact, both the maximum infected number and the average standard deviation have been significantly reduced, as confirmed in Table 3. On the other hand, DA has much less impact when selecting the individuals with the highest centrality, as shown by the dashed lines in Fig. 7. A reasonable explanation for this could be the phenomenon of brokerage [40]. The endpoints of the few inter-layer edges play an essential role in virus spread. These nodes, also known as “broker”, do not necessarily have a high degree in the graph. However, since many of the shortest paths pass by them from one layer to another, the betweenness centrality may peak at these nodes with or without adjusting {p_i,t}. This fact shows that when precise knowledge about inhomogeneous infectious probability is out of reach, proceeding with the highest centrality might be a robust choice. Nevertheless, both the dashed green line and the dashed red line are dominated by the solid green line (assimilated networks with the highest degree) in all 6 sub-figures.

• Download figure
• Open in new tab

Figure 7:

Evolution of infected number (average of 10 simulations) following initial conditions CI_a.. CI_f

We also note that for Fig. 7(a,b,d) where the five layers exhibit more variance for the initial probabilities, the assimilated curve is much closer to the optimal one. In fact, optimally vaccinating an inhomogeneous network requires less accurate knowledge of layer probabilities so long as the most infectious layers can be identified. For example, proceeding with (5%, 1%, 1%, 1%, 1%) and (7%, 0.5%, 0.5%, 0.5%, 0.5%) for vaccine priorities may lead to similar results.

The evolution of the normalized true layer probabilities is , while their posterior (analyzed) estimation is . The gap between the estimated and the true ratio of probabilities is rapidly reduced with the increasing of , which results in a more optimal vaccination strategy. Since vaccinating infected individuals is ineffective, the early phase (around the first 20 time steps) of the outbreak is crucial to delaying the COVID spread because the most active individuals (either in terms of degree or centrality) can be infected very quickly. Therefore, the DA correction at the start of the vaccination process plays an essential role in reducing the propagation speed. On another note, we also observe that a strong oscillation in the values of which implies high instability of the observation vector y_t = [ΔI_i,t]_i=1..5 due to sampling uncert ainties.

In summary, the assimilation-based vaccination strategy shows competitive performance in this multi-layer modelling even though the assimilated layer probabilities are just approximations. Using the assimilated temporal networks with “highest degree” dominates other approaches, with a smaller average infected number and lower standard deviation.