A computational model of epidemic process with three variants on a synthesized human interaction network

Virus mutations give rise to new variants that cause multiple waves of pandemics and escalate the infected number of individuals. In this paper, we develop both a simple random network that we define as a synthesized human interaction network and an epidemiological model based on the microscopic process of disease spreading to describe the epidemic process with three variants in a population with some features of social structure. The features of social structure we take into account in the model are the average number of degrees and the frequency of contacts. This paper shows many computational results from several scenarios both in varying network structures and epidemiological parameters that cannot be obtained numerically by using the compartmental model.


Scenario 1
The network structure of Scenario 1 with k = 7 and ω = 5 is shown in Fig. 1.The black histogram shows the degree distribution of the network.The blue one shows the contact distribution of the network.We can refer to the nodes with high degrees as hubs 24 .We will call the node a hub if its degree k ≥ 14 .The role of hubs in the epidemic process has been studied in 25 .
From the black histogram in Fig. 1, we can infer that there are fewer nodes with no degree.We call them isolated nodes.Almost all nodes in the network are connected.In addition, there are also many hubs in the network with degree 14-27.From the blue histogram in Fig. 1, we can infer that there are a lot of nodes with very high frequency of contact.Even, there are some nodes with contacts in the range 40-80.It makes sense if the disease spreads quickly as shown in Figure 2a.The figure is the time evolution of the number of individuals of all states (S, I, R) with 95% confidence interval.We divide this scenario into four cases.
In the first case we set µ 1→2 ∼ Bernoulli(0.005),µ 2→3 ∼ Bernoulli(0.005)and β 1 = β 2 = β 3 = 0.01 .The emergence of variants causes multiple waves as shown in Fig. 2a as the red line.We can see the emergence of each variant and their time evolution in detail in Fig. 2b.The figure is the average plot of 100 samples.The infection peak of Variant 1 is 4480 nodes at the time t = 46 .The infection of Variant 2 is 614 nodes at the time t = 94 .The infection of Variant 3 is 197 at the time t = 227 .Variant 1 ends at time t = 118 .Variant 2 starts to spread at time t = 37 and ends at time t = 227 .Variant 3 starts to spread at time t = 140 and ends at time t = 337 .We can see the final value of the total number of infected individuals for each variant in Fig. 2c.From the total number of infected individuals of Variant 1 AI T 1 400 = 8250 and the length of epidemic time of Variant 1 is t e 1 = 117 days, we can find the infection spread number of Variant 1 ρ 1 = 0.241 .For Variant 2 and Variant 3, we obtain AI T 2 400 = 2305 with t e 2 = 190 and AI T 3 400 = 992 with t e 3 = 197 .Then, we have ρ 2 = 0.109 and ρ 3 = 0.049 .We can compare the number of infected neighbors to the degree of each node in Fig. 2d.The higher the degree of a node the more neighbor will be infected.
In the second case we still set µ 1→2 ∼ Bernoulli(0.005),µ 2→3 ∼ Bernoulli(0.005)but β 1 < β 2 < β 3 .We choose β 1 = 0.01, β 2 = 0.02 and β 3 = 0.03 to show how different infection rates play roles in the epidemic process.The dynamic of the epidemic process in this case is shown in Fig. 3a.We can see that Variant To see how different infection rates from each variant play roles in a condition when the probability of the emergence of new variants is smaller, we set the fourth case where µ 1→2 ∼ Bernoulli(0.001),µ 2→3 ∼ Bernoulli(0.001)and β 1 < β 2 < β 3 .We still choose β 1 = 0.01, β 2 = 0.02 and β 3 = 0.03 .The dynamic of the epidemic process in this case is shown in Fig. 5a.The infected individuals are still dominated by Variant 1 as shown in Fig. 5b.The infection peak of Variant 1 is 4448 nodes at the time t = 45 .The infection peak of Variant 2 is 376 nodes at the time t = 70 .The infection peak of Variant 3 is 181 at the time t = 192 .Variant 2 starts to spread at time t = 30 and ends at time t = 147 .Variant 1 ends at time t = 117 .Variant 3 starts to spread at time t = 174 and ends at time t = 217 .The smaller probability of the emergence of new variants causes the final values of the total number of infected individuals for Variant 2 and Variant 3 to be smaller than the total number of infected individuals for Variant 1 as shown in Fig. 5c.The total number of infected individuals of Variant 1 AI T 1 400 = 8249 and the epidemic duration t e 1 = 116 days give the infection spread number of Variant 1 ρ 1 = 0.239 .For Variant 2 and Variant 3, we obtain AI T 2 400 = 939 with t e 2 = 117 and AI T 3 400 = 195 with t e 3 = 43 .We have ρ 2 = 0.027 and ρ 3 = 0.002 .Although a smaller probability of the emergence of new variants can reduce the total number of infected individuals of new variants, the higher infection rates still give a significant increase in the total number of infected individuals for each variant.

