Investigation of SIS epidemics on dynamic network models with temporary link deactivation control schemes

1.   Introduction

Infectious diseases are unique from other conditions in that the direct cause of the disease is another infected individual. A dichotomy to classify such communicable diseases is whether they confer a long-term or short-term immunity upon infection. The condition will eventually die out after one cycle in the former case. Short-term immunity can be of more significant concern, as witnessed by recent COVID-19. Reinfection by mutated diseases can also complicate the situation. In general, with the absence of any intervention, infectious diseases with low transmission rates or high recovery rates will progressively disappear. In contrast, those with high transmission rates or low recovery rates will reach some steady state.

The presence of infectious diseases causes harm to human society in myriad ways, including placing a high burden on healthcare services and disrupting economic activities. A control method can be implemented to reduce socioeconomic costs to eradicate the transmissible disease even quicker. Also, if the epidemic is difficult to eliminate, a control scheme can lower the spread rate and reduce unhappiness.

Various interventions are used in the real world and can be divided into two main types. First is pharmaceutical interventions, including vaccinations and medical treatments ^[1]. The second is non-pharmaceutical interventions (NPIs) like social distancing. Both can be considered control methods that attempt to reduce the spread of disease. However, the distribution of antidotes can be another complicating issue among countries. Hence, a more readily approach may be through personal and community efforts, including curfews and lockdowns ^[2].

Finding an effective means to control the disease spread can be achieved by taking mathematical approaches ^[3]. The susceptible-infectious-susceptible (SIS) model is one of the most popular models for short-term immunity. A simplistic view is taken from the viewpoint of disease transmission such that each individual in a randomly mixed population is susceptible (S) or infectious (I). A susceptible individual becomes an infectious one upon infection and back to susceptible upon recovery. Various modifications of compartmental models exist where extra classes are added to reflect disease-specific processes, including the standard SIR model that confers immunity upon recovery and SEIR that accounts for infected individuals who are not yet infectious. The flow between compartments can be described in ordinary differential equations (ODEs) in all these models ^[4,5].

Although compartmental epidemic models are intuitive and straightforward, their legitimacy has often been questioned for their oversimplification, assuming all-to-all connectivity among homogeneously well-mixed individuals ^[6,7]. Modeling epidemics on networks is an alternative and novel approach that emerged from renowned works ^[8,9]. Modern network-based models with explicit network structures can describe the individual-level interactions in more detail and thus capture the observed dynamics with higher integrity.

A more sophisticated investigation of disease control methods also followed ^[5,10]. Dynamic networks with changing links allow variations within the network structure and offer new modeling options for studying the influence of control measures. It is in many ways similar to our intuitive perception of forming and breaking human relationships ^[10]. Over the past decade, modeling and analyzing disease spread on temporal links have advanced significantly as studies continue to reveal the nature of networks existing in the real world, such as biological, social, and technological networks, are indeed dynamic ^[5].

This study examines how well an actual spread can be controlled within a set time by simulating social distancing on a network. Healthy individuals can temporarily stop interacting with infected neighbors. The non-linear model predictive control (NMPC) method is used to compute control signals to bring the number of infected closer to the target. These signals are applied to individual-based models (IBMs) having the same epidemic parameters as the ODE system. We further explore how the system can converge to the desired state even faster.

Many idealized networks have been suggested for modeling epidemiologically relevant contacts within a real population. Here we use the Erdős‑Rényi (ER) network, where the degree distribution follows a Poisson distribution. Our results demonstrate that more frequent updates of control signals do drive the system to converge to the target under epidemics with higher transmission rates. We note that the success of control is contingent on system parameters, control parameters, and damping parameters. Also, network-based epidemic simulations are yet subject to drawbacks where the sensitivity of outcomes to the properties of networks is difficult to ascertain ^[10]. Our approach focuses on closing this gap and gaining an intuitive understanding of connecting control theory with modern network-based epidemics.

2.   Models

In this section, we present the setup of the ODEs and the stochastic IBM of the dynamic network model. An SIS disease spreads on a fixed network, in which links can temporarily deactivate ^[11].

2.1.   Pairwise ordinary differential equation model

