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 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 is 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 neighbours. 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.
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.
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 k (k = 1, 2, · · · , n). Individuals with degree k are divided into susceptible (S k ), exposed (E k ), asymptomatically infected (A k ), symptomatically infected (I k ), hospitalized (H k ), recovered (R k ), and dead (D k ) states. Our model is formulated as follows
and d k = D k /N k represent the fractions of susceptible, exposed, asymptomatically infected, symptomatically infected, hospitalized, recovered, and dead individuals with degree k, respectively. Here, N k is the total number of individuals with degree k, and N k = S k +E k +A k + I k +H k + R k +D k , and s k +e k +a k + i k +h k + r k +d k = 1. P (k ′ |k) represents the probability that an edge from a node with degree k connects to a node with degree k ′ . For uncorrelated networks, P (k ′ |k) = k ′ P (k ′ )/ k [19] . Since the node with degree k ′ shares an edge with the node degree k, and only has (k ′ − 1) free edges, a fraction k ′ −1 k ′ 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. βks k k ′ k ′ −1 k ′ P (k ′ |k)i k ′ represents the fraction of nodes with degree k infected by symptomatically infected nodes, and σβks k k ′
k ′ P (k ′ |k)σa k ′ represents the fraction of nodes with degree k infected by asymptomatically infected nodes. Here, k ′ −1 k ′ P (k ′ |k)i k ′ represents the probability that an edge from a degree k node connects to a symptomatically infected node with degree k ′ , and k ′ −1 k ′ P (k ′ |k)σa k ′ represents the probability that an edge from a degree k node connects to an asymptomatically infected node with degree k ′ .
In Model (2.1), the term (1 − k ′ −1 k ′ P (k ′ |k)i k ′ ) represents the probability of not being infected by a symptomatically infected node with degree k ′ , and (1 − k ′ −1 k ′ P (k ′ |k)σa k ′ ) represents the probability of not being infected by an asymptomatically infected node with degree k ′ . 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 k ′ , and
is the probability of being infected by a symptomatically infected or an asymptomatically infected neighbor with degree k ′ . Therefore, the susceptible individuals are infected at rate
and enter the exposed state. After incubation period with a mean time of 1/ǫ days, exposed individuals become symptomatically infected and asymptomatically infected with probabilities δ and 1−δ, respectively. Symptomatically infected individuals are hospitalized at rate ξ, and die at rate µ. Asymptomatically infected individuals, symptomatically infected individuals, and hospitalized individuals recover at rates γ a , γ, and γ h , respectively. Both the hospitalized individuals and symptomatically infected individuals die at rate µ.
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.
To estimate the transmission potential of the epidemic, we derive the important epidemic threshold, R 0 , 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, (s 1 , · · · , s n , e 1 , · · · , e n , a 1 , · · · , a n , i 1 , · · · , i n , h 1 , · · · , h n , r 1 , · · · , r n , d 1 , · · · , d n ) T = (1, · · · , 1, 0, · · · , 0, 0, · · · , 0, 0, · · · , 0, 0, · · · , 0, 0, · · · , 0, 0, · · · , 0) T =: E 0 .
We compute R 0 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, e k , a k and i k , and rewrite the equations as the difference between vectors F k and V k following the notations in [21] [
Here, F ik represents the rate at which new infections are produced and V ik represents the rate at which individuals transfer between compartments, i = 1, 2, 3 and k = 1, · · · , n for Model (2.1).
The Jacobian matrix F is
where z = (z j ) = (e 1 , · · · , e n , a 1 , · · · , a n , i 1 , · · · , i n ) and
The matrices V and V −1 are
where I n is the n × n identity matrix, and
The next generation matrix is
Since the rank of matrix F ′ is 1, the spectral radius of F ′ is its trace, i.e.,
It follows from (2.2) that the basic reproduction number R 0 becomes
where β 1−δ γa σ 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
Hence, Model (2.1) can be simplified as
For a homogeneous network where all nodes have identical degree k, Model (2.3) can be reduced to the following model
By Model (2.4) and a direct calculation, we have
(2.5)
By the first equation in Model (2.4), we further have
where s(+∞) = lim t→∞ s(t). To determine the final size of susceptible individuals, s(+∞), we set
where y 0 = s 0 + e 0 . By (2.5), (2.6) and the definition of f (x), we have
It is clear that f (x) is a positive, increasing, strictly convex function, and f (s 0 ) < s 0 . Thus, f has a unique fixed point s + in the interval (0, s 0 ), which can be calculated numerically by using the iteration method and
where f m denotes composition of f for m times. Then, the final size of susceptible individuals for a homogeneous network, s(+∞), can be determined by s + .
We now derive the final size for heterogeneous networks. Integrating the first equation in Model (2.3) from 0 to t, we have
By summing and integrating the equations in Model (2.3),
where y k (t) = s k (t) + e k (t). We set
By Equations (2.7),(2.8) and (2.9), we have
where g k (0) ≥ 0, ∀k. Therefore, for all k = 1, ..., n, the final size of susceptible individuals satisfies
We define a map G :
To analyze the properties of G(x), we shall introduce some notations. For
The above definition defines a partial order in R n . For later use, we could extend this partial order to n × n matrices as follows. For any n × n matrices A, B, we have
. . , w n (0)] T , by the definition of G(x) and partial order defined in (2.10), we have
where G m is the composition function of G for m times. By the monotone criterion, we obtain
Due to the continuity of G, G(s) = s and G(s) = s. Therefore, we have the following property [23] . Due to the continuous differentiability of G,
for any x ∈ R n and 1 ≤ i, j ≤ n. Moreover, we shall simply write (2.11) in terms of the matrix form by
. By the monotony of G, DG is also monotonous, i.e., DG(x) ≤ DG(y) for any 0 ≤ x ≤ y ≤ s(0). By utilizing the properties of w(x) and G(x), we can obtain the following theorem. (2) when w(0) > 0, G has a unique fixed point s ++ satisfying 0 ≪ s ++ < s(0).
The proof of Theorem 2.1 directly follows the proof in [23] . Hence, the final size of susceptible individuals for a heterogeneous network, s(+∞), can be determined by s ++ to quantify the number of susceptible individuals left theoretically.
We parameterized the model with reported data on COVID-19 cases and presented forecasts of the epidemic trends for the three areas.
We simulated the spread of COVID-19 in Wuhan, Toronto, and Italy on the Watts-Strogatz network with degree k min = 1 and k max = 10.
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 3. Therefore, the range of the node degrees is assumed to be between 1 and 10.
The Watts-Strogatz model starts with a ring of N vertices in which each vertex is connected to its 2m nearest neighbors (m vertices clockwise and m counterclockwise). Each edge is connected to a clockwise neighbor with probability p and preserved with probability 1 − p [19] , where the degree distribution is
When p → 1, 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 11081000, 5928000, and 59430000 as shown in Table 1 , Table 3 , and The rate at which the fraction of the cumulative number of cases changes is dc k /dt = ξi k , where c k (t) represents the fraction of the cumulative number of infected individuals with degree k. The number of newly infected can be expressed as
where P k represents the number of new cases with degree k, and N k represents the total number of individuals with degree k. We run the MCMC simulation for 20000 iterations to fit the value of P k . 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 15 and 22 days [27] . 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 γ a = γ = γ h = 1/22 per day, and the mortality rate, µ is 1/18.5 per day. The incubation period of COVID-19 is around 7 days [4] , resulting in the progression rate ǫ = 1/7. Qiu et al. reported that around 30% − 60% of people infected with COVID-19 are asymptomatic or only have mild symptoms, and their transmissibility is lower, but still significant [28] . Thus, we assume that the probability that an infected individual is asymptomatic is 1 − δ = 0.6, and σ = 1 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 doorto-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.
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 5000 realizations of the basic reproduction numbers derived for Wuhan using the parameter values listed in Table 1 are shown in Table 2 .
From Jan 11 to Mar 31, we estimate that the mean reproduction number on the WS network decreases from 3.41 in the first phase to 5.34 × 10 −3 in the fourth phase. The epidemic on the WS network is shown in Figure 1 . Up to Jan 23, 2020 when Wuhan lockdown started, the estimated epidemic size is 3.96 × 10 6 . During the second stage, after the lockdown of Wuhan and before the TGMH and FGMH were put into use, the predicted final size is 2.17 × 10 6 . Thus, the lockdown of Wuhan reduced the expected final size by 45.22%. During the third stage, after TGMH and FMGH were put into use, the final size is 1.02 × 10 5 . Hence, the city lockdown and the usage of TGMH and FGMH reduced the final size by 97.42%. During the fourth stage, after MCH was put into use, the predicted final size is 51269, and the expected final size of infection is reduced by 98.70% due to the increase of healthcare capacity.
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 1, and the numbers of infected cases increase. In the last two phases, the spread of disease is controlled, and the reproduction numbers are smaller than 1. 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 spread of the disease 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 5000 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 and Table 5 . Simulation results are shown in Figure 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 2712 (95%CI: 1603-3820).
The parameters and initial condition of simulations for Italy is shown in Table 6 . The 5000 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 and Table 8 . Figure  3 shows that the peak number of new cases is 5492 (95%CI: 5277-5708) on Mar 26 (95%CI: Mar 24-Mar 27), and the final size is 2.59 × 10 5 (95%CI: 2.10 × 10 5 − 3.08 × 10 5 ).
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 sq k = SQ k /N k represents the fraction of quarantined individuals with degree k. The parameter q 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 λ = 1/14 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 Figure 4 . The results show that the lockdown and the increase in healthcare capacity are effective in controlling the numbers of confirmed cases.
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 Figure 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. Reducing the transmission rate β, by x% also reduces R 0 by x%. When β is reduced by 20% by personal protection or social distancing, the peak occurs one day earlier, and the final epidemic size is reduced by around 18%. When β is reduced by 40%, the peak occurs two days earlier, and the final epidemic size is reduced by around 33.3%. When q = 1/8, the peak occurs four days earlier, and the final epidemic size is reduced by 45.21%. When q = 1/4, the peak appears five days earlier, and the final epidemic size is reduced by 58.22%. When the node degrees of symptomatically infected individuals are reduced by 1, 2, and 3, the number of new cases produced per day at the peak is reduced by 13.74%, 26.93%, and 39.18%. The final epidemic size is reduced by 15.15%, 29.65%, and 43.55% when the node degrees of symptomatically infected individuals are reduced by 1, 2, and 3, respectively.
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 Figure 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.
When the probability of contact transmission coefficient β, is reduced by 20% by personal protection or social distancing, the peak occurs one day earlier, and the final epidemic size is reduced by 21.56%. When β is reduced by 41.44%, the peak occurs one day earlier, and the final epidemic size is reduced by around 40%.
When q = 1/8, the peak occurs six days earlier, and the final epidemic size is reduced by 52.87%. Yet, when q = 1/4, the peak occurs eight days earlier, and the final epidemic size is reduced by 67.12%. When the node degree of symptomatically infected individuals is reduced by 1, 2, and 3, the number of new cases produced per day at the peak is reduced by 16.50%, 32.93%, and 49.11%, respectively. The final epidemic size is reduced by 17.90%, 34.70%, and 50.51% when the node degrees of symptomatically infected individuals are reduced by 1, 2, and 3, respectively.
Modelling 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.
The flexible network model framework can simulate 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. 