A data-driven network model for the emerging COVID-19 epidemics in Wuhan, Toronto and Italy.

The ongoing Coronavirus Disease 2019 (COVID-19) pandemic threatens the health of humans and causes great economic losses. Predictive modeling and forecasting the epidemic trends are essential for developing countermeasures to mitigate this pandemic. We develop a network model, where each node represents an individual and the edges represent contacts between individuals where the infection can spread. The individuals are classified based on the number of contacts they have each day (their node degrees) and their infection status. The transmission network model was respectively fitted to the reported data for the COVID-19 epidemic in Wuhan (China), Toronto (Canada), and the Italian Republic using a Markov Chain Monte Carlo (MCMC) optimization algorithm. Our model fits all three regions well with narrow confidence intervals and could be adapted to simulate other megacities or regions. The model projections on the role of containment strategies can help inform public health authorities to plan control measures.

Introduction

The development of international trade and tourism has accelerated the spatial spread of infectious diseases. The limited data available on emerging epidemics adds to the challenge of mitigating the spread of emerging infections [1]. The unprecedented Coronavirus Disease 2019 (COVID-19) outbreak began at the end of 2019. The number of reported cases keeps rising worldwide and thousands of lives have been claimed. This pandemic is having an enormous impact on world health, disturbing the stability of the societies, and triggers great economic losses. Predicting the future of the pandemic, assessing the impacts of current interventions, and evaluating the effectiveness of alternate mitigation strategies are of utmost importance for saving lives.

Mathematical models can be used to understand the dynamics of epidemics and help inform control strategies. A numerous number of models are being used to project the current COVID-19 pandemic. Ziff and Ziff analyzed the number of reported cases for Wuhan (China) and showed that the growth of the daily number of confirmed new cases indicates an underlying fractal or small-world network of connections between susceptible and infected individuals [2]. Wang et al. developed an SEIR model to estimate the epidemic trends in Wuhan, assuming the prevention and control measures were either sufficient or insufficient to control the epidemic [3]. Kucharski et al. combined a stochastic transmission model with data on cases of COVID-19 in Wuhan and international cases to estimate how the transmission had varied over time between January and February in 2020 [4].

Kraemer et al. analyzed the impact of interventions on the spread of COVID-19 in China using transportation data [5]. Chinazzi et al. used a global meta-population disease transmission model to project the impact of travel limitations on the national and international spread of the epidemic. They showed that the travel restriction of Wuhan, China had a more marked effect on the international scale than that on Mainland China [6].

Ferguson et al. found that optimal mitigation policies (combining home isolation of suspected cases, home quarantine of those living in the same household as suspected cases, and social distancing of the elderly and others at most risk of severe disease) might reduce peak healthcare demand by 2/3 and deaths by half [7]. Likewise, Hellewell et al. developed a stochastic transmission model and found that highly effective contact tracing and case isolation are enough to control a new outbreak of COVID-19 within three months in most scenarios [8].

Zhang et al. fitted the reported serial interval (mean and standard deviation) with a gamma distribution to estimate the basic reproduction number at the early stage of a COVID-19 outbreak, indicating the potential of second outbreaks [9]. Maier et al. developed a compartmental model dividing individuals into susceptible, exposed, removed, and quarantined symptomatically infected and showed that the distinctive subexponential increase of confirmed cases in mainland China could be explained as a direct consequence of containment policies that effectively deplete the susceptible population [10].

Most of these models are based on assuming the population is homogeneously mixing, that is, the contacts between people are random and uniformly distributed throughout the population. However, different individuals may have varying numbers of acquaintances and contacts in the real world. The important role that heterogenous contact networks play in the transmission dynamics of infectious diseases is often overlooked [11]. Models that take into account contact heterogeneity better represent the actual transmission network through a population and are more likely to capture the true epidemic dynamics.

Disease propagation is closely linked with the structure of social contact networks [12]. The ubiquitous diversity in contact patterns and heterogeneity among individuals depends on differences in social structures, spatial distances, and behavior [13]. The heterogeneity exists at a wide range of scales and leads to highly variable transmission dynamics of infectious diseases [14], [15].

