An implementation of a multilayer network model for the Covid-19 pandemic: A Costa Rica study

1.   Introduction

Throughout the Covid-19 pandemic, mathematical and statistical models have been implemented to understand disease transmission, forecast the spread of the disease, and plan for strategic sanitary measures to minimize the burden on the public health system ^{[1,2,3,4,5,6]}. Deterministic models have been widely used to model infectious diseases in the past ^[7,8,9]. However, due to human mobility, unpredictable human behavior, and other factors, adaptable models are better suited to evolve with the ever-changing behavior of the population during the pandemic as sanitary measures are implemented in a particular region of interest ^[10,11].

Agent-based and network models have been extensively used to study the transmission of SARS-CoV-2 and to determine optimal and feasible measures to diminish the impact of the Covid-19 pandemic. For example, in ^[12,13], the authors estimate where, when, and how many transmission events happened in urban areas during the first and second waves of the pandemic for particular regions. Similarly, an agent-based model that describes age-specific risk factors for exposure and transmission is analyzed in ^[14], where they include different layers in their model. Moreover, ^[15] studies a stochastic, spatially structured individual-based simulation with other spaces (such as households, schools, workplaces, and the wider community) and age structure for an influenza pandemic. In general, agent-based and network models allow for a level of flexibility in which different scenarios can be formulated to provide valuable input to health authorities for the decision-making process; see, e.g., ^{[12,13,16,17,18]}.

Recent studies on Covid-19 with similar methodologies include ^[19,20], but their experiments are for networks with at most $110,000$ nodes. In ^[20], authors used six layers and an SEIR model for Brazil and concluded that the increase of isolation is the best option to keep the situation under the healthcare system capacity by removing layers and reducing social contacts. In ^[19], authors tested several what-if scenarios on school openings during the vaccination campaigns in France, with a multilayer network similar to ours (schools instead of social contacts). Moreover, in ^[10] authors constructed a two-layer multiplex network for the coupled spread of disease and conflicting opinions to determine how negative/positive opinions affect the transmission. In ^[11], authors analyzed the risk of the Covid-19 importation to Mexico through the Air Transportation Network with a multilayer network approach.

One of the main advantages of a network model is its flexibility when considering individual and social behavior, which allows us to mimic sanitary measures and evaluate possible impacts on the population. As the pandemic progresses, it is possible to mimic early or ongoing behavioral changes such as reductions in community interactions, mask use, mobility restrictions, network contacts, or lockdowns. Timely access to these types of tools contributes to a better reaction in case of emergencies. Moreover, we believe it is essential to share the experience in constructing such models for practical purposes and future studies.

In this paper, we consider a multilayer temporal network (see ^[21,22]) capable of efficiently managing large populations. A probability of infection is independently computed daily for each individual according to its interactions, permitting a parallel implementation for real-life applications. We count daily interactions per individual and determine if the virus is transmitted depending on the attributes of each interaction.

The model we present in this paper has been taken into consideration by health authorities in Costa Rica, with slight modifications related to the use of initial conditions for infected individuals and the parameters related to public health policies pertinent to the country. The complete model includes the entire population (approximately five million individuals) and the 81^* counties of the country, and considers different scenarios with gradual re-openings, mobility restrictions, different percentages for mask use, and short and long-term projections of the disease. Results were delivered to the public health authorities as part of the decision-making process ^{[23,24,25,26]}. It also has been used to project the impact of Covid-19 variants and vaccination strategies; see ^[27].

^*There are 82 counties in the country. However, the 82nd was established in 2018, and no historical census data exists to model mobility.

In this work, we provide insight into the development of multilayer network models with an easy-to-follow implementation that allows generalizing the transmission of infectious diseases. We also include a Matlab app with the algorithm implementation that can be found in ^[28], where several attributes and parameters can be modified, and the Matlab code.

2.   Materials and methods

2.1.   Network model

Consider a population that lives in a fixed number of counties. Mathematically, we consider a temporal three-layer network (a set of nodes that are connected by edges) where each node represents one individual; see ^[10,22]. We assume that all individuals are present in all layers and we couple layers by connecting nodes that correspond to the same individual; see the mathematical definition of a multiplex network in [29,Section 5.3]. Intra-layer edges correspond to edges within a layer, and it represents a possible interaction between two individuals. Each layer is an undirected, unweighted, simple graph, an on each layer, each node can have a different number of neighbors, which can change over time.

