Spread of variants of epidemic disease based on the microscopic numerical simulations on networks

Viruses constantly undergo mutations with genomic changes. The propagation of variants of viruses is an interesting problem. We perform numerical simulations of the microscopic epidemic model based on network theory for the spread of variants. Assume that a small number of individuals infected with the variant are added to widespread infection with the original virus. When a highly infectious variant that is more transmissible than the original lineage is added, the variant spreads quickly to the wide space. On the other hand, if the infectivity is about the same as that of the original virus, the infection will not spread. The rate of spread is not linear as a function of the infection strength but increases non-linearly. This cannot be explained by the compartmental model of epidemiology but can be understood in terms of the dynamic absorbing state known from the contact process.

Scientific Reports | (2022) 12:523 | https://doi.org/10.1038/s41598-021-04520-0 www.nature.com/scientificreports/ with a very large number of connections k, and the distribution of k is often scale-free, with a power distribution such as k −a13,14 . Barabási and Albert 14 proposed an algorithm to create networks that exhibit scale-free properties.
It is an algorithm that grows networks by selectively combining them. This mechanism is called "preferential attachment". The emergence of the Barabási-Albert network has led to rapid progress in the study of complex networks across many fields. The importance of the network structure in the analysis of epidemics was emphasized, and studies in several directions have been done [15][16][17][18][19][20][21] . In connection with COVID-19, stochastic simulations of the epidemic model have been performed 22,23 . The existence of absorbing states is one of the vital non-equilibrium processes on the complex network. Once the state has fallen into an absorbing state, the dynamics cannot escape from it [24][25][26] . The non-equilibrium absorbing phase transition is related to the directed percolation 27 . In epidemic spreading processes 6 , a fully healthy state can be regarded as an absorbing state in this sense. The contact process 28 and the susceptible-infected-susceptible (SIS) model 29 are often used to describe epidemic dynamics. It is well known that the SIS model can be mapped onto the logistic equation. In the epidemic scenario of the SIS model, individuals can be infected or susceptible. A phase transition between a disease-free (absorbing) phase and an active stationary phase is separated by an epidemic threshold. In the latter phase, a fraction of the population is infected. In the framework of the SIS model, the absence of an epidemic threshold was discussed for the spreading of infections on scale-free networks 15 .
The present authors 30 performed simulations for the microscopic SIR model on networks, and in particular, the relationship between the SIR model for macroscopic quantities and the corresponding microscopic SIR model was discussed. For the network, we discussed the difference between random networks and scale-free networks, and the role of hubs in scale-free networks was elucidated. The simulation method follows the method of Herrmann and Schwartz 22 . They discussed the role of the absorbing state in the SIR model.
In this paper, we use a microscopic model of infectious disease transmission on networks to simulate the spread of variants. The simulation method is an extension of the method used in the absence of variants 30 . A related problem is that of "competing epidemics" [31][32][33] . When there are two competing diseases and the relative infectivity of the two diseases is different, a phase diagram of the behavior of infection was investigated 32 .
In the following, we describe the simulation method of microscopic SIR model including variants in "Method of simulation of the microscopic SIR model" section of Methods. As for networks, we treat two types of networks; the Erdös-Rényi (ER) network 11,12 , a random network, and the Barabási-Albert (BA) network 14 , a scale-free network. We will discuss the relationship between the infectious strength of the variant and the spread of the infection. We emphasize the relation to the non-equilibrium absorbing phase transition.