Many real-world social networks can be characterized by a random Watts–Strogatz (WS) small-world network [16], [17]. In a small-world network, most nodes can be reached from every other node by a small number of hops or steps, even if they are not immediate neighbors. This type of network model allows us to adapt changes to some realistic network structures and examine the effects of control and intervention countermeasures such as social distancing, self-isolation, and personal protection. The framework and analysis can be applied to study the transmission dynamics in different regions and many other infectious diseases.

The COVID-19 epidemic in Wuhan ended in April, while the epidemics in the Greater Toronto Area (GTA, Canada) and the Italian Republic are continuing to grow. We fit the parameters of our network model to the confirmed cases in each of these regions. Although Wuhan, Toronto, and Italy differ in some ways, the way that SARS-CoV-2 is transmitted from one person to another is quite similar. Individuals may acquire infection from other infectious individuals, even if they do not contact each other directly. The Watts–Strogatz model supplies an ideal tool to study the spread of epidemics among individuals even if their locations are not considered. We used the Watts–Strogatz model to generate random networks with the small world properties appropriate for infectious disease transmission in these cities [16], [17].

The epidemic curves are all fitted very well using the small-world network structure models, indicating that the typical small-world property is able to capture the contact patterns during COVID-19 epidemics. The differences in these fitted parameters and starting times reflect the differences in the underlying transmission mechanisms and potential spread in the regions. The model then projected the trends of COVID-19 spread by simulating epidemics in the Wuhan, Toronto, and Italy networks. Our findings can guide public health authorities to implement effective mitigation strategies and be prepared for potential future outbreaks.

The network model

We develop a network-based model by extending the network SIR model [18] by incorporating the characteristics of COVID-19 transmission to assess the spread of the disease in heterogeneous populations. We derive the explicit expression of the epidemic threshold and discuss the final epidemic size for the network model.

Model formulation

We classify individuals by their average number of contacts in a typical day (time unit for the modeling) represented on the network by their degree . Individuals with degree are divided into susceptible (), exposed (), asymptomatically infected (), symptomatically infected (), hospitalized (), recovered (), and dead () states. Our model is formulated as follows where , , , , , , and represent the fractions of susceptible, exposed, asymptomatically infected, symptomatically infected, hospitalized, recovered, and dead individuals with degree , respectively. Here, is the total number of individuals with degree , and , and . represents the probability that an edge from a node with degree connects to a node with degree . For uncorrelated networks,  [19]. Since the node with degree shares an edge with the node degree , and only has free edges, a fraction of nodes may acquire the infection.

We assume that the transmission rates of symptomatically infected individuals and asymptomatically infected individuals are and , respectively. The factor accounts for the different transmission rates between asymptomatically infected individuals and symptomatically infected individuals. represents the fraction of nodes with degree infected by symptomatically infected nodes, and = represents the fraction of nodes with degree infected by asymptomatically infected nodes. Here, represents the probability that an edge from a degree node connects to a symptomatically infected node with degree , and represents the probability that an edge from a degree node connects to an asymptomatically infected node with degree .

In Model (2.1), the term represents the probability of not being infected by a symptomatically infected node with degree , and represents the probability of not being infected by an asymptomatically infected node with degree . Thus, is the probability that a node will neither be infected by a symptomatically infected nor be infected by an asymptomatically infected neighbor with degree , and is the probability of being infected by a symptomatically infected or an asymptomatically infected neighbor with degree .

Therefore, the susceptible individuals are infected at rate and enter the exposed state. After incubation period with a mean time of days, exposed individuals become symptomatically infected and asymptomatically infected with probabilities and , respectively. Symptomatically infected individuals are hospitalized at rate , and die at rate . Asymptomatically infected individuals, symptomatically infected individuals, and hospitalized individuals recover at rates , , and , respectively. Both the hospitalized individuals and symptomatically infected individuals die at rate .

Mathematical analysis

We derive the epidemic threshold to predict whether the epidemic will spread or die out and derive final epidemic size to quantify the total number of infected individuals.

The epidemic threshold

To estimate the transmission potential of the epidemic, we derive the important epidemic threshold, , defined as the average number of secondary cases produced by an infected individual in a completely susceptible population [20]. There exists a disease-free equilibrium, We compute following the next generation matrix approach presented by van den Driessche and Watmough [21]. For simplicity, we only consider the compartments related to infection, namely, , and , and rewrite the equations as the difference between vectors and following the notations in [21] where Here, represents the rate at which new infections are produced and represents the rate at which individuals transfer between compartments, and for Model (2.1).