Interactions are classified into three layers, which represent different types of contacts with particular characteristics that account for the virus transmission. These layers are: (1) a household network (people that live in the same house), (2) a social network (known contacts such as friends and colleagues), and (3) a sporadic network (strangers that you may encounter in short periods when you visit random locations). These layers are randomly generated for each simulation; layers 1 and 2 are fixed, and layer 3 can change for each time step since usually a node has no control over sporadic encounters.

Mathematically, each layer (indexed by $\alpha\in\{1, 2, 3\}$) is a graph $G_\alpha = (V, E_\alpha),$ composed by the set of nodes $V$ that represents the population (which is the same for each layer), and the set of intra-layer edges $E_\alpha$ that connects nodes within layer $\alpha$. In this case, the full information about the multiplex network is encoded in three distinct adjacency matrices $a^{[\alpha]}$; see [29,Section 5.3.2].

Each node on the graph includes several attributes. We classify them into two classes: (1) fixed attributes, such as the county, household members, age group, degree, and graph connectivity for layer 2 (social network), and (2) variable attributes, such as connectivity for layer 3 (sporadic network), epidemiological state, number of days at current epidemiological state, number of interactions with social contacts per day, mask use, and self-care behavior. The former are defined at the beginning of each simulation, and the latter can take different values every time step.

For the first layer, it is assumed that we know the total number of households per county and the average number of individuals per household per county. We then group all the individuals into households. We assume a complete graph for each family and no edges between different families. It is possible to randomly assign an age group to each family member if data is available to include differences in interactions accordingly.

In this implementation, we create all the households and assign their inhabitants by using a Poisson distribution with a mean equal to the average number of individuals per household per county. Each individual is labeled with an ID from $1$ to $N$, where $N$ is the size of the population. We create the cells $\texttt{IDtoCounty}$, $\texttt{IDtoFamily}$ and $\texttt{IDtoAge}$ that return the county, family and age group of a given ID, and the cells $\texttt{familyToIDs}$ and $\texttt{countyToIDs}$ that return the IDs for a given family or county, respectively. These cells allow us to access the information required to build the graph directly.

The degree of each node in layer 2 (social network) is fixed for each simulation. Nevertheless, a subset of edges is selected on each time step because only a small number of daily contacts is expected to occur. As mentioned before, layer 3 (sporadic network) changes daily. In order to keep track of the transmission of the disease, we only need to include the edges between new exposed individuals $E$ and their contacts.

For the second and third layers, choosing the contacts of a node requires some assumptions. It is natural to expect different interactions depending on age, location, density, socioeconomic factors, and more. In our implementation, for every node, we need to: (1) define the degree of the node (number of contacts) per layer, given by a uniform distribution on an interval with endpoints equal to the minimum and maximum number of allowed contacts per county per layer, and (2) choose its contacts randomly from different counties based on a given connectivity distribution between counties.

We choose layer 2 contacts for new nodes as follows. We first create a list of eligible nodes based on a connectivity distribution between counties. For instance, density, demographic information, or mobility data could be used for this purpose. For simplicity, we consider a binary matrix where its entry $(m, n)$ is equal to $1$ if an individual from county $m$ can have contacts from county $n$, and $0$ otherwise; see a toy example in Table 1. This means that the eligible list for a node includes all IDs from counties that have connectivity with the node's county.

We then randomly choose the predefined number of contacts for each node from the eligible list, according to its degree. Further restrictions can be added when defining the eligible list, such as age group, density, and distance between counties. Edges on layer 3 are updated daily similarly, excluding preassigned contacts in layers 1 and 2. We remark that different approaches can be considered for the connectivity between counties, such as probability distributions as in ^[30], the use of mobility networks from mobile phone data as in ^[31], or by defining activity potentials as introduced in ^[32].

2.2.   Modeling the disease