Compartmentalizing individuals according to their disease status is the foundation of almost all mathematical infectious disease models ^[4,12]. The basic SIS model takes a simplified view of how disease propagates within a population by classifying individuals as either S or I. The assumption leads to the standard system of differential equations:

where $\left[S\right], \left[I\right]$ are mean numbers of susceptible and infectious individuals, respectively, where the total number of the population is $N$ ^[10]. $\lambda$ is the rate of infection, which depends on the expected number of infected individuals such that $\lambda = \tau n\left(I/N\right)$. Here, $\tau$ denotes the transmission rate per contact and $n$ is the total number of links in the system. $\gamma$ is the rate of recovery from infection. The system of ODEs describes the transition between compartments, which allows tracking the variation in the expected number of individuals in each compartment. In the simplest case, $\gamma$ and $\tau$ can be considered constant.

A pairwise model is an alternative approach that approximates the spreading process using linked pairs of individuals ^[13,14]. A system of pairwise ODEs consists of the expected number of links $\left[XY\right]$, which connects individuals of status $X$ and status $Y$ where $X$, $Y$ ∈ $S, I$. For instance, [SI] represents a connection between a susceptible and an infected. Furthermore, $\left[XYZ\right]$ indicates the number of triples, where $X$, $Y$, $Z$ ∈ $S, I$. Here, the central node $Y$ connects to the nodes of status $X$ and $Z$. Since the network is assumed to be undirected and symmetric, the relationship [SI] = [IS] holds. Also, every link shall be counted twice.

The system of pairwise differential equations for SIS epidemics on an undirected and unweighted network is shown in (2.2). A disease can propagate from an infected node to its susceptible neighbors with a constant rate $\tau$ per SI link. The infection rate of a susceptible individual is $\mathrm{k}{\rm{ \mathsf{ τ} }}$ where $\mathrm{k}$ denotes the number of infected nodes. An infected individual is assumed to recover independently, regardless of the status of its neighbors.

Triples in the system can be further simplified by adopting a moment closure approximation. This

assumption replaces each triple with singles and pairs as follows ^[10]:

In dynamic network-based epidemic models with a temporary link deactivation mechanism, each susceptible can adjust their contacts to their neighbors, whose infection statuses are assumed to be known at all times. The model consists of a deactivation rate $d$ at which SI links temporarily deactivate, and a reactivation rate $a$ at which deactivated links reactivate once both nodes become susceptible.

We use the ER random network model ^[15]. It is one of the most widely used networks for modeling complex processes including population dynamics models. It is a type of homogeneous network, in which each node is assumed to have nearly the same degree. The links are randomly created with a fixed probability such that the average degree distribution $\kappa = 2K/N$. The network topology is fixed with $N$ nodes and $K$ links.

The system of ODEs that approximates the change in the number of nodes and links can be expressed as follows ^[16,17]:

here, $\left[\widehat{XY}\right], \left[\widehat{XY}Z\right], \left[X\widehat{YZ}\right]$ where $X, Y, Z\in \mathrm{S}, \mathrm{I}$ are temporarily deactivated edges. Also, the following relation holds: $\left[\widehat{SI}\right] = \left[\widehat{IS}\right], \left[\widehat{SS}I\right] = \left[I\widehat{SS}\right], \left[\widehat{IS}I\right] = \left[I\widehat{SI}\right]$. The closure relation for deactivated links is as follows:

The system of ODEs can be further simplified by substituting the conservation of nodes and links:

The lower-dimensional approximations suitable for computing control methods are (2.4). The values of $\left[X\right], \left[XY\right], \left[\widehat{XY}\right]$ where $X$, $Y$ ∈ $S, I$ are assumed continuous.

2.2.   Stochastic individual-based model

IBMs simulate the interactions of multiple abstract entities computationally ^[18,19]. A variety of Python packages for running IBM simulations has come into existence over the last few years. In this study, we use Mesa (version 0.8.8.1), which provides an individual-based modeling framework, and NetworkX, (version 2.6.3), which facilitates creating and investigating the dynamics and structures of complex networks.

The IBM generates a random ER network with the same parameters used in the pairwise model. These include the total population, the average degree of individuals, the number of infected, the virus spread chance per SI link at each time step, the recovery rate per individual, the deactivation rate $d$, and the reactivation rate $a$. At each step, deactivated links are activated with the reactivation rate $a$, and SI links are deactivated with the deactivation rate $d$. Then, infected nodes spread the infection to all their susceptible neighbors with a transmission rate ${\rm{ \mathsf{ τ} }}$.