The Jacobian matrix is where and

The matrices and are where is the identity matrix, and The next generation matrix is Since the rank of matrix is , the spectral radius of is its trace, i.e.,

It follows from (2.2) that the basic reproduction number becomes where and represent the average numbers of secondary cases produced by an asymptomatically infected individual and a symptomatically infected individual in a homogeneously mixed population, respectively. The term represents the average excess degree of nodes in the network [22].

Final epidemic size

We shall derive the final size following the approach in [23]. The nonlinear term in the first and second equations of Model (2.1) can be rewritten as When and , . Hence, Hence, Model (2.1) can be simplified as

We first claim that . Otherwise, there exist such that , since . Since , . Therefore, , leading to a contradiction. Similarly, we can show that

For a homogeneous network where all nodes have identical degree , Model (2.3) can be reduced to the following model with the initial values , , , , , and .

By Model (2.4) and a direct calculation, we have

By the first equation in Model (2.4), we further have where . To determine the final size of susceptible individuals, , we set where . By (2.5), (2.6) and the definition of , we have

It is clear that is a positive, increasing, strictly convex function, and . Thus, has a unique fixed point in the interval , which can be calculated numerically by using the iteration method and where denotes composition of for times. Then, the final size of susceptible individuals for a homogeneous network, , can be determined by .

We now derive the final size for heterogeneous networks. Integrating the first equation in Model (2.3) from to , we have

By summing and integrating the equations in Model (2.3), and where . We set

By Eqs. (2.7),(2.8), (2.9), we have where . Therefore, for all , the final size of susceptible individuals satisfies where

We define a map , by

To analyze the properties of , we shall introduce some notations. For , , we denote

Moreover, we shall claim if and . The above definition defines a partial order in . For later use, we could extend this partial order to matrices as follows. For any matrices , we have

When and , by the definition of and partial order defined in (2.10), we have

Since each is a increasing function, we have where is the composition function of for times. By the monotone criterion, we obtain

Due to the continuity of , and . Therefore, we have the following property [23].

All the fixed points of in the interval are contained in .

Due to the continuous differentiability of , for any and . Moreover, we shall simply write (2.11) in terms of the matrix form by where and .

By the monotony of , is also monotonous, i.e.,  for any . By utilizing the properties of and , we can obtain the following theorem.

Assume that the network is connected, we have

if and only if ;

when , has a unique fixed point satisfying .

The proof of Theorem 2.1 directly follows the proof in [23]. Hence, the final size of susceptible individuals for a heterogeneous network, , can be determined by to quantify the number of susceptible individuals left theoretically.

Parameter estimation and model based forecasting

We parameterized the model with reported data on COVID-19 cases and presented forecasts of the epidemic trends for the three areas.

Fitting reported confirmed cases

We simulated the spread of COVID-19 in Wuhan, Toronto, and Italy on the Watts–Strogatz network with degree and .

The study period for Wuhan starts from January 11, 2020, after the confirmed cases were reported, the public becomes aware of the infection and most people are trying to avoid gathering. The study period starts from January 26 for Toronto and from January 31 for Italy. In Toronto and Italy, usually people do not gather, especially after lockdown on Wuhan city, the awareness of avoiding exposure to the virus is increasing. Most people stay home during the study period, and the family sizes in Wuhan, Toronto, and Italy on average are all around . Therefore, the range of the node degrees is assumed to be between and .

The Watts–Strogatz model starts with a ring of vertices in which each vertex is connected to its nearest neighbors ( vertices clockwise and counterclockwise). Each edge is connected to a clockwise neighbor with probability and preserved with probability  [19], where the degree distribution is When , the expression reduces to a Poisson distribution as follows In the simulations, we used this degree distribution.

The total number of nodes for Wuhan, Toronto, and Italy are , , and as shown in Table 1, Table 3, and Table 6, respectively. We parameterized the model using the MCMC approach [24] by MATLAB R2016a according to the number of newly confirmed cases and the cumulative number of cases reported by the Health Commission of Hubei Province [25] and WHO [26].

Table 1

The parameter values and initial condition of four-phase simulations for Wuhan.