The propagation in our model is based on an SEIR type model with seven compartment states defined by the following epidemiological variables: susceptible ($S$), exposed ($E$), diagnosed or observed ($O$), undiagnosed or not observed ($U$), hospitalized ($H$), recovered ($R$), and dead ($D$) states; see Figure 1. We consider a discrete-time stochastic model, for which each individual belongs to one compartment, and can move to a feasible next state on a daily basis, according to probabilities that depend on interactions with infected contacts (for moving from $S$ to $E$) or the number of days at a current state (for any other transition).

A susceptible node $i$ in layer $\alpha$ may become exposed if it contacts an observed or undiagnosed node on its network, depending on a probability of infection $p_{i\alpha}$. Once exposed, it becomes observed or undiagnosed, which might require hospitalization (including intensive care units). Eventually, a node can recover or die. We assume that a fixed fraction of diagnosed cases does not isolate and that hospitalized individuals cannot infect susceptible nodes. Modifications with additional states and transitions are straightforward to include.

For simplicity, we assume an initial graph given by all the exposed individuals at time zero and their contact networks. If data is available, we can include chains of transmission for a fixed time (depending on data) as an initial condition by creating edges among identified clusters (nodes belonging to a chain of transmission). Data includes anonymous patient ID, ID from an individual that infected him/her (if it was not an imported case), relationship (that we interpret as one of the three layers), county, age group, and date of diagnosis. This information is translated into the multilayer network, applying it to the model's randomly created families and social networks. Before community transmission in Costa Rica, the public health authorities were able to keep track of chains of transmission for approximately four months. This data allowed us to create accurate contact networks for the early stages of the pandemic.

We model the transmission of the virus in the multilayer network by computing probabilities of infection $p_{i\alpha}$ for interactions between a susceptible node $i$ and its infectious contacts. We consider their contacts via the different layers for a given initial set of exposed nodes; adjacent nodes of infectious nodes can become exposed. We then increase the network by including new exposed nodes and their contacts. In this way, we are considering a growing network that increases daily.

If the probability that node $j$ infects a susceptible node $i$ at layer $\alpha$ at a given day is $\beta_{ij}$, then the probability of infection $p_{i\alpha}$ for node $i$ at layer $\alpha$ is given by

where $j$ includes all the indices of nodes that can infect node $i$ that particular day at layer $\alpha$. Mathematically, following the notation introduced in [29,Section 13.2], we can write

since we count daily interactions between a susceptible node $i$ and its exposed contacts. Let $X_{i\alpha}(t)\in\{0, 1\}$ be a variable that denotes whether a susceptible node $(i, \alpha)$ is infected ($1$) or not ($0$). We can write then

The value $\beta_{ij}$ may depend on the layer, epidemiological state, social distancing, mask use, and self-care behavior, among other factors. Transmission of the virus can also include imported cases where infection occurs in an external location. Non-periodic events can also be considered (as massive events or the creation of clusters due to super spreading events) to study their impact and possible scenarios, where different probabilities $\beta_{ij}$ can be defined for each particular case.

Changes in social behavior are a tremendous challenge when modeling the transmission of infectious diseases. We have included the following ideas, allowing us to mimic the effect of different measures that were adopted in several countries:

$1.$ Bubble size: in periods of strict measures, such as lockdowns or mobility restrictions, we assume a reduction in the number of daily contacts for each individual in layers 2 and 3 (separately), which can depend on its county, emulating differences at least between rural and urban areas due to different population densities.

$2.$ Mask use, which was mandatory in Costa Rica on July 2020.

$3.$ Social distance between individuals, which depends on the interaction layer.

In this paper, we consider two binary random variables $\texttt{useMask}$ and $\texttt{selfCare}$ that indicate whether or not an encounter between two nodes includes mask use and self-care, respectively. We then have different values for $\beta_{ij}$ depending on the interaction between the two nodes and their attributes. We only need to compute $p_{i\alpha}$ for every susceptible node $i$ that has infectious contacts. Once $p_{i\alpha}$ is determined, a uniformly distributed pseudo-random number $r_i \in (0, 1)$ is generated. If $r_i < p_{i\alpha}$, node $i$ is marked as an exposed individual at the current time step.

Any other event between epidemiological compartments depends on a transition probability. Every node has a different probability of moving to a feasible next state (according to Figure 1) depending on the period of time at its current epidemiological state. If such an event should happen, we then store the new state and day; if not, the state remains unchanged. Similar expressions as (2.2) can be written, but for the sake of brevity we omit such equations.