3.   Control

This section describes the process for computing control signals for the pairwise ODEs and the stochastic IBMs. These models are placed in a feedback loop, where control signals are updated at each time step. Individuals are assumed to be willing to maintain their active contacts when there is no threat to infection. Extreme measures, such as deactivating all connections pausing the whole propagation, are considered unrealistic and therefore neglected.

The purpose of applying control measures is to eradicate the epidemic and at the same time keep the number of activated connections as many as possible before the end of the period. We incorporate the time-dependent control, a method that iteratively tunes the control signals at each time step. The deactivation rate $d$ and the reactivation rate $a$ become piecewise constant functions as they can change every step. These rates are bounded as follows:

A fixed time step $\Delta t > 0$ determines how often the control signals are updated. In realistic scenarios, the value of $\Delta t$ translates to the period of updating control policies every $\Delta t$ days or every $\Delta t$ weeks.

As carried out in ^[2], the target numbers to be achieved within the finite time horizon $T > 0$ are denoted as $I\left(T\right) = {I}^{*}$ for infected individuals and $n\left(T\right) = {n}^{*}$ for the average number of active links. We set these values as ${I}^{*} = 0$ and ${n}^{*} = n\left(0\right)$, which translates to zero infected individuals and fully activated social contacts. ${\rm{ \mathsf{ ε} }}$ is defined to measure the effectiveness of the control as follows:

An ideal control scheme would drive the system such that ${\rm{ \mathsf{ ε} }} = 0$.

3.1.   Controlling the pairwise model

Substituting the bounds of the rates into the system yields:

The piecewise constant functions $a$ and $d$ are updated at each step so that these signals are applied to the system such that the target values reach lower or equal to ${\rm{ \mathsf{ ε} }}$.

The NMPC method is an adaptation of the model predictive control method, where sequences of control values are repeatedly computed over the finite prediction range to minimize the objective functional ^[20,21].

To implement the NMPC method, the ODEs are first discretized using the Runge-Kutta-Fehlberg method (RKF45), a method of 4^th order, and an error estimator of 5^th order. Then the system is vectorized to ease matrix computations:

where $x\left(k\right) = X\left(k\Delta t\right), y\left(k\right) = Y\left(k\Delta t\right)$, and $c\left(k\right) = C\left(k\Delta t\right)$. $P$ denotes the length of the prediction horizon, which shortens at each time step. $P$ is used to calculate the sequence of control actions $C$ at each time step $k\Delta t$ until the time horizon is reached, such that ${C}_{i}(k+j)$, where $i = \mathrm{1, 2}$, $j = \mathrm{0, 1}, ...., P-1$. These control signals produce solutions of the system at each time step:

The objective functional $J$ is defined for computing series of control signals $C\left(k\right), C\left(k+1\right), \dots, C\left(k+P-1\right)$ to let the targets $Y\left(k+1\right), Y\left(k+2\right), \dots, Y\left(k+P\right)$ converge to $\left({I}^{*}, {n}^{*}\right)$:

here, ${\lambda }_{1}, {\lambda }_{2}, {\lambda }_{3}, {\lambda }_{4}$ are damping parameters and ${y}_{1}, {y}_{2}, {C}_{1}, {C}_{2}$ are the first and second elements of $Y$ and $C$, respectively. ${\Delta C}_{1}(k+j) = a(k+j)-a(k+j-1)$ and ${\Delta C}_{2}(k+j) = d(k+j)-d(k+j-1)$ measures the cost occurring from changing control signals. ${\lambda }_{1}{\left({y}_{1}\left(k+j\right)-{I}^{*}\right)}^{2}$ accounts for the difference between the current number of infected individuals and the target, and ${\lambda }_{3}{\left({\Delta C}_{1}\left(k+j\right)\right)}^{2}$ accounts for that of the average active degree and the target.