Parameter	Mean value	Std	95% CI	References
β (Phase 1)	0.04644	2.83×10^−3	[0.0409, 0.0520]	MCMC
β (Phase 2)	0.01597	4.79×10^−3	[6.59×10^−3, 0.0254]	MCMC
β (Phase 3)	2.90×10^−4	2.24×10^−4	[0, 7.29×10^−4]	MCMC
β (Phase 4)	7.30×10^−5	7.27×10^−5	[0, 2.16×10^−4]	MCMC
ξ (Phase 1)	0.8550	0.1220	[0.6159, 1.0942]	MCMC
ξ (Phase 2)	0.1499	0.1557	[0, 0.4552]	MCMC
ξ (Phase 3)	0.3369	0.1089	[0.1234, 0.5504]	MCMC
ξ (Phase 4)	0.8764	0.0836	[0.7124, 1.0403]	MCMC
∑_kN_k(0)	11081000	–	–	[27]
∑_kS_k(0)	11080770	–	–	Calculated
∑_kE_k(0)	200.26	18.88	[163.25, 237.26]	MCMC
∑_kA_k(0)	17.81	4.13	[9.71, 25.90]	MCMC
∑_kI_k(0)	11.78	6.06	[0, 23.66]	MCMC
∑_kH_k(0)	41	–	–	[25]
∑_kR_k(0)	0	–	–	Estimated
∑_kD_k(0)	0	–	–	Estimated

Phase 1 is between Jan 11 and Jan 23, Phase 2 is between Jan 23 and Feb 1, Phase 3 is between Feb 1 and Feb 12, and Phase 4 is between Feb 12 and Mar 31.

Table 3

The parameter values and initial condition of simulations for Toronto.

Parameter	Mean value	Std	95% CI	References
β (Phase 1)	7.95×10^−3	1.12×10^−3	[5.76×10^−3,0.01]	MCMC
β (Phase 2)	1.38×10^−4	1.79×10^−4	[0,4.90×10^−4]	MCMC
ξ (Phase 1)	0.1421	0.0734	[0, 0.2858]	MCMC
ξ (Phase 2)	0.1140	0.0709	[0, 0.2530]	MCMC
∑_kN_k(0)	5928000	–	–	[30]
∑_kS_k(0)	5927990	–	–	Calculated
∑_kE_k(0)	4.12	4.83	[0, 13.58]	MCMC
∑_kA_k(0)	2.43	2.22	[0, 6.78]	MCMC
∑_kI_k(0)	3.22	2.89	[0, 8.88]	MCMC
∑_kH_k(0)	1	–	–	[26]
∑_kR_k(0)	0	–	–	Estimated
∑_kD_k(0)	0	–	–	Estimated

Phase 1 is from Jan 26 to Mar 18, and Phase 2 is from Mar 18 to Mar 29.

Table 6

The parameter values and initial condition of simulations for Italy.

Parameter	Mean value	Standard derivation	95% CI	References
β (Phase 1)	0.0179	1.19×10^−3	[0.0156, 0.0203]]	MCMC
β (Phase 2)	4.45×10^−4	2.30×10^−4	[0,8.95×10^−4]	MCMC
ξ (Phase 1)	0.0996	0.0458	[9.87×10^−3,0.1894]	MCMC
ξ (Phase 2)	0.1312	0.0250	[0.0823, 0.1801]	MCMC
∑_kN_k(0)	59430000	–	–	[31]
∑_kS_k(0)	59429892	–	–	Calculated
∑_kE_k(0)	69.01	56.92	[0, 180.58]	MCMC
∑_kA_k(0)	22.60	18.24	[0, 58.35]	MCMC
∑_kI_k(0)	32.11	24.23	[0, 79.59]	MCMC
∑_kH_k(0)	2	–	–	[26]
∑_kR_k(0)	0	–	–	Estimated
∑_kD_k(0)	0	–	–	Estimated

Phase 1 is between Jan 31 and Mar 8, and Phase 2 is between Mar 8 and Mar 26.

The rate at which the fraction of the cumulative number of cases changes is , where represents the fraction of the cumulative number of infected individuals with degree . The number of newly infected can be expressed as where represents the number of new cases with degree , and represents the total number of individuals with degree . We run the MCMC simulation for iterations to fit the value of .

