Microscopic Numerical Simulations of Epidemic Models on Networks

Abstract

Mathematical models of the spread of epidemic diseases are studied, paying special attention to networks. We treat the Susceptible-Infected-Recovered (SIR) model and the Susceptible-Exposed-Infectious-Recovered (SEIR) model described by differential equations. We perform microscopic numerical simulations for corresponding epidemic models on networks. Comparing a random network and a scale-free network for the spread of the infection, we emphasize the role of hubs in a scale-free network. We also present a simple derivation of the exact solution of the SIR model.

1. Introduction

The Coronavirus Disease 2019 (COVID-19) [1], first identified in Wuhan, China on December 2019 [2,3], has spread globally. The World Health Organization (WHO) declared the outbreak as a public health emergency of international concern on 30 January 2020. Furthermore, the WHO announced the COVID-19 outbreak as a ’pandemic’ on 11 March 2020 [4,5]. Various epidemiological models are used to predict and analyze the spread of COVID-19. Most commonly, the epidemiological model is a compartmental model in which each individual is considered to be in one of the possible types: for example, susceptible (S), infected (I), or recovered (R). The proportions of individuals in each type are regarded as continuous variables, and their rate equations (time derivatives) are set up as a function of these proportions with appropriate rate constants. By solving these differential equations, the fraction of each type as a function of time is obtained. The Susceptible-Infected-Recovered (SIR) model, proposed in 1927 by Kermack and McKendrick [6], is a basic model of such a framework.

The SIR model [6,7] has been used as a mathematical model of epidemics in the case of COVID-19 [8,9]. Some extended models have also been developed, such as the Susceptible-Infectious-Recovered-Deceased (SIRD) model [10,11,12]. This model differentiates between recovered and deceased individuals. Another example is the Susceptible-Exposed-Infectious-Recovered (SEIR) model [13,14,15,16], where an ’exposed’ category is added. It considers the situation involving a significant incubation period, during which individuals have been infected but are not yet infectious.

It is obvious that a given infective individual does not have an equal probability of infecting all others. Each individual only has contact with a small fraction of the total population, and the number of contacts that people have can vary greatly from one person to another. The connection between individuals can be described by a network. The nodes (vertices) in the network are connected by edges, which represent the interaction between those nodes. The number of edges of a node is referred to as the degree of that node. In the context of a network, the existence of ’hub’ is important. A hub is a node that has connections with many other nodes. Many biological interactions are not homogeneous and not random [17,18,19,20]. The degree distribution of nodes in a biological network can often be described by a power-law [21]. Networks with this kind of property are known as scale-free.

The importance of the network structure in the analysis of epidemics was pointed out, and a number of studies have been done [22,23,24,25,26,27,28]. In connection with COVID-19, stochastic simulations of epidemic models on a scale-free network were studied aiming at showing effective mitigation strategies [29]. The spread of COVID-19 under local quarantine measures (K-quarantine) in South Korea was modeled using the SEIR model on complex networks [30].

In the present paper, we study the mathematical models of epidemics, paying special attention to networks. We perform microscopic numerical simulations of the SIR and SEIR models on random and scale-free networks. We compare the results of the microscopic simulations with those of the macroscopic model described by differential equations using the same parameters, rate constants. We elucidate the role of hubs in a scale-free network and discuss the effective method of mitigation strategy. We also show a simple derivation of the analytical solution of the SIR model, and the computer code in Python for the analytical solution will be presented. The present authors previously published a tutorial for students to describe a mathematical model of epidemics, mainly focusing on differential equations [31]. This article is a continuation of the tutorial, so that the description is planned to be understandable for students.

This paper is organized as follows. In Section 2, the outline of network theory is described. The Erdös-Rényi network [32,33], a random network, and the Barabási-Albert network [21], a scale-free network, are introduced. In Section 3, we explain the SIR model of differential equations and present a simple derivation of the exact solution of the SIR model. We perform the microscopic simulation of the SIR model on networks. In Section 4, the SEIR model is introduced, and the microscopic simulations on networks is performed. The final section is devoted to a summary and discussion. The Python code of the exact solution of the SIR model is given in Appendix A.

2. Network

2.1. Outline of Network Theory

Network science is an academic field which studies complex networks, and it has applications in many disciplines including statistical physics, computer science, electrical engineering, biology, economics, finance, public health, etc. [34,35,36]. It covers a wide range of topics from social networking services (SNS) to movie-actor networks.