Scenario 2
The network structure of Scenario 2 where k = 3 and ω = 5 is shown in Fig. 6.In this scenario, we contain the spread of infection in the network by implementing self-quarantine without physical distancing.We reduce individual connectivity without reducing their contact frequency.
The black histogram in Fig. 1 shows that there is no hub in this network.From the blue histogram in Fig. 1, there are fewer nodes with contacts in the range 30 − 40 which is the high contacts in the network.The dynamic of the epidemic process is shown in Fig. 7a.The figure is the time evolution of the number of individuals of all states (S, I, R) with 95% confidence interval.Here, we still divide this scenario into four cases.
In the first case of Scenario 2, we set µ 1→2 ∼ Bernoulli(0.005),µ 2→3 ∼ Bernoulli(0.005)and β 1 = β 2 = β 3 = 0.01 .The number of infected individuals from all variants indicated by the red line is shown 2 400 = 688 with t e 2 = 336 and AI T 3 400 = 319 with t e 3 = 317 .Then, we have ρ 2 = 0.058 and ρ 3 = 0.025 .In this case, reducing the connectivity of each individual by reducing the average degree prolongs the length of the epidemic time of each variant.We can compare the number of infected neighbors to the degree of each node in Fig. 7d.The higher the degree of a node the more neighbor will be infected.
In the second case of Scenario 2, we still set µ 1→2 ∼ Bernoulli(0.005),µ 2→3 ∼ Bernoulli(0.005)but β 1 < β 2 < β 3 .We choose β 1 = 0.01, β 2 = 0.02 and β 3 = 0.03 to show how different infection rates play roles in the epidemic process when the average degree is reduced.The dynamic of the epidemic process in this case is shown in Fig. 8a.There is an increase in the number of infections due to Variant 2 and Variant 3 as shown in  3 400 = 1 .We have ρ 2 = 0.019 and ρ 3 = 0.
To see how different infection rates from each variant play roles in a condition when the probability of the emergence of new variants is smaller, we set the fourth case of Scenario 2 where µ 1→2 ∼ Bernoulli(0.001),µ 2→3 ∼ Bernoulli(0.001)and β 1 < β 2 < β 3 .We still choose β 1 = 0.01, β 2 = 0.02 and β 3 = 0.03 .The dynamic of the epidemic process in this case is shown in Fig. 10a