Zhou et al. showed that the median time from illness onset (i.e., before admission) to discharge was 22 days (IQR 18–25), whereas the median time to death was 18.5 days with IQR between and days [28]. We assume an exponential distribution for the time to recovery for asymptomatically infected individuals, symptomatically infected individuals, and hospitalized individuals. This results in the recovery rates per day, and the mortality rate, is per day. The incubation period of COVID-19 is around days [4], resulting in the progression rate . Qiu et al. reported that around of people infected with COVID-19 are asymptomatic or only have mild symptoms, and their transmissibility is lower, but still significant [29]. Thus, we assume that the probability that an infected individual is asymptomatic is , and for simulations.

We divided the Wuhan epidemic into four phases according to the reported data [3]. The first phase is before lockdown on Jan 23, 2020. The second phase is between Jan 24, 2020 and Feb 1, 2020 when the hospitals were short of beds. The third phase is between Feb 2, 2020 and Feb 6, 2020 when the Thunder God Mountain Hospital (TGMH) and Fire God Mountain Hospital (FGMH) were put into use. The fourth phase began when door-to-door screening was implemented on Feb 7, 2020 and TGMH, FGMH, and Mobile Cabin Hospitals (MCH) were put into use.

The study period for Toronto (Canada) was decomposed into two phases, namely, the period before Mar 18 and the period after Mar 18 when the city announced the emergence and schools and universities in Toronto were closed on Mar 18.

The study period for Italy was divided into two phases. The early epidemic phase was between Jan 31, 2020 and Mar 8, 2020 when the infection was spreading through the northern provinces. The second period begins on Mar 9, 2020 when the national lockdown started.

Predicting future epidemic trends

The parameters and initial conditions of simulations for Wuhan on the WS network are shown in Table 1. The probability of transmission through adequate contact is estimated by MCMC. The realizations of the basic reproduction numbers derived for Wuhan using the parameter values listed in Table 1 are shown in Table 2.

Table 2

Basic reproduction numbers computed by MCMC on the WS network.

Location	Period	Mean value	Standard derivation	95% CI
Wuhan	Jan 11–Jan 23	3.4074	0.2099	[2.9959, 3.8188]
	Jan 23–Feb 1	1.3065	0.3976	[0.5273, 2.0858]
	Feb 1–Feb 12	0.0221	0.0170	[0, 0.0555]
	Feb 12–Mar 31	5.35×10^−3	5.34×10^−3	[0, 0.0158]
Toronto	Jan 26–Mar 18	0.6416	0.0867	[0.4716, 0.8116]
	Mar 18–Mar 29	0.0115	0.0151	[0, 0.0412]
Italy	Jan 31–Mar 8	1.4763	0.0984	[1.2834, 1.6691]
	Mar 8–Mar 26	0.0359	0.0185	[0, 0.0721]

From Jan 11 to Mar 31, we estimate that the mean reproduction number on the WS network decreases from in the first phase to in the fourth phase. The epidemic on the WS network is shown in Fig. 1. Up to Jan 23, 2020 when Wuhan lockdown started, the estimated epidemic size is . During the second stage, after the lockdown of Wuhan and before the TGMH and FGMH were put into use, the predicted final size is . Thus, the lockdown of Wuhan reduced the expected final size by . During the third stage, after TGMH and FMGH were put into use, the final size is . Hence, the city lockdown and the usage of TGMH and FGMH reduced the final size by . During the fourth stage, after MCH was put into use, the predicted final size is , and the expected final size of infection is reduced by due to the increase of healthcare capacity.

Fig. 1

Fitting the number of reported new cases and the cumulative number of reported cases between Jan 11, 2020 and Mar 31, 2020 for Wuhan on Watts–Strogatz network.(A) Fitting the number of reported new cases on the Watts–Strogatz network. (B) Fitting the cumulative number of reported cases on the Watts–Strogatz network.

The variability of the numbers of confirmed new cases is consistent with the variability of the reproduction numbers listed in Table 2. In the first two phases, the epidemic spread rapidly with larger reproduction numbers that are larger than one, and the numbers of infected cases increase. In the last two phases, the spread is contained as the reproduction numbers falling below . In the third phase, because a large number of cases are confirmed by door-to-door screening and expanded healthcare capacity, the cumulative number of confirmed cases increased. On the other hand, the epidemic will die out because the reproduction number is less than one. In the fourth phase, the epidemic has been under control with the reproduction number being less than one. Hence, the number of new cases decreases.