The series of control values $c\left(k\right), c(k+1), ..., c(k+P-1)$ are obtained by computing the optimization equation using the nonlinear least-squares solver lsqnonlin from the MATLAB Optimization toolbox. The first value of the sequence $c\left(k\right)$ is applied to $x\left(k\right)$, updating the system to $x(k+1)$ for the time period $[k\Delta t, (k+1)\Delta t]$. The entire process is repeated until the system reaches the end of the time horizon $T$.

3.2.   Controlling the stochastic individual-based model

Two different cases are considered for computing the control inputs for the IBM ^[2]. In case A, the entire sequences of control signals ${C}_{1}$ and ${C}_{2}$ for $\left[\right(k-1)\Delta t$, $k\Delta t]$ where $k = \mathrm{1, 2}, ...., K$ and $K\Delta t = T$, are computed by the pairwise ODEs. The sequences of control signals are then applied to the IBM at their respective time step over the finite time horizon. In case B, the entire sequences of control signals are computed in the same way as done in case A. However, only the first values in the sequences are applied to the IBM. The same process is repeated until the system reaches the end of the time horizon. We note that the model parameters and the control parameters of the IBM are the same for pairwise ODEs.

4.   Simulations

This section presents the results of the simulations where the stochastic IBM is controlled using cases A and B.

4.1.   Case A

We begin by running extensive simulations on the pairwise ODEs and use the Monte Carlo algorithm to find the maximum value of ${\rm{ \mathsf{ τ} }}$ to two decimal places that drives the solution to desired values $({I}^{*} = 0$ and ${n}^{*} = n\left(0\right))$. The sequence of control signals can then be obtained and applied back to the stochastic IBM. Figure 2(a) presents the effect of those signals on the prevalence and the mean active degree of both the pairwise ODEs and the stochastic IBM. Figure 2(b) shows that the target values of the system of pairwise ODEs diverge such that ${\rm{ \mathsf{ ε} }} > 10\%$ when ${\rm{ \mathsf{ τ} }}$ is higher than 0.15, which is the maximum value of ${\rm{ \mathsf{ τ} }}$ that drives the system such that $\varepsilon \le 10\mathrm{\%}$.

This experiment is repeated with different values of ${\rm{ \mathsf{ τ} }}, {a}_{max}$, and ${d}_{max}$. Table 1 compares the target average active degree and the number of infected nodes at the terminal with those of simulation results. We observe that the number of prevalence and the mean active degree for stochastic simulations converge to the target values when those for pairwise models do so as well, and vice versa. Also, both models tend to be controllable when the transmission rate is not too large.

4.2.   Case B

Figure 3(a) shows that case B is an effective method to control the stochastic IBMs up to $\tau = 0.18$. Figure 4 shows the control signals obtained from the pairwise ODEs in case B. Control signals are piecewise constant functions applied to the stochastic IBMs in Figure 3(a). Figure 3(b) shows the temporal dependence of the infected and the mean active degree of the pairwise ODEs and the stochastic IBM, where the same parameters are used as in Figure 3(a) but with ten times longer time step ∆t.

5.   Discussion

Regarded as one of the simplest disease dynamics models, SIS disease is simulated so that the study can readily be extended to its variants. SIS epidemics involve connections between susceptible and infectious individuals, which are just sufficient for transmitting infection ^[22].

Exact equations encompassing the probability of realistic size networks being in every possible state are beset with fundamental difficulties, resulting in too many equations for practical applications. For example, there are ${2}^{N}$ distinct state spaces for N individuals for SIS epidemic models.

As done by one of the earliest developed dynamic network-based epidemic models, we use the pairwise modeling framework to relax the high dimensionality. Mean-field models that approximate stochastic processes can be formulated using pairwise approximations, having the potential to be exact in large systems. However, its focus on expected population-scale quantities makes it difficult to analyze variations of epidemic dynamics in different network structures ^[10,13].

Modern IBMs simulating simultaneous interactions of multiple entities on explicit structures have shown increasingly superior to traditional compartmental models in regard to accuracy in predictions ^[18,19]. On the other hand, an increase in uncertainty is inevitable when parameters are added for building more realistic IBMs.

Computing control signals of the stochastic IBMs using the system of pairwise ODEs requires much less computational cost than computing control signals directly from the stochastic model. That is, the rates at which the SI links temporarily deactivate and SS links reactivate are regulated by the system of pairwise ODEs. As a result, we examine how well an actual spread can be controlled using the pairwise ODEs.