Results
Simulation of the microscopic model for the variants on the ER network. We first consider the case of the ER network, a random network. The total number of nodes (individuals) is N = 10,000 and the average number of degrees is �k� = 8 . To perform the simulation of the microscopic model for the spread of infection, we choose the probability p of infection as p = 9/200 , which leads to β = �k�p = 0.36 in terms of a rate constant of the SIR model, which will be explained in "Differential equation of SIR model" section of Methods. The average infected period, 1/γ , is chosen as 5.0 days, which is realized by a Poisson distribution with an average value of 5.0. The basic reproduction number R 0 = β/γ , Equation (7), becomes R 0 = 1.8.
As a reference system, we deal with the case where there are no variants. The initial condition ( t = 0 ) is that 10 randomly selected individuals are infected. The time evolution of the number of individuals of the three types (S, I, and R) is shown in Fig. 1. We performed simulations for 100 samples, and the time evolution of all samples are plotted in Fig. 1a. We observe a variation in the time evolution for each sample. On the other hand, Fig. 1b is an average plot of 100 samples. 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. The microscopic SIR model on a random network is considered to reproduce the macroscopic SIR model almost quantitatively. The final value of the total number of infected individuals ( R(∞) ) is 0.732 from Eq.   34 , in the case of R 0 = β/γ = 1.8 , and the measured value for the present microscopic model on the ER network is about 0.77. As a related argument, the final size was discussed in the framework of a stochastic epidemic model 35 . In the Reed-Frost model 36,37 , which is a typical example of a stochastic epidemic model, the final size can be treated analytically. Britton et al. discussed the final size of stochastic epidemic models on networks 38 . In the previous study 30 , for a small number of initially infected individuals, some examples showed the behavior such that the infection vanishes quickly and does not spread throughout the network. This behavior is regarded as the absorbing state [24][25][26] in the contact process 28 . We chose 10 initially infected cases to avoid the situation of the absorbing state. We next consider the effects of variants. Suppose that 10 susceptible individuals are infected with the variant due to external factors at t = 21 . We choose 10 individuals randomly among the non-infected. The variant is assumed to be 3.0 times more infectious with β ′ = 1.08 , and the average infection period is chosen as the same value as that of the original virus, 1/γ ′ = 5.0 . Then, the basic reproduction number of the variant becomes R ′ 0 = 5.4 . The situation at t = 21 is that about 900 individuals are infected, 500 individuals are recovered, and 8600 individuals are not infected. In    To summarize what we have shown in this subsection, we can say the following. When the number of infected individuals reaches about 900 (the total number is 10,000), a small number of individuals infected with the variant are added. If the variant is highly infectious, it will spread to the wide space from the 10 additional infected individuals. By systematically examining the variants whose basic reproduction numbers are 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, and 3.0 times that of the original species, we showed that the infection increases non-linearly with the basic reproduction number. The lack of infection to spread in the case of variants whose infectivity is not much different from that of the original species is due to the dynamic effect of infection, that is, the absorbing state.
Simulation of the microscopic model for the variants on the BA network. Next, we consider the case of the BA network, a scale-free network. We choose the network which is similar to the case of the ER network. The total number of nodes (individuals) is N = 10,000 and the average number of degrees is �k� = 8 . The spread of infection is more rapid in the case of scale-free networks because of a hub structure 30 . Thus, for the probability of infection p, we choose a value smaller than that for the ER network, that is, p = 1/25 . This leads to β = �k�p = 0.32 in terms of a rate constant of the SIR model. The average infected period, 1/γ , is again chosen as 5.0 days. Then, the basic reproduction number becomes R 0 = β/γ = 1.6.
We treat the case where there are no variants as a reference system. The conditions are the same as the case of the ER network. The time evolution of the number of individuals of the three types (S, I, and R) is shown in Fig. 6. We performed simulations for 100 samples. We plot the time evolution of all 100 samples in Fig. 6a, whereas the average of 100 samples is plotted in Fig. 6b.
We turn to the investigation of the effects of variants. Suppose that 10 susceptible individuals are infected with the variant due to external factors at t = 10 . The variant is assumed to be 3.0 times more infectious with β ′ = 0.96 . The average infection period is chosen as 1/γ ′ = 5.0 , which leads to R ′ 0 = 4.8 . The situation at t = 10 is that about 620 individuals are infected, 100 individuals are recovered, and 9280 individuals are not infected. This is similar to the case of ER network at t = 21 shown in Fig. 2