The parameters and initial conditions of simulations for the GTA are shown in Table 3. The realizations of the basic reproduction numbers derived for Toronto using the parameter values listed in Table 3 are shown in Table 2. The reproduction numbers are much smaller due to social distancing policy, school closure, as well as behavior changes. The summary of the simulations is shown in Table 4, Table 5 and simulation results are shown in Fig. 2. The peak size is 60.19 (95% CI: 47.42–72.97), the peak time is Apr 2 (95% CI: Mar 29-Apr 7), and the final size is (95% CI: 1603–3820).

Table 4

The peak number of new cases, peak time, and final epidemic size after containment strategies are implemented in Toronto.

Scenarios	Peak size (95% CI)	Peak time (95% CI)	Final size (95% CI)
β	60.19 (47.42, 72.97)	Apr 2 (Mar 29, Apr 7)	2712 (1603, 3820)
0.8β	50.30 (41.03, 59.57)	Apr 1 (Mar 28, Apr 6)	2217 (1451, 2984)
0.6β	40.94 (34.55, 47.32)	Mar 31 (Mar 26, Apr 4)	1751 (1239, 2262)
q=0	60.19 (47.42, 72.97)	Apr 2 (Mar 29, Apr 7)	2712 (1603, 3820)
q=1∕8	44.66 (39.12, 50.21)	Mar 28 (Mar 26, Mar 31)	1486 (1111, 1861)
q=1∕4	38.79 (33.86, 43.72)	Mar 27 (Mar 25, Mar 29)	1133 (878, 1388)

Table 5

The peak number of new cases, peak time, and final epidemic size for Toronto when varying node degrees of symptomatically infected individuals.

Degree	Peak size (95% CI)	Peak time (95% CI)	Final size (95% CI)
k	60.19 (47.42, 72.97)	Apr 2 (Mar 29, Apr 7)	2712 (1603, 3820)
k−1	51.92 (42.08, 61.77)	Apr 1 (Mar 28, Apr 6)	2301 (1476, 3125)
k−2	43.98 (36.76, 51.20)	Mar 31 (Mar 27, Apr 5)	1908 (1310, 2506)
k−3	36.61 (30.96, 42.26)	Mar 30 (Mar 25, Apr 3)	1531 (1114, 1949)

Fig. 2

Fitting the number of reported new cases and the cumulative number of reported cases for Toronto on the Watts–Strogatz network. (A) Fitting the number of reported new cases on the Watts–Strogatz network. (B) Fitting the cumulative number of reported cases on the Watts–Strogatz network.

The parameters and initial condition of simulations for Italy is shown in Table 6. The realizations of the basic reproduction numbers derived for Italy using the parameter values listed in Table 6 are shown in Table 2. The reproduction numbers in the second phase are much smaller than that in the first phase due to the awareness of the severity of the epidemic. The summary of the simulation results is shown in Table 7, Table 8. Fig. 3 shows that the peak number of new cases is (95% CI: 5277–5708) on Mar 26 (95% CI: Mar 24-Mar 27), and the final size is (95% CI: ).

Table 7

The peak number of new cases, peak time, and final size after containment strategies are implemented in Italy.

Scenarios	Peak size (95% CI)	Peak time (95% CI)	Final size (95% CI)
β	5492 (5277, 5708)	Mar 26 (24, 27)	2.59×10^5(2.10×10^5,3.08×10^5)
0.8β	4340 (4235, 4564)	Mar 25 (24, 26)	2.03×10^5(1.75×10^5,2.31×10^5)
0.6β	3323 (3197, 3450)	Mar 25 (24, 26)	1.52×10^5(1.37×10^5,1.66×10^5)
q=0	5492 (5277, 5708)	Mar 26 (24, 27)	2.59×10^5(2.10×10^5,3.08×10^5)
q=1∕8	3413 (3279, 3547)	Mar 20 (19, 21)	1.22×10^5(1.14×10^5,1.30×10^5)
q=1∕4	2609 (2434, 2783)	Mar 18 (17, 19)	8.51×10^4(7.91×10^4,9.11×10^4)

Table 8

The peak size, peak time and final size for Italy when varying node degrees of symptomatically infected individuals.