Notably, the outcomes are contrasting when the method used in case A is repeated on dynamic networks with contact-conserving rewiring mechanisms. In ^[2], case A is deemed ineffective since control signals obtained from the controllable system of pairwise ODEs do not drive the stochastic model to the target values, even for lower values of $\tau$. The effectiveness of the methods used in case A may highlight the variation in responses between different dynamic network models to model predictive control methods. On the other hand, the difference may be simply because our model considers a larger population than that used in ^[2]. We recall that the sensitivity of the network-based simulation outcomes to the property of the networks is difficult to assess ^[10].

Our outcome agrees with that of ^[2] in that the process used in case B is a more effective approach than that of case A when using the NMPC method to control the stochastic simulations of dynamic networks. Also, when the time step ∆t is increased as in ^[2], the probability of stochastic simulations being effectively controlled decreases.

In Figure 4, the deactivation rate reaches the maximum at earlier stages of the process. Link activation rate remains low at this stage so that the number of prevalence can decrease. However, the link activation rate increases steadily so that the mean active degree can reach the target value. These observations can be seen as control measures focusing on the early eradication of diseases and steady relaxation of social distancing. Seeking an optimal balance between minimizing the number of infected and maintaining active social contacts is intentionally planned through the damping parameters of the model predictive control.

We highlight that the controllability of our simulations highly depends on the control parameters, including the time horizon, the step size, epidemic parameters, and damping parameters that significantly impact the optimization processes of the model predictive control.

6.   Conclusions

A limited amount of research exists where pairwise ODEs are used to compute control signals for controlling real networks modeled with IBMs. This paper considers two different cases for computing control inputs for dynamic networks with a temporary link mechanism. The ER network is assumed. In the first case, little computational resources are used to compute the control signals. That is, control inputs are calculated only once at the beginning of the control period. In the second case, we collect the epidemic data of the stochastic IBM at each time step and use the data as the initial condition of the pairwise model for computing the control signals for stochastic IBM. Our results show that the prevalence and the mean active degree converging to the target values increased when control signals were updated more frequently.

The goal of the study is to explore the opportunity arising from linking the modern network-based epidemic models and control, thereby gaining an intuitive understanding. The scope is limited to dynamic networks with a temporary link deactivation mechanism, investigating the degree of effectiveness for controlling stochastic IBMs using the results obtained from the pairwise ODEs. We compute control signals from ODE approximations, which requires fewer costs and efforts than calculating those from stochastic simulations. Our method provides a stepping stone towards the applications of controlling real-world complex systems.

When the transmission rate is low, the effectiveness of control measures predicted by pairwise ODEs coincides with stochastic simulations applied with the same control signals used for the pairwise model. On the other hand, when the transmission rate is high, case B outperforms case A for controlling stochastic simulations using pairwise ODEs. As with the limitations of all network-based simulations, where the sensitivity of the epidemiological outcomes to the specifications of the network is unclear, we highlight that the controllability strongly depends on the parameters that define the models ^[10].

Our approach can be extended to more refined settings, complementing the simplifications and limitations. While our approach concentrates on the end target of the prevalence and active degree, studies can look into minimizing the peak or the total number of infections. Such an approach may be more practical in the context of epidemiology when maintaining the intensive care units (ICUs) capacity rises as an issue. Moreover, studies can account for the costs occurring from the activation and deactivation of links. They do not focus on network control, and this cannot be examined with simple compartmental models. Finally, we used the homogeneous mean-field pairwise approximations, which do not fully capture the system behavior of networks in which degree distributions are heterogeneous. Since many observed networks are far from homogeneous, the next step shall be incorporating degree heterogeneity, significantly affecting disease dynamics ^[10].

Acknowledgments

This research was a part of the project titled “Understanding the deepsea biosphere on seafloor hydrothermal vents in the Indian Ridge” (No. 20170411), funded by the Ministry of Oceans and Fisheries, Republic of Korea. The initial stage of the work was also supported by the research programs of the Korea Institute of Ocean Science and Technology as part of the Shared Use Program (No. PE99656 and PE99796) and the Fellowship for Fundamental Academic Fields by Seoul National University (No. 2020-86).

Conflict of interest

All authors declare no conflicts of interest in this paper.