Scenario 3
We assume that the network structure of Scenario 3 where k = 7 and ω = 1 is the epidemic process in the net- work with measured physical distancing but no self-quarantine as shown in Fig. 11.In this scenario, we contain the spread of infection in the network by implementing physical distancing without self-quarantine.We reduce the individual frequency of contact without reducing the connectivity.Here, we still divide this scenario into four cases.The black histogram in Fig. 11 is the same as the blue histogram.There are quite a lot of individuals who are hubs and have frequency contact of 14 − 28 .The dynamic of the epidemic process is shown in Fig. 12a.The figure is the time evolution of the number of individuals of all states (S,I, R) with 95% confidence interval.
In the first case of Scenario 3, we set µ 1→2 ∼ Bernoulli(0.005),µ 2→3 ∼ Bernoulli(0.005)and β 1 = β 2 = β 3 = 0.01 .The number of infected individuals from all variants indicated by the red line is shown in Fig. 12a.Although the infected individuals are still dominated by the infection of Variant 1, the total number    12c.From the total number of infected individuals of Variant 1 AI T 1 400 = 1925 and the length of epidemic time of Variant 1 is t e 1 = 363 days, we can find the infection spread number of Variant 1 ρ 1 = 0.175 .For Variant 2 and Variant 3, we obtain AI T 2 400 = 592 with t e 2 = 362 and AI T 3 400 = 90 with t e 3 = 195 .Then, we have ρ 2 = 0.054 and ρ 3 = 0.004 .In this case, reducing the frequency of contact also prolongs the length of the epidemic time of each variant.We can compare the number of infected neighbors to the degree of each node in Fig. 12d.
In the second case of Scenario 3, we still set µ 1→2 ∼ Bernoulli(0.005),µ 2→3 ∼ Bernoulli(0.005)but β 1 < β 2 < β 3 .We still choose β 1 = 0.01, β 2 = 0.02 and β 3 = 0.03 to show how different infection rates play roles in the epidemic process when the frequency of contact is reduced.The dynamic of the epidemic process in this case is shown in Fig. 13a.There is an increase in the number of infections due to Variant 2 and Variant 3 as shown in Fig. 13b  1 400 = 1617 and the length of epidemic time t e 1 = 399 days, we can find the infection spread number of Variant 1 ρ 1 = 0.161 .For Variant 2 and Variant 3, we obtain AI T 2 400 = 3472 with t e 2 = 381 and AI T 3 400 = 1794 with t e 3 = 341 .We obtain ρ 2 = 0.331 and ρ 3 = 0.153 .In this case, if the infection rate of the new variant is larger than the previous one, we can see that not only does the total number of infected individuals increase but also the length of the epidemic time is longer.Now we set the third case of Scenario 3 where µ 1→2 ∼ Bernoulli(0.001),µ 2→3 ∼ Bernoulli(0.001)and β 1 = β 2 = β 3 = 0.01 .In this case, we want to see what happens if we reduce the probability of the emergence of new variants.The dynamic of the epidemic process in this case is shown in Fig. 14a.We can see that the time 1 400 = 1599 and the epidemic duration t e 1 = 303 days give the infection spread number of Variant 1 ρ 1 = 0.0121 .Since there is no spread of Variant 2, we write ρ 2 = 0 .For Variant 3, we obtain AI T 2 400 = 51 with t e 2 = 141 and then ρ 3 = 0.002.To see how different infection rates from each variant play roles in a condition when the probability of the emergence of new variants is smaller and the frequency of contact is reduced, we set the fourth case of Scenario 3 where µ 1→2 ∼ Bernoulli(0.001),µ 2→3 ∼ Bernoulli(0.001)and β 1 < β 2 < β 3 .We still choose β 1 = 0.01, β 2 = 0.02 and β 3 = 0.03 .The dynamic of the epidemic process in this case is shown in Fig. 15a.The time evolution of infected individuals for each variant is shown in Fig. 15b.The infection peak of Variant 1 is 258 nodes at the time t = 118 .The infection peak of Variant 2 is 187 nodes at the time t = 256 .The infection peak of Variant 3 is 100 at the time t = 226 .Variant 1 ends at time t = 323 .Variant 2 starts to spread at time t = 46 and ends at time t = 375 .Variant 3 starts to spread at time t = 73 and ends at time t = 173 .The total number of infected individuals for each variant is shown in Fig. 15c.The total number of infected individuals of Variant 1 AI T 1 400 = 1751 and the epidemic duration t e 1 = 322 days give the infection spread number of Variant 1 ρ 1 = 0.141 .For Variant 2 and Variant 3, we obtain AI T 2 400 = 1357 with t e 2 = 329 and AI T 3 400 = 328 with t e 3 = 100 .We have ρ 2 = 0.111 and ρ 3 = 0.008.