Degree	Peak size (95% CI)	Peak time (95% CI)	Final size (95% CI)
k	5492 (5277, 5708)	Mar 26 (24, 27)	2.59×10^5(2.10×10^5,3.08×10^5)
k−1	4583 (4410, 4755)	Mar 25 (24, 27)	2.13×10^5(1.81×10^5,2.44×10^5)
k−2	3683 (3546, 3821)	Mar 25 (24, 26)	1.69×10^5(1.50×10^5,1.88×10^5)
k−3	2795 (2680, 2910)	Mar 24 (23, 26)	1.28×10^5(1.17×10^5,1.39×10^5)

Fig. 3

Fitting the number of reported new cases and the cumulative number of reported cases for Italy on the Watts–Strogatz network. (A) Fitting the number of reported new cases on the Watts–Strogatz network. (B) Fitting the cumulative number of reported cases on the Watts–Strogatz network.

The impact of mitigation strategies

The close contacts identified by contact tracing will be quarantined due to exposure to COVID-19 to see if they become sick. To evaluate the impact of mitigation strategies on the spread of COVID-19, Model (2.1) is rewritten as follows where represents the fraction of quarantined individuals with degree . The parameter represents the rate at which susceptible individuals are quarantined, and represents the rate at which the quarantined and uninfected close contacts transfer to the susceptible compartment again. In the simulations, we let to approximate a mean time of 14 days in the exposed state.

For Wuhan, the cumulative number of infected individuals after lockdown and TGMH, FGMH, as well as MCH were put into use are shown in Fig. 4. The results show that the lockdown and the increase in healthcare capacity are effective in controlling the numbers of confirmed cases.

Fig. 4

The impact of the variability on the healthcare capacity on the spread of the epidemic in Wuhan on the Watts–Strogatz network.

For Toronto, the number of newly infected individuals and the cumulative number of infected individuals produced on the WS network after implementing additional containment strategies besides school closure are shown in Fig. 5. We simulated the scenarios of implementing various containment strategies for Toronto. Simulation results showed that personal protection, reducing the node degrees of symptomatically infected individuals, and quarantine of close contacts are effective in reducing the peak epidemic size and final epidemic size. When is reduced by by personal protection or social distancing, the peak occurs one day earlier, and the final epidemic size is reduced by around . When is reduced by , the peak occurs two days earlier, and the final epidemic size is reduced by around . When , the peak occurs four days earlier, and the final epidemic size is reduced by . When , the peak appears five days earlier, and the final epidemic size is reduced by . When the node degrees of symptomatically infected individuals are reduced by , , and , the number of new cases produced per day at the peak is reduced by , , and . The final epidemic size is reduced by , , and when the node degrees of symptomatically infected individuals are reduced by , , and , respectively.

Fig. 5

The impact of mitigation strategies on the spread of COVID-19 epidemic in Toronto on the Watts–Strogatz network. In this figure and the following figure, the dashed lines represent confidence intervals. In (A) and (B), the red, purple, and green lines represent that the transmission rates are unchanged, reduced by , and reduced by , respectively. In (C) and (D), the red, purple, and green lines represent the rate of quarantine, , 1/8, and 1/4, respectively. In (E) and (F), the red, purple, and green, and light blue lines represent that the node degrees of symptomatically infected individuals are reduced by , , , and , respectively. (A) The number of newly infected individuals after reducing the transmission rates by personal protection and social distancing. (B) The cumulative number of infected individuals after reducing the transmission rates by personal protection and social distancing. (C) The number of newly infected individuals after close contacts are quarantined. (D) The cumulative number of infected individuals after close contacts are quarantined. (E) The number of newly infected individuals after the node degrees of symptomatically infected individuals are reduced. (F) The cumulative number of infected individuals after the node degrees of symptomatically infected individuals are reduced.

For Italy, the number of newly infected individuals and the cumulative number of infected individuals simulated on the WS network after implementing hypothetical containment strategies are shown in Fig. 6. Various scenarios of implementing mitigation strategies showed that the peak epidemic size and final epidemic size in Italy are greatly reduced by personal protection, social distancing, behavior change of symptomatically infected individuals, and quarantine. The simulations show that the peak would have arrived earlier if the containment had been intensified.

Fig. 6