Discussion
We studied the spread of variants of the virus by performing numerical simulations of the microscopic model on the network. In the middle of a simulation of infectious disease on the network, we added a variant that is more transmissible than the original lineage. When a highly infectious variant is added, the variant spreads quickly. It is noteworthy that the rate of spread is not linear in the infectivity of the variant but it is non-linear. If a variant is slightly more transmissible than the original lineage, the virus does not spread widely. It is related to a non-equilibrium phase transition of the epidemic dynamics on the complex network between a disease-free (absorbing) state and an active stationary phase [24][25][26] . This cannot be accounted for by the compartmental model of epidemiology. We performed simulations both on the ER model, a random network, and the BA model, a scalefree network. It was shown that the existence of the hub for a scale-free network stimulates the rapid increase of the infection 30 . But the effect of the variants for the BA network is essentially similar to the case of the ER network. It should also be noted that there is a lot of fluctuation in the spread of the more infectious variants. The greater fluctuation is due to the network structure. In summary, the reason why the more infectious variants spread, while the less infectious variants do not, has been clarified.
Although the present work is for a model system with some parameters, we may learn some lessons for a realworld strategy to fight against the pandemic. The spread of a slightly more infectious variant is not significant. Variants that are more infectious than the original lineage will spread the infection in a non-linear fashion. Of course, whether or not the variant is severe is an important factor. It can be said that the major role of overseas quarantine is to identify variants. For highly infectious variants, it is necessary to identify them by genetic testing and to isolate infected individuals. To conduct special contact tracing is of particular importance in the case that they spread throughout the city. Those variants must be monitored more carefully. Measures need to be taken to control the spread of the virus and its variants. Travel restrictions can be effective in reducing the spread of the virus. Contact with people over long distances is associated with the presence of hubs in the network. As shown in the previous paper 30 , isolation of hub sites can lead to transmission control. This is also the case for variants, for which travel restrictions are effective. However, to ultimately control the epidemic, it may be necessary to greatly accelerate vaccine roll-out.

Methods
Differential equation of SIR model. We perform simulation of the microscopic model of epidemic disease, which corresponds to the compartmental model of macroscopic variables. We consider SIR model 4-7,30 as a compartmental model. 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. Let S(t), I(t), and R(t) denote the number of individuals in three states S, I, and R as a function of time t. In the SIR model, the time evolution is described by the following simultaneous differential equations.
where β is the rate of infection and γ 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, β is given as kp. We assume that the total number of individuals is set to be constant such that Here, we do not consider the birth and death processes.
The SIR Eqs. (1)∼ (3), is said to have been given an exact solution by Harko et al. in 2014 39 . However, from the standpoint of differential equation theory, it was known to be solvable before that 40 . From Eqs. (1) and (2), we obtain The integration with respect to S yields where C is an integral constant. If we insert Eq. (6) into Eq. (1), we obtain a closed expression for S(t). We note that the ratio of β and γ appears for rate constants. The number R 0 defined by is known as the basic reproduction number 6,41,42 , and the number of infected individuals increases when R 0 > 1 , whereas it decreases when R 0 < 1.
It is noteworthy that I(t) is the number of infected individuals, not the number of new infected individuals announced every day. It should be noted that the number of new infections can be discussed in the framework of the SIR model 7 .
We here make a comment on the final size equation. We have the relation for the final value of R(∞) . In the limit of N 1 → N , we obtain This relation is known as the final size equation 34 . Method of simulation of the microscopic SIR model. Let us consider a simulation of microscopic infectious disease propagation in a network, in which each node assumes the state of infectious disease and changes its state stochastically. The rules for updating the state correspond to the macroscopic SIR model given by the differential equations, and the parameters are chosen to be equivalent so that the macroscopic and microscopic SIR models can be compared. The actual procedure, which is similar to that by Herrmann and Schwartz 22 , is as follows 30 : 1. Generate a network. 2. At t = 0 , an individual or individuals are initially infected (I). 3. A susceptible individual (S) will be infected (I) with a probability p if a connecting individual (one of k) is infected (I). We do this for all nodes (non-infected) and for all connections (k) of each node. In terms of the SIR model, the parameter β is β = kp.      www.nature.com/scientificreports/ 4. An infected individual (I) will recover (R) in 1/γ days on average. The infected period is chosen by a Poisson distribution with an average of 1/γ. 5. At each time t, the processes 3, 4 are repeated. 6. The time sequence obtained from the above procedure is regarded as a single sample. Simulations are performed for several samples.
In the simulation to study the effect of the variant, we will assume that some individuals are infected by the variant at time t = t var . We will denote by I′ those who are infected by the variant and by R′ those who recover from the variant. Modify the processes 2-4 of the simulation as follows: 2′ At t = t var , an individual or individuals are infected with the variant (I′). 3′ A susceptible individual (S) will be infected by the original virus (I) with a probability p if a connecting individual is infected by the original virus (I). If a connecting individual is infected by the variant (I′), a susceptible individual (S) will be infected by the variant (I′) with a probability p ′ . 4′ An infected individual by the original virus (I) will be recovered (R) in 1/γ days on average. The infected period is chosen by a Poisson distribution with an average of 1/γ . An infected individual by the variant of the virus (I′) will be recovered (R′) in 1/γ ′ days on average.