Scenario 4
We assume that the network structure of Scenario 4 where k = 3 and ω = 1 is the epidemic process in the net- work with measured physical distancing and self-quarantine as shown in Fig. 16.In this scenario, we contain the spread of infection in the network by implementing both self-quarantine and physical distancing.We reduce both the connectivity and the individual frequency of contact.The black histogram in Fig. 16 is the same as the blue histogram.There is no hub in this network.The dynamic of the epidemic process is shown in Fig. 17a.The figure is the time evolution of the number of individuals of all states (S, I, R) with 95% confidence interval.Here, we still divide this scenario into four cases.
For the first case of Scenario 4, we set µ 1→2 ∼ Bernoulli(0.005),µ 2→3 ∼ Bernoulli(0.005)and β 1 = β 2 = β 3 = 0.01 .The number of infected individuals from all variants indicated by the red line is shown in Fig. 17a.The emergence of new variants is shown in Fig. 17b.The figure is the average plot of 100 samples.The infection peak of Variant 1 is 8 nodes at the time t = 16 .There is no spread of Variant 2 and Variant 3. The spread of Variant 1 ends at time t = 39 .We can see the final value of the total number of infected individuals for each variant in Fig. 17c.From the total number of infected individuals of Variant 1 AI T 1 400 = 11 and the length of epidemic time of Variant 1 is t e 1 = 38 days, we can find the infection spread number of Variant 1 ρ 1 = 0.0001 .We can compare the number of infected neighbors to the degree of each node in Fig. 17d.
In the second case of Scenario 4, we still set µ 1→2 ∼ Bernoulli(0.005),2→3 ∼ Bernoulli(0.005)but β 1 < β 2 < β 3 .We still choose β 1 = 0.01, β 2 = 0.02 and β 3 = 0.03 to show how different infection rates play roles in the epidemic process when the two containment measures are implemented.The dynamic of the epidemic process in this case is shown in Fig. 18a.There is an increase in the number of infections due to Variant 3 as shown in Fig. 18b.The infection peak of Variant 1 is 9 nodes at time t = 16 .There is no spread of Variant 2. The infection peak of Variant 3 is 34 at time t = 147 .Variant 1 ends at time t = 47 .Variant 3 starts to spread at time t = 76 and ends at time t = 225 .The final value of the total number of infected individuals for each variant is shown in Fig. 18c.From the total number of infected individuals of Variant 1 AI T 1 400 = 12 and the length of epidemic time t e 1 = 46 days, we can find the infection spread number of Variant 1 ρ 1 = 0.0001 .Since there is no spread of Variant 2, we write ρ 2 = 0 .For Variant 3, we obtain AI T 3 400 = 134 with t e 3 = 149 .We obtain ρ 3 = 0.005 .In this case, reducing the connectivity and the frequency of contact is not strong enough to contain the spread of variants if the infection rate of the new variant is more infectious.Now we set the third case of Scenario 3 where µ 1→2 ∼ Bernoulli(0.001),µ 2→3 ∼ Bernoulli(0.001)and β 1 = β 2 = β 3 = 0.01 .In this case, we want to see what happens if we reduce the probability of the emergence of new variants in the network with reduced connectivity and frequency of contact.The dynamic of the epidemic process in this case is shown in Fig. 19a.There is no infection by Variant 2 and Variant 3 as shown in Fig. 14b   To see how different infection rates from each variant play roles in a condition when the probability of the emergence of new variants is smaller and all contained measures are implemented, we set the fourth case of Scenario 3 where µ 1→2 ∼ Bernoulli(0.001),µ 2→3 ∼ Bernoulli(0.001)and β 1 < β 2 < β 3 .We still choose β 1 = 0.01, β 2 = 0.02 and β 3 = 0.03 .The dynamic of the epidemic process in this case is shown in Fig. 20a

Discussion
We developed a computational model of disease spreading with three variants that can capture the microscopic process in a synthesized human interaction network.In generating the network we need some social features in human interaction such as connectivity and contact frequency or interaction intensity.Both features determine the network structures such that every individual has an unequal chance of coming into contact with any other individual in the population.Of course, it is more realistic to describe a human population although not a real representation of the social network of the human population.Therefore, discussing the spread of contact-based diseases in a human population is more relevant in smaller-scale environments.However, we need to develop some critical aspects.Because it is a microscopic process and we brought the model to an individual scale, as a consequence, unlikely the classical SIR model, our model does not use merely infection rate β .we also do not use a recovery rate as usual we use it in the common SIR model.Incorporating more epidemiological parameters such as the average duration of the infectious period, the probability of the emergence of new variants, the length of the infection period, the length of the incubation period, and many more, has also consequences.In our model, we perform some tricks to the parameters such that they can be incorporated into the model and more applicable for both mathematics and computation purposes.We can also generalize our model to any number of variants and incorporate more epidemiological parameters such that the model is more complicated and realistic.Hopefully, there will be a study to analyze our model rigorously.
Our model has rich computational results by varying several parameters such as the number of individuals in the network, the average degree, the range of contact frequency, the values of epidemiological parameters, the initial conditions, the probability of the emergence of the new variants, and the infection rate of each variant.In the simulation, we chose to vary the average degree and contact frequency to elucidate the role of social network structures in disease spreading by keeping several epidemiological parameters and initial conditions constant.By varying the average degree and contact frequency we also address the containment measures implemented into the network to contain the spreads.The simulation results gave us a particular insight into how the social network structures determine epidemic behaviors such as the peak values of infection curves, the time to reach them, the increasing and decreasing rate of the infection curves, the infection spread numbers, and the final size of the epidemic.Here, we try to describe the disease spreading between the population that implements the containment measures and the one that does not.Next, from the model, we can develop a model that incorporates more interventions to derive efficient and realistic strategies for curbing the spread of disease throughout the network.