Networks can be used to represent connections between individuals. The nodes in the network are connected by edges, which represent the relation between those nodes. Two nodes sharing the same edge are considered to directly interact with one another and are referred to as neighboring nodes. The number of edges of a node is referred to as the degree of that node. Nodes with many edges that greatly exceed the average are referred to as `hubs’. In computer science and network science, network theory is a part of graph theory: a network can be defined as a graph in which nodes and/or edges have attributes. A network can be an exceedingly complex structure, as the connections among the nodes can exhibit complicated patterns.

The people who study the compartmental model of differential equations are not necessarily familiar with the network theory. To fill in the gaps, the history of network theory will be surveyed. Most of the descriptions are following a book by Caldarelli and Catanzaro [36].

Graph theory is considered to begin when, in 1735, Euler [37] solved a popular puzzle about bridges [38]. Euler started a whole branch of mathematics, built on graph analysis. His intuition can be considered as the first, foundational moment of network science. After him, many mathematicians studied the formal properties of networks, while scientists applied them to a wide range of problems [39,40,41]. Moreno [42] was one of the first researchers to apply the idea of network to society. His work started one of the most important lines of network science: the analysis of social networks.

The next important moment in the foundation of network science came with Erdös and Rényi [32,33]. They studied a mathematical model representing a graph where vertices are connected to each other completely at random. This model, lately known as random graph, is a simplified model, and its properties are very different from those of real networks. However, the random graph model is very important. It can be used in the probabilistic method to prove the existence of graphs satisfying various properties or to provide a rigorous definition of what it means for a property to hold for almost all graphs. The simplest procedure for building a random graph is the following. We take all the possible pairs of vertices. For each pair, we draw a link or pass to the next pair with some probability p, until all the pairs are finished.

Real-world networks are not homogeneous and not random. In 1967, Milgram and Travers [43,44] undertook a set of experiments examining the average path length for social networks of people. The research suggested that human society is a small-world-type network characterized by short path-lengths. A chain of “a friend of a friend” statements can be made to connect any two people in a maximum of six steps, which is called the notion of “six degrees of separation”. The small-world property is something intrinsic to networks. By contrast, regular grids or lattices do not display it. In 1998, Watts and Strogatz [45] tried to produce graphs with small-world properties. They started by considering a very simple, regular structure. As the initial structure, they construct a regular ring lattice, a graph with N nodes each connected to k neighbors. They cut one of the links in the initial structure and rewired it with another node, chosen at random. A small number of shortcuts is enough to bring all the elements of the system much closer to each other. The key ingredient that transforms a structure of connections into a small world is the presence of a little disorder.

The appearance of hubs, heterogeneity, scale-invariance, and  fat-tailed distributions are observed in a set of different networks. The major limitation of the Watts–Strogatz network is that it produces an unrealistic degree distribution, that is, it cannot produce scale-free properties. In contrast, Barabasi and Albert [21] put forward a mathematical model of the growth of a network, having hubs and scale-free degree distribution. They imagined a graph that starts with a small set of vertices, connected at random. New nodes are added at a steady rate to this initial nucleus, each of them carrying a given number of links. A simple rule establishes how new nodes are linked: incoming vertices prefer old ones that already have many links. This mechanism is called “preferential attachment”.

In 2001, Pastor-Satorras and Vespignani [22] studied the problem by modeling and simulating the spread of a disease in a social network. They treated a minimal model of epidemic, Susceptible-Infected-Susceptible (SIS) model. They found that the epidemic threshold depends crucially on the features of the underlying network: a random network or a scale-free network.

Real-world networks are scale-free? It has been a problem. In 2019, Broido and Clauset [46] claimed that a strongly scale-free structure is empirically rare, while for most networks, log-normal distributions fit the data as well or better than power laws. Even if the distribution does not follow a pure power law, the importance of hubs is obvious.

2.2. Network Generation

Here, we treat two types of networks; random networks and scale-free networks. For random networks, each link is added between two randomly selected nodes. Thus, the degrees of each node have a Poisson distribution. Because  this model was first proposed by Erdös-Rényi, it is often called the ER network [32,33]. The degrees of each node are heterogeneous for scale-free networks, following the power law. This implies that a few nodes have large degrees, but the remaining nodes have small ones. The nodes with many neighbors are called hubs. Because this model was first proposed by Barabási-Albert, it is often called the BA network [21]. The BA network is built with an iterative growth based on preferential attachment.

We use igraph (Version 0.8.2) package of Python (Version 3.7.3) to generate networks.

2.2.1. Erdös-Rényi Network

An example of the ER network, a random network, and its histogram of the distribution of degrees are shown in Figure 1. The number of nodes is $N=80$ and the average number of degrees is $\langle k\rangle =4$. The command to generate the graph is igraph.Graph.Erdos_Renyi().

2.2.2. Barabási-Albert Network

Figure 2 shows an example of the BA network, a scale-free network, and its histogram of the distribution of degrees. The number of nodes is $N=80$ and the average number of degrees is $\langle k\rangle =4$, which is the same as the ER network. We see from Figure 2b that there are nodes with a large number of many degrees. The command to generate the graph is igraph.Graph.Barabasi(). Because the BA network is an evolving network, we use blue colors for the nodes of the first three steps in Figure 2a.

To check the scale-free property of the BA network, we show the log-log plot of the degree distribution for a large BA network in Figure 3. The number of the nodes is $N=$ 10,000, and the average connecting number is $\langle k\rangle =8$. An average over 40 samples was taken. This sample of BA networks will be used in the simulation of microscopic epidemic models on the network later. For convenience, we plot $k^{-3}$ by red dotted curve. We observe a power law of the form, $P\left(k\right)\propto k^{-3}$, for the degree distribution resulting from the BA network. This power dependence was derived by the mean-field theory [21,47] and the master equation method [48].

3. SIR Model

3.1. Differential Equation of SIR Model

We first explain the SIR model. (See Ref. [31].) We consider a closed society of N individuals and classify individuals into three types: susceptible (S), infected (I), and recovered (R). Infected individuals can only transmit the virus to susceptible individuals. Once infected individuals have recovered (or have passed away), they can no longer infect others and cannot be reinfected. The spread of epidemic diseases is described by the changes in the numbers of the three types of individuals. We denote the number of susceptible (S) individuals by $x\left(t\right)$, the number of infected (I) individuals by $y\left(t\right)$, and the number of recovered or removed (R) individuals by $z\left(t\right)$, as a function of time t. We do not consider the birth and death processes. We illustrate the scheme of the SIR model in Figure 4.

The SIR model is expressed as simultaneous differential equations of variables $x\left(t\right)$, $y\left(t\right)$, and  $z\left(t\right)$, such as

$$\begin{array}{ccc}\hfill \frac{dx\left(t\right)}{dt}& =& -\frac{\beta }{N}x\left(t\right)y\left(t\right),\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \frac{dy\left(t\right)}{dt}& =& \frac{\beta }{N}x\left(t\right)y\left(t\right)-\gamma y\left(t\right)=\left(\frac{\beta }{N}x\left(t\right)-\gamma \right)y\left(t\right),\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \frac{dz\left(t\right)}{dt}& =& \gamma y\left(t\right),\hfill \end{array}$$

(3)

where $\beta$ is the rate of infection and $\gamma$ the rate of recovery (the rate of quarantine). We consider that each person is in contact with k persons per unit time (day), and the probability of infection for each contact is set as p. Then, $\beta$ is given as $kp$. We assume that the total number of individuals is set to be constant such that

$$x\left(t\right)+y\left(t\right)+z\left(t\right)=N.$$

(4)

The SIR differential Equations (1)–(3) can be solved numerically by using a standard Runge–Kutta method or a Gillespie algorithm of event-driven type. The exact analytical solution of the differential equations was presented by Harko et al. in 2014 [49]. They demonstrated that these simultaneous equations can be expressed as a differential equation of a single variable $x\left(t\right)$. Then, the problem is reduced to a Bernoulli differential equation, which can be solved exactly. However, there was a detour in the derivation of the exact solution in Ref. [49]; they also used $z\left(t\right)$ in the process of derivation. Only two variables $x\left(t\right)$ and $y\left(t\right)$ are independent because of Equation (4), and  $x\left(t\right)$ and $y\left(t\right)$ form closed Equations (1) and (2). Okabe and Shudo presented a more straightforward derivation using only $x\left(t\right)$ and $y\left(t\right)$ in a tutorial paper for epidemics [31].

The SIR model is treated as a typical system of nonlinear differential equations and has been studied from various mathematical points of view. The nullclines analysis was described in a textbook by Hirsch, Smale, and Devaney [50]. There, the flow in the $xy$(SI)-plane is discussed.

In the present paper, we present a simple derivation of the exact analytical solution of the SIR model referring to the textbook by Hirsch, Smale, and Devaney [50]. We also give a Python code for the exact analytical solution of the SIR model, for convenience.

3.2. Exact Analytical Solution

We here show a simple derivation of the analytic solution of the SIR model. From Equations (1) and (2), we obtain

$$\frac{dy}{dx}=\frac{-(\beta /N)x+\gamma }{(\beta /N)x}=-1+\frac{\gamma }{\beta }\left(\frac{N}{x}\right).$$

(5)

The integration with respect to x yields

$$y=-x+\frac{\gamma }{\beta }Nlnx+C,$$

(6)

where C is an integral constant. In the textbook by Hirsch, Smale, and Devaney [50], this integral form was used in the discussion of flow in the $xy$(SI)-plane.

As initial conditions, we consider

$$\begin{array}{ccc}\hfill x\left(0\right)& =& N_1,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill y\left(0\right)& =& N-N_1,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill z\left(0\right)& =& 0.\hfill \end{array}$$

(9)

Using these initial conditions, we have

$$C=N-\frac{\gamma }{\beta }Nln{N}_1.$$

(10)

Then, $y\left(t\right)$ is expressed as

$$y\left(t\right)=-x\left(t\right)+\frac{\gamma }{\beta }Nln(x\left(t\right)/N_1)+N.$$

(11)

Substituting Equation (11) into Equation (1), we obtain

$$dt=\frac{dx}{(\beta /N)x[x-(\gamma /\beta )Nln(x/N_1)-N]}.$$

(12)

By integrating this equation, t is obtained as a function of $x\left(t\right)$ as follows:

$$t={\int }_1^{x\left(t\right)/N_1}\frac{du}{\beta u[(N_1/N)u-(\gamma /\beta )lnu-1]},$$

(13)

where we changed a variable such that $u=x/N_1$. This form is essentially the same as Equation (26) of Ref. [49] and Equation (52) of Ref. [31]. However, we do not use the Bernoulli equation. We have obtained the exact analytical form simply and directly. When $x\left(t\right)$ is obtained, $y\left(t\right)$ is calculated using Equation (11). Moreover, $z\left(t\right)$ is calculated through

$$z\left(t\right)=N-x\left(t\right)-y\left(t\right).$$

(14)

We note that for rate constants only the ratio of $\beta$ and $\gamma$ appears. The number $R_0$ defined by

$$R_0=\frac{\beta }{\gamma }$$

(15)

is known as the basic reproduction number [51,52,53], and the number of infected individuals increases when $R_0>1$, whereas it decreases when $R_0<1$. It is also noteworthy that by introducing scaled variables $\tilde{x}=x/N,\tilde{y}=y/N,\tilde{z}=z/N$, and  $\tilde{t}=\beta t$, we can rewrite the Equations (1)–(3) in scaled forms.

The Python program is given in the Appendix A. Numerical integration of Equation (13) is performed by summing up an integrand of the right-hand side with a sufficiently small step size. Then, t is calculated as a function of $x\left(t\right)$. For numerical values, we set $N=$ 10,000 and $N_1=9900$. The calculated results using parameters $\beta$ = 0.4 (solid), 0.45 (dashed), 0.5 (dotted) are shown in Figure 5; $\gamma$ is fixed as 0.2.

Because we have the Python code for the exact solution, we can demonstrate various quantities easily. The directions of flow in the $xy$-plane are shown in Figure 6. The parameters are the same as in Figure 5. Every 5 in time (t), the directions of flow are depicted by arrows. The length of the arrow is measured for the period of 1.2 in t.

We here make a comment on the final size equation. Putting $y(\infty )=0$ in Equation (11), we have the relation for the final value of $x(\infty )$:

$$x(\infty )=\frac{\gamma }{\beta }Nln(x(\infty )/N_1)+N.$$

(16)

In the limit of $N_1\to N$, we obtain

$$x(\infty )/N-1=\frac{\gamma }{\beta }ln(x(\infty )/N).$$

(17)

Then, we have

$$1-z(\infty )/N=exp\left[-\frac{\beta }{\gamma }(z(\infty )/N)\right].$$

(18)

This relation is known as the final size equation [54]. This derivation of the final size equation through Equation (6) was discussed by Miller [55]. For $R_0=\beta /\gamma =2.0\text{,} 2.25$ and 2.5, we have $z(\infty )/N\to$ 0.797, 0.853, and 0.893, respectively.

3.3. Method of Microscopic Simulation of SIR Model

We perform a microscopic simulation of the SIR model on a network. We follow a procedure, of which parameters correspond to those for the differential equation of the macroscopic variables. The actual procedure, which is similar to that by Herrmann [29], is as follows:

• Generate a network.
• At $t=0$, an individual or individuals are infected (I).
• A susceptible individual (S) will be infected (I) with a probability p if a connecting individual (one of k) is infected (I). In terms of the SIR model, the parameter $\beta$ is $\beta =kp$.
• An infected individual (I) will be recovered (R) in $1/\gamma$ days on average. The infected period is chosen by a Poisson distribution with the average of $1/\gamma$.
• At each time t, the processes 3, 4 are repeated.
• The time sequence obtained from the above procedure is regarded as a single sample. Simulations are performed for several samples.

3.4. Results of SIR Model on the Erdös-Rényi Network

We show the simulational results of the microscopic SIR model on the ER network. First, we treat the case that only a single individual is infected at $t=0$. The average number of the connecting nodes per a node is $\langle k\rangle =8$, and the probability p of infection is 1/20, which leads to $\beta =\langle k\rangle p=0.4$. The average infected period, $1/\gamma$, is chosen as 5 days, which is realized by a Poisson distribution with the average value of 5. The total number of nodes (individuals) is 10,000. The temporal variations of the number of individuals for three types (S, I, R) are plotted in Figure 7. We simulated 40 samples. All the results of $x\left(t\right),y\left(t\right),z\left(t\right)$ for 40 samples are plotted in Figure 7a. There are several sources of randomness in the temporal variation of each sample; the generation of a network, the choice of the node of an initial infected individual, the probability p of transmission, and  the choice of the infected period $1/\gamma$ with a Poisson distribution. The temporal evolution is classified to two types: (i) The number of infected individuals (I) increases with time, reaches a peak, and gradually decreases, which is the same behavior as that for the SIR model of the differential equation, as shown in Figure 5. (ii) The infection vanishes quickly and does not spread throughout the network; the number of susceptible individuals (green curve) remains almost the same. The behavior of the second type was observed in the simulation by Herrmann and Schwartz [29] and  is regarded as the absorbing state [56,57] in the contact process [58]. The transition to the absorbing state is related to the percolation phenomena [59]. For the parameters of Figure 7, about 14% of samples followed the time variation of the second type. Starting from a single infected individual, a disease is transmitted to another individual with a probability of $\beta =0.4$. However, some infected individuals disappear shortly in the infected period with the Poisson distribution. Then, the infection stops with some probability.

The averaged values of the time evolution of 40 samples are plotted in Figure 7b. Because there are two types of time evolution, the average proportion of the final infected individuals, $z(\infty )/N$, becomes smaller than that of the macroscopic SIR model, 0.797 for $R_0=2.0$.

In order to compare with the macroscopic SIR model, we perform other simulations. That is, 10 individuals are set to be infected initially at $t=0$. Then, we did not observe the type 2 behavior, absorbing state. Even if some initially infected individuals are caught in an absorbing state, others transmit infection and the disease will spread. We simulate the same situation as the macroscopic SIR model, where there is no absorbing state for continuous variables. In Figure 8, we plot the time variation of the macroscopic SIR model with initial 10 infected individuals. Initially infected individuals are chosen randomly. In Figure 8a, we plot all the results of 40 samples, whereas the average over them is plotted in Figure 8b. The microscopic model on the random network can be considered to reproduce the behavior of the macroscopic SIR model. The solid curves in Figure 5 are the exact solution of the SIR model of differential equations with the same rate constants as Figure 8; $\beta =0.4$ and $1/\gamma =5.0$. The time of the peak for the infected individuals, the peak value, and the final number of recovered individuals ( $z(\infty )$) or the individuals who experienced infection, take almost the same values as those of the macroscopic SIR model, but the quantitative coincidence is not complete.

Because we mainly focus on the comparison of the microscopic simulation with the epidemic model of differential equations, in the rest of the paper we do not deal with the cases of the absorbing states.

3.5. Results of SIR Model on the Barabási-Albert Network

We plot the simulational results of the microscopic SIR model on the BA network, a scale-free network, in Figure 9. We set the number of the nodes N and the average number of the connecting nodes $\langle k\rangle$ as 10,000 and 8, respectively, which are the same as the case of the ER network. Because there is a hub structure in scale-free networks, the spread of infection may depend on the origin of the initially infected individuals. We consider two cases. One is the case that initially 10 infected individuals are within the hub, and the other is the case that initially 10 infected individuals are outside the hub. We define the hub by the node which has connecting nodes of $k\ge 14$. Figure 9a,b are the results within the hub and outside the hub, respectively. An average over 40 samples was taken for both cases. First, we compare the behavior of the BA network with that of the ER network. We observe the increase of the infected individuals, reaching a peak, and the gradual decrease. The burst of the infection is rapid compared with the case of the ER network. The peak is higher. The decrease of the infection is also fast. If we compare the final number of the recovered, $z(\infty )$, this number of the BA network is smaller than that of the ER network. The average number of the connecting nodes, $\langle k\rangle$, is the same for the BA network and the ER network. The existence of the hub for the BA network, a scale-free network, stimulates the rapid increase of the infection. At the same time, the existence of a large amount of nodes of very small connecting nodes leads to the smaller final size of $z(\infty )$. For the difference of initially infected individuals, if the initially infected individuals are within the hub, the increase of the infection is more rapid and the peak of the infected individuals is higher than the case of the initial individuals outside the hub. However, the difference is not so large.

To elucidate the role of the hub in the BA network, we investigate the contribution of the individuals within the hub for the time variation of the infection. We plot the time sequence up to $t=40$ in Figure 10, and the rate constants are the same as in Figure 9. The solid curves in Figure 10 represent the time variation of the total individuals, whereas the dashed curves represent that of the individuals within the hub (10%). Initial 10 infected individuals are chosen randomly. The results of the average over 40 samples are shown. The hub is defined by the node which has connecting nodes of $k\ge 14$. For the convenience of comparison, the ten times of the values of the hub are shown in dotted curves. We observe that the number of infected individuals within the hub shows significantly rapid growth. They reach the maximum very fast, and the peak value is high. The decrease speed is also fast, and almost all the individuals within the hub will be infected. That is, there are almost no noninfected individuals ( $x(\infty )=0$). On the other hand, the nodes outside the hub grow rather slowly, and the decrease is also slow. Some amount of nodes remain noninfected.

We examine the role of the hub from a different point of view. We plot the distribution of degrees of the final noninfected individuals $x(\infty )$ for the BA network in Figure 11. The parameters of the rate constants are the same as before. The initial 10 individuals were chosen randomly, and the average over the 40 samples was taken. The gray bar represents the distribution of total nodes. For the BA network, the distribution of small k is large, and it has a tail with power dependence for large k; the nodes of the hub are shown by an arrow. The log-log plot of the degree distribution was shown in Figure 3. The distribution for the final noninfected individuals was represented by green bars. The absolute value of the final noninfected individuals is larger for the BA network compared to the ER network. (See Figure 8 and Figure 9) We observe in Figure 11 that the proportion of nodes of small k is especially large for the final noninfected individuals. That is, before the infection reaches the nodes of small connections, the overall spread of the infection ceases.

3.6. Mitigation Strategy

Next, we investigate the mitigation strategy of infection for the BA network. We focus on the role of the hub. Various mitigation strategies have been taken to slow down the rapid spread of COVID-19, such as complete or partial lockdown, travel bans, mass gathering restrictions, home quarantines within communities, social distancing measures, etc. Here, as a model of mitigation, we consider the isolation of individuals by restricting the contact.

Let us consider the trial that we quarantine (isolate) 10% of individuals for the time $t\ge 10$. We do not perform the process 3 of the method of simulation for quarantined individuals. We compare the results of two cases in Figure 12. In Figure 12a, we quarantine randomly chosen 10% individuals for $t\ge 10$, whereas in Figure 12b, 10% individuals within the hub are quarantined for $t\ge 10$. Initial 10 infected individuals are chosen randomly, and an average over 40 samples was taken for both cases. If we compare the case of the nonquarantine, shown in Figure 9, the quarantine curbs the spread of the infection, and the final number of recovered individuals ( $z(\infty )$), or the individuals who experienced infection, becomes smaller. The effects of quarantine are predominant for the 10% individuals within the hub. The number of the final recovered individuals becomes drastically reduced.

4. SEIR Model

4.1. Differential Equation of SEIR Model

For many important infections, there is a significant latent (incubation) period during which individuals have been infected but are not yet infectious themselves. In the Susceptible-Exposed-Infectious-Recovered (SEIR) model [13,16], an Exposed (E) category for these individuals is added to the SIR model.

We denote the number of susceptible (S) individuals by $x\left(t\right)$, the number of infectious (I) confirmed individuals by $y\left(t\right)$, the number of recovered (R) individuals by $z\left(t\right)$, and the number of exposed (E) individuals by $w\left(t\right)$, as a function of time t. The spread of infection is then described by the following system of nonlinear ordinary differential equations:

$$\begin{array}{cc}\hfill \frac{dx\left(t\right)}{dt}& =-\frac{\beta }{N}x\left(t\right)y\left(t\right)\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \frac{dw\left(t\right)}{dt}& =\frac{\beta }{N}x\left(t\right)y\left(t\right)-\sigma w\left(t\right)\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \frac{dy\left(t\right)}{dt}& =\sigma w\left(t\right)-\gamma y\left(t\right)\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \frac{dz\left(t\right)}{dt}& =\gamma y\left(t\right).\hfill \end{array}$$

(22)

Here, contact rate $\beta$ controls the rate of spread, which represents the probability of transmitting disease between a susceptible and an infectious individual. Incubation rate $\sigma$ is the rate of latent individuals becoming infectious (i.e., the average incubation period is ${\sigma }^{-1}$). Recovery rate $\gamma$ is the rate of infected individuals becoming recovered (i.e., the average recovery period is ${\gamma }^{-1}$). We do not consider birth and death processes. The total number of individuals is constant such that

$$x\left(t\right)+w\left(t\right)+y\left(t\right)+z\left(t\right)=N.$$

(23)

We illustrate the scheme of the SEIR model in Figure 13.    

For the differential equations of the SEIR model, we cannot obtain the solution by employing the transformation applied to the SIR model [49]. Thus, we show the numerical solution using the second-order Runge–Kutta method in Figure 14. The parameters were chosen as $\beta =0.5,1/\sigma =3.0$ and $1/\gamma =5.0$. The number of individuals (nodes) is $N=$ 10,000, and the initial conditions are $x\left(0\right)=N_1=9990,w\left(0\right)=N-N_1=10$, and  $y\left(0\right)=z\left(0\right)=0$.

The ratio of the final infected individuals, $\zeta =z(\infty )/N$, is obtained as the final size equation. In the limit of $N_1/N\to 1$, the final size equation becomes $R_0=-(ln(1-\zeta ))/\zeta$, which is the same form as the SIR model. For $R_0=\beta /\gamma =2.5$, it can be shown that $z(\infty )/N\to 0.893$, which is reproduced in the numerical solution shown in Figure 14.

We briefly make a remark on the final size equation for the SEIR model. Considering the sum of Equations (20) and (21), we obtain a relation, similar to Equation (5) for the SIR model,

$$\frac{d(y+w)}{dx}=\frac{-(\beta /N)x+\gamma }{(\beta /N)x}=-1+\frac{\gamma }{\beta }\left(\frac{N}{x}\right).$$

(24)

Then, the integration with respect to x yields

$$y+w=-x+\frac{\gamma }{\beta }Nlnx+C,$$

(25)

where C is an integral constant. Thus, by replacing y in the SIR model by $y+w$, we have the final size equation for the SEIR model in the limit $N_1/N\to 1$ as

$$\frac{\beta }{\gamma }=\frac{ln\left(x\right(\infty )/N)}{x(\infty )/N-1}.$$

(26)

Or,

$$1-z(\infty )/N=exp\left[-\frac{\beta }{\gamma }(z(\infty )/N)\right],$$

(27)

which is the same form as the case of the SIR model. It is noteworthy that only $\beta /\gamma$ appears and $\sigma$ not. The analytical solution of the SIR model is obtained by inserting $y\left(t\right)$, Equation (6), to the $dx\left(t\right)/dt$ in the differential equation, which results in the closed equation of only $x\left(t\right)$. However, in the case of the SEIR model, the left-hand side of Equation (25) is $y\left(t\right)+w\left(t\right)$, and we cannot insert $y\left(t\right)$ in $dx\left(t\right)/dt$. Thus, we cannot obtain the exact solution of the SEIR model.

4.2. Method of Microscopic Simulation of SEIR Model

We can extend a microscopic simulation on a network to the case of the SEIR model. The actual procedure is as follows:

• Generate a network.
• At $t=0$, an individual or individuals are infected as the latent state (E).
• A susceptible individual (S) will be infected in the latent state (E) with a probability p if a connecting individual (one of k) is symptomatically infected (I). In terms of the SEIR model, the parameter $\beta$ is $\beta =kp$.
• An exposed individual in the latent state (E) will be symptomatically infected (I) in $1/\sigma$ days on average. The latent period is chosen by a Poisson distribution with an average of $1/\sigma$.
• A symptomatically infected individual (I) will be recovered (R) in $1/\gamma$ days on average. The infected period is chosen by a Poisson distribution with an average of $1/\gamma$.
• At each time t, processes 3–5 are repeated.
• The time sequence obtained from the above procedure is regarded as a single sample. Simulations are performed for several samples.

4.3. Results of SEIR Model on the Erdös-Rényi Network

As an example of network, we performed a simulation of the microscopic SEIR model on the ER network. In Figure 15, the temporal variations of the number of individuals for four types (S, E, I, R) are plotted. The average number of the connecting nodes per a node is $\langle k\rangle =8$, and the probability p in the process 3 is 1/16. The parameters $\beta$, $\gamma$, $\sigma$ are chosen to be equivalent to the case of the numerical solution of the macroscopic SEIR model shown in Figure 14. The total number of nodes (individuals) is 10,000. Initial 10 infected individuals in the latent state are chosen randomly. All the results of 40 samples are plotted in Figure 15a, whereas the averaged values of 40 samples are given in Figure 15b.

The obtained results coincide with the results of the differential equations shown in Figure 14. We can perform the microscopic simulation of the SEIR model on other networks.

5. Summary and Discussion

We described the mathematical models of the spread of epidemic diseases, paying special attention to networks. We gave a survey of network theory with an emphasis on the relation to epidemics. Two types of networks were introduced. One is the Erdös Rényi network, a random network, and  the other is the Barabási Albert network, a scale-free network.

A simple derivation of the exact analytical solution of the SIR model [49] was presented. The Python program of the exact solution was given in Appendix A. We do not have to use numerical methods, such as a Runge–Kutta method or a Gillespie algorithm, to obtain the solution of the differential equations.

We systematically studied the microscopic numerical simulations of the epidemic models on networks and the corresponding macroscopic epidemic models described by differential equations. We compared the results of the microscopic and macroscopic models with the same (equivalent) parameters (rate constants). In the case of the ER network, a random network, we obtained the results of the microscopic simulation of the SIR model which reproduce the solution of the SIR model of differential equations. For the BA network, a scale-free network, we observed a rapid burst of infection and a high peak of the number of infected individuals. However, the decrease of the infection was also fast, and the final number of infected individuals is smaller than the case of the random network. All the findings in the scale-free network were interpreted by the effects of hubs. We elucidated the role of hubs in the scale-free network from various points of view. We also investigated the mitigation strategy for the BA network and emphasized that the isolation of individuals in the hub is important to curb the spread of infection. We also studied the macroscopic and microscopic models of the SEIR model.

The present observations are natural consequences from the hub structure of a network. We explicitly showed these results by a systematic comparison of different models. We may learn lessons for a real-world strategy to fight against the pandemic. Because the infection originates from contact, the social network structure, especially the existence of a hub, is important. A quick response to an initial outbreak is more crucial for the hub structure. To develop control strategies, such as the shutdown of individual connections, and the complete or partial isolation of cities, the hub structure should be taken into account. When discussing the cases per 100,000 population, the possibility of secondary infection should be considered. In many countries, the role of the contact tracing is emphasized. However, there is sometimes a gap between the tracing and the analysis based on the compartmental model. The network theory bridges these two. In this context, the work by Choi et al. [30] to analyze the K-quarantine strategy based on the complex network theory is interesting. The network theory may suggest that networks with a weak hub structure are better to design a social system that is robust against pandemics.