To explain how data is stored and read, we show a toy example of transition probabilities for the diagnosed class $O$ in Table 2. For this example, a diagnosed individual $O$ has a $40\%$ probability of being hospitalized on day 1. Otherwise, it will remain in the same compartment on day 2. On day three, there are three possibilities: remain at $O$ ($50\%$), move to $H$ ($10\%$) or recover $R$ ($40\%$). Note that each row adds up to 1 (one state should always be chosen according to these probabilities), and eventually, all nodes that have been exposed should end up as $R$ or $D$. This approach allows us to include variable stay periods in each class if data is available; average stay periods can be considered otherwise. Pseudocode for the modeling of Covid-19, as described previously, is presented in Pseudocode. 1.

2.3.   Fictitious case: effect of behavior changes

We first present a simple experiment to explain how data is used and read in the model and to measure the effect of the variables $\texttt{useMask}$ and $\texttt{selfCare}$. We consider a population of one million individuals distributed in 15 counties; its demographic distribution and the connectivity matrix are shown in the appendix. We consider initially $E_0$ exposed nodes randomly chosen. Parameters used in these simulations appear in Table 3. We assume that $70\%$ of the population wears a face mask, $35\%$ respects self-care, and $60\%$ of exposed nodes become diagnosed. Values for the probability of infection are shown in Table 3 for layer 2; probabilities are halved for layer 3. Finally, the transition probabilities between states are shown in the appendix. For this particular set of parameters and population, the average running time for a simulation with a serial code is about 40 minutes. We then impose restrictions on social behavior by increasing the percentage of individuals wearing a mask and decreasing the number of daily interactions over two months, starting on day 20, emulating the effect of mild regulations.

2.4.   Costa Rica case study

In this section, we present an application of the model for Costa Rica. As stated earlier, we include chains of transmission for a fixed time as an initial condition by creating edges among identified clusters according to data from the Ministry of Health. To give an insight into reported cases in Costa Rica, we illustrate how diagnosed (confirmed) cases evolved from March 8, 2020, to June 24, 2020; see Figure 2. For simplicity, we only show edges between reported cases and not their entire networks and include a random position for each node within its county.

Regarding the first layer, families are created first from data from the "National Household Survey" ^[33]. Layers two and three are constructed as described above. By this time, the degree of nodes in layer two was taken as a uniform distribution in the interval ^[5,20], and daily, a subset of contacts is considered, depending on the county of residence (greater values for the metropolitan area and lower values for rural zones). The degree of nodes in layer three was taken as a uniform value in the interval $[0, 20]$; see, e.g., ^[20,34]. Connectivity between counties for Costa Rica is obtained from ^[33], where a matrix with the household and workplace is available.

Most of the parameters for transition between states are essentially provided by the Costa Rica Social Security Fund, related to the average number of days in each state, percentages of individuals that require hospitalization or ICU, and mortality rates. On the other hand, we assume in this particular example that 10% of diagnosed individuals do not isolate (still have daily interactions with their contacts) and that there is a relation of 75%/25% for observed and undiagnosed cases in the population. The probability of infection for each encounter between nodes given in (2.1) is modified by considering a calibration parameter $\gamma \in (0, 1]$ such that

For approximating $\gamma$, we considered a period where we calibrated the model by using a standard bisection method for finding an optimum value for $\gamma$, such that the square error between observed and predicted cases is minimum. The calibration period includes the real number of observed cases according to data provided by health authorities. The value of $\gamma$ is then fixed for all simulations to obtain short-term projections of Covid-19 cases. An Approximate Bayesian computation was used to estimate the initial value of the transmission rate $\beta$ for the beginning of the pandemic in Costa Rica considering an SEIR type model; see more details in ^[35]. We then compute $\beta_{ij} = c_{ij} \beta$, where $c_{ij}$ is a reduction parameter that depends on the use of mask and self care for each interaction; values for $\beta_{ij}$ are as shown in Table 3.

3.   Results and discussion

We have presented an implementation for a multilayer network to support public health authorities during the Covid-19 pandemic. Capturing the transmission dynamics of the SARS-CoV-2 virus is essential to provide information about disease projections and specific data-driven scenarios. We remark that in some cases, it may be important to consider the whole network, for example, if the probability of contact varies among susceptible individuals or in modeling vector-borne diseases where vectors play an essential role in the transmission.