Method
The first step to a synthesized human interaction network with size N and an average degree k is using the simplest examples of a random network in which we fix only the number of nodes N and the number of edges L = kN/2 .We take N nodes and place L edges among them at random 19 .The next step is to attach each edge with a natural number which is randomly chosen from set W = {1, 2, . . ., ω} ⊆ N , where ω is the maximum of contact frequency.Here, we obtain the weighted adjacency matrix A = ω ij with ω ij ∈ W ∪ {0} to represent the network.We assume that the network is undirected.If ω ij = 0 , then there is no interaction or link between node i and node j.From A , we can obtain N i as a set containing all neighbors of node i.
To simplify the discrete-time process we assume that one day is a unit time step and we write t = 0, 1, 2, . . ., t f as a time discretization.We denote t f as the duration of the epidemic process.At time t, a node i can be in a state X i t belonging to a finite set of states = S i t , I i The state X i t ∈ {0, 1} and Since we put three variants into the network, we write v = 1, 2, 3 to denote the variants, and the original virus is Variant 1.All states are explained in the following table (Table 1).
To introduce the infections due to the Variant 1 at the initial time, we choose randomly N 0 individuals.Variant 2 will emerge from Variant 1 by moving the infected individual in the state I i 1 t to the state I i 2 t which is denoted by a Bernoulli random variable µ 1→2 with a small probability.Variant 3 will emerge from Variant 2 with a similar (1) for all nodes i and for all time t procedure.The infected individual in the state I i 2 t will move to the state I i 3 t which is denoted by a Bernoulli random variable µ 2→3 with a small probability.
For each variant v = 1, 2, 3 , we let ξ v t as a Bernoulli random variable with a probability β v to denote viral transmission between two nodes.We define that β v is the average rate of infection in one contact for Variant v.The process of viral transmission from one of the infected neighbors of node i to node i can be illustrated in Fig. 21.The viral transmission to node i at time t for each variant can be expressed by the following equation.
Note that the Eq. ( 2) describes how viral transmission occurs at the individual level or microscopic scale.Not only does the equation consider all neighbors for each individual but also considers the frequency of contact with all neighbors.It also tells us that the bigger the number of neighbors and weight, the more likely a node i to get an infection.
To make clear when the neighbors of node i contain the infected individuals with more than one variant, the viral transmission to node i occurs when the first infected individual succeeds in transmitting the infection.With this rule, we avoid an individual getting more infections at one time.Thus, we write We assume that each individual will gain permanent immunity after she or he recovers from one variant but lose immunity if she or he gets infected by other variants.
Unlike the common SIR model, we do not use the recovery rate.Instead, we use the parameters t i v , τ v and the unit-step functions u t, t i v and u t, t i v + τ v for each variant v.The parameter t i v is the time when a node i gets infected by Variant v for the first time.Meanwhile, the parameter τ v is the average duration of the infectious period for each variant.The functions u t, t i v and u t, t i v + τ v will move the individual from the susceptible state to the infected state and from the infected state to the recovery state respectively.The process of state transition for each node in the network is shown in Fig. 22 and the process is governed by the following system of equations. (2) for all node i and for all time t   In epidemiology, we used to consider basic reproduction number R 0 to describe the ability of the disease to spread in a population.It is an average number of additional infected individuals produced by such individual passes the disease onto before it recovers 19 .Referring to the classic SIR model, when R 0 > 1 , infected individu- als will grow geometrically, causing an epidemic 26 .When R 0 < 1 , the disease will die out.It can be calculated straightforwardly by finding the ratio of infection and recovery rate.However, it is not easy to find the number from our model.We do not have a recovery rate for each variant.In addition, our model uses three variants and relies on the microscopic process in discrete time.We consider the epidemic process on an individual scale.The size of the population is limited to a certain number which is much smaller than the size of the population in the common SIR model.Trying to find the basic reproduction number of the model in this research is no longer relevant.Thus, we do not deal with the basic reproduction number for our model because the model we developed is fully different from the common SIR model.Instead of using it, we define a number that represents how big and long an epidemic spreads.We call the number the infection spread number and denote it with ρ .We define it explicitly for each variant as the following.
where t e v is the length of epidemic time which is defined as when the infection starts to spread and when it ends.The infection spread number has no dimension.We can interpret the number to infer how severe the epidemic is.The larger its value takes, the bigger and the longer, the size and duration of the epidemic will be.( 16)