The impact of mitigation strategies on the spread of COVID-19 epidemic in Italy on the Watts–Strogatz network. (A) The number of newly infected individuals after reducing the transmission rates by personal protection and social distancing. (B) The cumulative number of infected individuals after reducing the transmission rates by personal protection and social distancing. (C) The number of newly infected individuals after close contacts are quarantined. (D) The cumulative number of infected individuals after close contacts are quarantined. (E) The number of newly infected individuals after the node degrees of symptomatically infected individuals are reduced. (F) The cumulative number of infected individuals after the node degrees of symptomatically infected individuals are reduced.

When the probability of contact transmission coefficient , is reduced by by personal protection or social distancing, the peak occurs one day earlier, and the final epidemic size is reduced by . When is reduced by , the peak occurs one day earlier, and the final epidemic size is reduced by around .

When , the peak occurs six days earlier, and the final epidemic size is reduced by . Yet, when , the peak occurs eight days earlier, and the final epidemic size is reduced by . When the node degree of symptomatically infected individuals is reduced by , , and , the number of new cases produced per day at the peak is reduced by , , and , respectively. The final epidemic size is reduced by , , and when the node degrees of symptomatically infected individuals are reduced by , , and , respectively.

Summary and discussions

Modeling the dynamics of COVID-19 epidemics and assessment of mitigation strategies could be instrumental to public health agencies for surveillance and healthcare planning. For the models to be reliable, the simulated epidemic must account for the stochastic and heterogeneous contact among individuals. Hence, we developed a network model that captured the contact heterogeneity among individuals. We applied the model to analyze the transmission potential, and mitigation strategies for curbing the spread of COVID-19 epidemics in the cities of Wuhan, China and Toronto, Canada, and in the Italian Republic. The epidemic threshold derived from our network model can be used to predict the risks of spreading scenarios. We also provided an explicit expression of the final epidemic size, which facilitates estimating the scale of an outbreak for any region of interest. Our results provide insights in defining a mathematical framework for the analysis and containment of epidemic transmission in the real world.

A wide range of mitigation strategies can be examined by the flexible model framework. It can be extended to quantify the effectiveness of personal protection, social distancing, reducing the node degree of infected individuals, and quarantine on the dynamics of epidemics in different regions. When the mitigation strategy is intensified, the model predicts that the number of new cases peaks earlier and the final epidemic size is greatly reduced.

The social contact network structure and parameter values determine the transmission and epidemic course of such an emerging infectious disease. We choose the Watts–Strogatz to approximate real social networks, when the exact contact tracing data is unavailable. We assumed that the range of the node degree is between one and ten for each network in the absence of real contact tracing data, that is, on average each day an infected person would have between one and ten contacts where they could transmit the infection to another person. In the real world, the range of the degree will depend on the distribution of the household sizes of the region and time being studied. Moreover, the network structure can be altered by behavior change of individuals during epidemics. When this happens, the network structure can be adapted in our model to predict the impact of these changes on the epidemic threshold, epidemic peak value, peak time, stopping time, and final size of infected population.

The epidemics for the three places under study were fitted very well by our model with a small confidence interval. Hence, the forecasts by the model can be reliable. We did not provide the stopping time since too many uncertainties may affect the duration of the epidemics. As shown in the simulations, the transmission dynamics for four phases in Wuhan are quite different due to the variability on the intensity of interventions, the availability of healthcare facilities, as well as the utilization of personal protective equipment (PPE). The dynamics in the first phase is quite different from that in the second phase for Toronto. The same phenomenon is observed in Italy.

At the early stage, almost no interventions were implemented, and the public was not aware of or did not pay much attention to the severity of the highly contagious disease. With the increase of the number of reported confirmed cases and with the aid of social media, the public becomes aware of the severe consequence and has increased the level of personal protection and have avoided gathering, so that the reproduction number decreases and the estimated epidemic size declines by reducing the node degree of the network. Similarly, after applying the mitigation measures in Italy on March 8 and closing all schools in Toronto on March 18, the epidemics tend to be under control.

Hence, social distancing, self-isolation, quarantine, the utilization of PPE, and other measures of avoiding exposure to the virus can greatly reduce the size of infection during the COVID-19 outbreak. Therefore, it is essential to raise the awareness of these countermeasures to avoid contact between individuals. The possibility of recurrent outbreaks of the disease cannot be overstated. Even if the number of new cases is declining, it is still necessary to continue taking protective measures to prevent the occurrence of future outbreaks. The social media should warn the public not to relax their vigilance against the contagion of such a highly infectious disease.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.