The model exhibits a flexible structure that allows us to include particular attributes for each node and contact network, describing the complexities of social behavior and the specific disease. In this article, the implementation of the general structure of the model is described, which can be used in other social and epidemic contexts.

We first present results for 1,000 simulations in Figure 3 for a fictitious example, where we compared the effect of mild restrictions as described above. We observed a reduction in the peak of hospitalized individuals, essential to reduce the stress on health systems, and a reduction in the peak of observed cases. We also showed that uncertainty among simulations is insignificant, even though the model is not deterministic. This experiment gives insight into how different measures can be translated into effective ways to flatten the curve. Moreover, we observe that even though the final value of cumulative cases is slightly reduced by day 140, the relevance of these results is related to the fact that the peak on the hospitalized individuals curve is reduced, allowing to reduce the stress on the public health system. It is relevant that this kind of result relies on the availability of data to compute times in each compartment. The model is sensitive to the proportion of diagnosed/undiagnosed individuals, but the effect on hospitalizations is expected to be negligible since most undiagnosed individuals have no symptoms.

Regarding the case study of Costa Rica, in Figure 4, we present some numerical results (500 simulations) where we took five days for calibration (we created the network until June 18, 2020), optimized the value for $\gamma$ according to observed cases between June 19 and June 24, and then projected results for two weeks, until July 7, 2020. We compared the projected cumulative number of cases and hospitalized individuals with real data after June 24, 2020. We observed how the model captures the expected growth in the cumulative cases in a short period and the tendency in the number of hospitalized individuals. We observed a small relative variation in all the simulations. We refer to ^{[23,24,25,26,27]} for further scenarios and results of this implementation. We remark that data before June 24, 2020, is exact for cumulative cases, and they are used to calibrate the probability of infection for the following days. Short-term projections are accurate because a five-day calibration window is used. Retrospectively, we observe that the simulated and real data are a good fit for observed and undiagnosed individuals. In the case of hospitalizations, we used average stay periods for hospitalized individuals according to data from the Ministry of Health. Moreover, including age groups has helped to improve results for hospitalizations, and in general terms, the behavior is captured, and our projections provide a worst-case scenario.

Compared to recent studies on Covid-19 with similar methodologies such as ^[19,20], our temporal network accommodates larger populations. Furthermore, our results mimic some restricted collective behavior changes in the population and can include further details regarding sanitary measures.

Implementing and executing these models can be challenging. The availability of information is closely related to the specificity of the question that needs to be answered. Information on the initial conditions for the model and parameter value estimates can be a limiting factor when data is unavailable, or its quality is not optimal. Efficient and fast algorithms are necessary due to the on-demand health scenarios that must be simulated frequently, and computational resources limit execution times.

In a rapidly evolving epidemic such as Covid-19, the model must be calibrated frequently, as social behavior is constantly changing and, in many cases, unpredictable. The model discussed in this paper allows us to adapt and include frequent changes in variables that mimic social behavior. Additionally, the model can be generalized to have other factors, such as the effect of herd immunity and the role of vaccination, as countries take the necessary steps to get back to the new normal.

Public health authorities were forced to make short-term decisions during a health emergency such as the one caused by the SARS-CoV-2 virus. Mathematical models help authorities systematize the virus transmission mechanisms and project the expected short- and medium-term behavior under some assumptions. In particular, a multilayer network model is a valuable tool as an input in the decision-making process for health authorities. It provides the flexibility to include social factors and an epidemiological structure, enabling a greater level of detail in the study of the disease. However, despite the model's flexibility, limitations mostly pertain to data availability and quality and the ability of health authorities to execute policies based on scientific evidence.

Acknowledgments

The authors would like to thank support by the Research Center in Pure and Applied Mathematics and the Department of Mathematics at Universidad de Costa Rica.

Conflict of interest

The authors declare there is no conflict of interest.

Supplementary data

Data for Example 1

We present the fictitious demographic data and parameters used in the first experiment; see Tables 1, 3 and 4. Transition probabilities between epidemiological states used in the toy example depend on the age group; see Tables 2, 5, 6, 7, 8.