Figure 1 .
Figure 1.The network structure of Scenario 1.The black histogram shows the degree distribution of the network with k = 7 and ω = 5 .It shows that the network has a lot of hubs and fewer isolated nodes.The blue histogram shows the contact distribution.More than half of the population have a contact frequency of 20 − 90.

Figure 2 .
Figure 2. The dynamics of the epidemic process with three variants on a synthesized human interaction network for the first case of Scenario 1.(a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of all variants are indicated by I 1 46 = 4480 , I 2 94 = 614 , and I 3 227 = 197 respectively.(c) The total number of infections for all variants are AI T 1 400 = 8250 , AI T 2 400 = 2305 , and AI T 3 400 = 992 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread number are ρ 1 = 0.241 , ρ 2 = 0.109 , and ρ 3 = 0.049 respectively.
Fig. 8b compared to Fig. 7b.The emergence of new variants is earlier than the emergence of new variants in the first case.The infection peak of Variant 1 is 337 nodes at the time t = 105 .The infection peak of Variant 2 is 451 nodes at the time t = 101 .The infection peak of Variant 3 is 251 at the time t = 145 .Variant 1 ends at time t = 332 .Variant 2 starts to spread at time t = 29 and ends at time t = 306 .Variant 3 starts to spread at time t = 42 and does not end until the last time.The final value of the total number of infected individuals for each variant is shown in Fig. 8c.From the total number of infected individuals of Variant 1 AI T 1 400 = 2195 and the length of epidemic time t e 1 = 331 days, we can find the infection spread number of Variant 1 ρ 1 = 0.182 .For Variant 2 and Variant 3, we obtain AI T 2 400 = 2610 with t e 2 = 277 and AI T 3 400 = 2078 with t e 3 = 358 .We obtain ρ 2 = 0.181 and ρ 3 = 0.186.Now we set the third case of Scenario 2 where µ 1→2 ∼ Bernoulli(0.001),µ 2→3 ∼ Bernoulli(0.001)and β 1 = β 2 = β 3 = 0.01 .The dynamic of the epidemic process in this case is shown in Fig. 9a.The infection of Variant 2 and Variant 3 is smaller than the two previous cases as shown in Fig. 9b.The infection peak of Variant 1 is 454 nodes at the time t = 106 .The infection peak of Variant 2 is 29 nodes at the time t = 334 .There is no spread of Variant 3. Variant 1 ends at time t = 270 .Variant 2 starts to spread at time t = 115 and does not end until the last time.The final value of the total number of infected individuals for Variant 2 and Variant 3 are

Figure 3 .
Figure 3.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the second case of Scenario 1.(a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of all variants are indicated by I 1 46 = 4241 , I 2 87 = 977 , and I 3 152 = 519 respectively.(c) The total number of infections for all variants are AI T 1 400 = 8181 , AI T 2 400 = 4675 , and AI T 3 400 = 2510 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread number are ρ 1 = 0.272 , ρ 2 = 0.218 , and ρ 3 = 0.122 respectively.
. The time evolution of infected individuals for each variant is shown in Fig. 10b.The infection peak of Variant 1 is 429 nodes at the time t = 102 .The infection peak of Variant 2 is 241 nodes at the time t = 140 .The infection peak of Variant 3 is 56 at the time t = 169 .Variant 1 ends at time t = 303 .Variant 2 starts to spread at time t = 59 and ends at time t = 303 .Vari- ant 3 starts to spread at time t = 125 and ends at time t = 226 .The total number of infected individuals for each variant is shown in Fig. 10c.The total number of infected individuals of Variant 1 AI T 1 400 = 2547 and the epidemic duration t e 1 = 302 days give the infection spread number of Variant 1 ρ 1 = 0.192 .For Variant 2 and Variant 3, we obtain AI T 2 400 = 1000 with t e 2 = 244 and AI T 3 400 = 143 with t e 3 = 101 .We have ρ 2 = 0.061 and ρ 3 = 0.003.

Figure 4 .Figure 5 .
Figure 4.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the third case of Scenario 1.(a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of all variants are indicated by I 1 46 = 4482 , I 2 104 = 158 , and I 3 265 = 99 respectively.(c) The total number of infections for all variants are AI T 1 400 = 8256 , AI T 2 400 = 496 , and AI T 3 400 = 166 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread number are ρ 1 = 0.243 , ρ 2 = 0.022 , and ρ 3 = 0.003 respectively.

Figure 6 .
Figure 6.The network structure of Scenario 2. The black histogram shows the network's degree distribution with k = 3 and ω = 5 .The highest degree is 13.The network is dominated by the nodes with degree 2. The blue histogram shows the distribution of the total number of contacts.More than half of the population have a contact frequency of 10 − 40.
compared to Fig. 12b.The emergence of new variants is earlier than the emergence of new variants in the first case.The infection peak of Variant 1 is 196 nodes at the time t = 139 .The infection peak of Variant 2 is 540 nodes at the time t = 100 .The infection peak of Variant 3 is 297 at the time t = 211 .Variant 1 spreads until the last time.Variant 2 starts to spread at time t = 19 and has no end.Variant 3 starts to spread at time t = 59 and has no end either.The final value of the total number of infected individuals for each variant is shown in Fig. 13c.From the total number of infected individuals of Variant 1 AI T

Figure 7 .
Figure 7.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the first case of Scenario 2. (a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of all variants are indicated by I 1 115 = 381 , I 2 212 = 74 , and I 3 384 = 30 respectively.(c) The total number of infections for all variants are AI T 1 400 = 2426 , AI T 2 400 = 688 , and AI T 3 400 = 319 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread numbers are ρ 1 = 0.198 , ρ 2 = 0.058 , and ρ 3 = 0.025 respectively.

Figure 8 .
Figure 8.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the second case of Scenario 2. (a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of all variants are indicated by I 1 105 = 337 , I 2 101 = 451 , and I 3 145 = 251 respectively.(c) The total number of infections for all variants are AI T 1 400 = 2195 , AI T 2 400 = 2610 , and AI T 3 400 = 2078 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread numbers are ρ 1 = 0.182 , ρ 2 = 0.181 , and ρ 3 = 0.186 respectively.
. The infection peak of Variant 1 is 9 nodes at the time t = 16 .There is no spread of Variant 2 and Variant 3. Thus ρ 2 = ρ 3 = 0 .Variant 1 ends at time t = 59 .The total number of infected individuals of Variant 1 AI T 1 400 = 14 and the epidemic duration t e 1 = 58 .We have ρ 1 = 0.0002.

Figure 9 .Figure 10 .
Figure 9.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the third case of Scenario 2. (a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of all variants are indicated by I 1 106 = 454 , I 2 334 = 29 , and I 3 139 = 1 respectively.(c) The total number of infections for all variants are AI T 1 400 = 2498 , AI T 2 400 = 277 , and AI T 3 400 = 1 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread numbers are ρ 1 = 0.168 , ρ 2 = 0.019 , and ρ 3 = 0 respectively.

Figure 11 .
Figure 11.The network structure of Scenario 3. The black histogram shows the network's degree distribution with k = 7 and ω = 1 .It is the same distribution as the blue histogram because the contact frequency is 1 for each link in the network.There are many hubs with a contact frequency of 14 − 28.
. The time evolution of infected individuals for each variant is shown in Fig. 20b.The infection peak of Variant 1 is 9 nodes at time t = 16 .The infection peak of Variant 2 is 5 nodes at time t = 158 .There is no infection of Variant 3. Variant 1 ends at time t = 41 .Variant 2 starts to spread at time t = 114 and ends at time t = 224 .The total number of infected individuals for each variant is shown in Fig. 20c.The total number of infected individuals of Variant 1 AI T 1 400 = 12 and the the length of epidemic time of Variant 1 t e 1 = 40 days give the infection spread number of Variant 1 ρ 1 = 0.0001 .For Variant 2 we obtain AI T 2 400 = 26 with t e 2 = 110 .We have ρ 2 = 0.0007.

Figure 12 .
Figure 12.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the first case of Scenario 3. (a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of all variants are indicated by I 1 122 = 271 , I 2 177 = 44 , and I 3 244 = 18 respectively.(c) The total number of infections for all variants are AI T 1 400 = 1925 , AI T 2 400 = 592 , and AI T 3 400 = 90 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread numbers are ρ 1 = 0.175 , ρ 2 = 0.054 , and ρ 3 = 0.004 respectively.

Figure 13 .
Figure 13.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the second case of Scenario 3. (a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of all variants are indicated by I 1 139 = 196 , I 2 100 = 540 , and I 3 211 = 297 respectively.(c) The total number of infections for all variants are AI T 1 400 = 1617 , AI T 2 400 = 3472 , and AI T 3 400 = 1794 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread number of all variants are ρ 1 = 0.161 , ρ 2 = 0.331 , and ρ 3 = 0.153 respectively.

Figure 14 .
Figure 14.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the third case of Scenario 3. (a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of all variants are indicated by I 1 126 = 226 , I 2 127 = 1 , and I 3 314 = 10 respectively.(c) The total number of infections for all variants are AI T 1 400 = 1599 , AI T 2 400 = 1 , and AI T 3 400 = 51 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread numbers are ρ 1 = 0.121 , ρ 2 = 0 , and ρ 3 = 0.002 respectively.

Figure 15 .
Figure 15.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the fourth case of Scenario 3. (a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of all variants are indicated by I 1 118 = 258 , I 2 256 = 187 , and I 3 226 = 100 respectively.(c) The total number of infections for all variants are AI T 1 400 = 1751 , AI T 2 400 = 1357 , and AI T 3 400 = 328 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread numbers are ρ 1 = 0.141 , ρ 2 = 0.111 , and ρ 3 = 0.008 respectively.

Figure 16 .
Figure 16.The network structure of Scenario 4. The black histogram shows the network's degree distribution with k = 3 and ω = 1 .The black histogram is the same as the blue one.There is no hub in this network.

Figure 19 .
Figure 19.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the third case of Scenario 4. (a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of Variant 1s is I 1 16 = 9 .(c) The total number of infections for all variants are AI T 1 400 = 14 , AI T2 400 = 0 , and AI T 3 400 = 0 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread numbers are ρ 1 = 0.0002 , ρ 2 = 0 , and ρ 3 = 0 respectively.

Figure 20 .
Figure 20.The dynamics of the epidemic process with three variants on a synthesized human interaction network for the fourth case of Scenario 4. (a) The time evolution of all states with 95% confidence interval.The red lines indicate all infected nodes from all variants.(b) The average of 100 infection plots for each variant.The infection peaks of Variant 1 and Variant 2 are I 1 16 = 9 and I 2 158 = 614 .(c) The total number of infections for all variants are AI T 1 400 = 12 , AI T 2 400 = 26 , and AI T 3 400 = 0 . (d) Comparison between the number of infected neighbors and each node's degree for all variants.The infection spread numbers are ρ 1 = 0.0001 , ρ 2 = 0.0007 , and ρ 3 = 0 respectively.

Figure 21 .
Figure 21.The viral transmission to node i from its neighbours.The red node represents an infected one.The blue nodes can be susceptible or recovered from previous variants.

Figure 22 .
Figure 22.The process of all possibilities of state transition for each node in the network at a single time t.

Table 1 .
The Description of states.State at time t when a node i is infected by Variant j ( I i vt = 1 ) or not ( I i vt = 0) R i vt State at time t when node i is recovered from Variant j ( R i vt = 1 ) or not ( R i vt = 0) t State at time t when a node i is susceptible ( S i t = 1 ) or not ( S i t