A Data-Driven Network Model for the Emerging COVID-19 Epidemics in Wuhan, Toronto and Italy
Ling Xue, Shuanglin Jing, Joel C. Miller, Wei Sun, Huafeng Li, Jose Guillermo Estrada-Franco, James M Hyman, Huaiping Zhu
aa r X i v : . [ q - b i o . P E ] M a y A Data-Driven Network Model for the EmergingCOVID-19 Epidemics in Wuhan, Toronto and Italy
Ling Xue a , , Shuanglin Jing a , Joel C. Miller b , Wei Sun a , , Huafeng Li a ,Jos´e Guillermo Estrada-Franco c, ∗ , James M Hyman d , Huaiping Zhu e, ∗ a College of Mathematical Sciences, Harbin Engineering University, Harbin, Heilongjiang,150001, China b School of Engineering and Mathematical Sciences, Melbourne, La Trobe University,3086, Australia c Instituto Polit´ecnico Nacional, Centro de Biotecnolog´ıa Gen´omica, Cd. Reynosa,Tamaulipas, 88710, M´exico d Department of Mathematics, Tulane University, New Orleans, LA, 70118, USA e Lamps and Center of Disease Modelling (CDM), Department of Mathematics andStatistics, York University, Toronto, ON, M3J 1P3, Canada
Abstract
The ongoing Coronavirus Disease 2019 (COVID-19) pandemic threatens thehealth of humans and causes great economic losses. Predictive modellingand forecasting the epidemic trends are essential for developing counter-measures to mitigate this pandemic. We develop a network model, whereeach node represents an individual and the edges represent contacts betweenindividuals where the infection can spread. The individuals are classifiedbased on the number of contacts they have each day (their node degrees)and their infection status. The transmission network model was respectivelyfitted to the reported data for the COVID-19 epidemic in Wuhan (China),Toronto (Canada), and the Italian Republic using a Markov Chain MonteCarlo (MCMC) optimization algorithm. Our model fits all three regions wellwith narrow confidence intervals and could be adapted to simulate othermegacities or regions. The model projections on the role of containmentstrategies can help inform public health authorities to plan control measures. ∗ Corresponding author, Contributed equally
Email addresses: [email protected] (Jos´e Guillermo Estrada-Franco), [email protected] (Huaiping Zhu)
Preprint submitted to Elsevier June 1, 2020 eywords:
COVID-19, mitigation strategies, network model,heterogeneity, control measures
1. Introduction
The development of international trade and tourism has accelerated thespatial spread of infectious diseases. The limited data available on emerg-ing epidemics adds to the challenge of mitigating the spread of emerginginfections [1]. The unprecedented Coronavirus Disease 2019 (COVID-19)outbreak began at the end of 2019. The number of reported cases keepsrising worldwide and thousands of lives have been claimed. This pandemicis having an enormous impact on world health, disturbing the stability ofthe societies, and triggers great economic losses. Predicting the future of thepandemic, assessing the impacts of current interventions, and evaluating theeffectiveness of alternate mitigation strategies are of utmost importance forsaving lives.Mathematical models can be used to understand the dynamics of epi-demics and help inform control strategies. A numerous number of modelsare being used to project the current COVID-19 pandemic. Ziff and Ziff an-alyzed the number of reported cases for Wuhan (China) and showed that thegrowth of the daily number of confirmed new cases indicates an underlyingfractal or small-world network of connections between susceptible and in-fected individuals [2]. Wang et al. developed an SEIR model to estimate theepidemic trends in Wuhan, assuming the prevention and control measureswere either sufficient or insufficient to control the epidemic [3]. Kucharski etal. combined a stochastic transmission model with data on cases of COVID-19 in Wuhan and international cases to estimate how the transmission hadvaried over time between January and February in 2020 [4].Kraemer et al. analyzed the impact of interventions on the spread ofCOVID-19 in China using transportation data [5]. Chinazzi et al. used aglobal meta-population disease transmission model to project the impact oftravel limitations on the national and international spread of the epidemic.They showed that the travel restriction of Wuhan, China had a more markedeffect on the international scale than that on Mainland China [6].Ferguson et al. found that optimal mitigation policies (combining homeisolation of suspected cases, home quarantine of those living in the samehousehold as suspected cases, and social distancing of the elderly and others2t most risk of severe disease) might reduce peak healthcare demand by 2/3and deaths by half [7]. Likewise, Hellewell et al. developed a stochastictransmission model and found that highly effective contact tracing and caseisolation is enough to control a new outbreak of COVID-19 within threemonths in most scenarios [8].Zhang et al. fitted the reported serial interval (mean and standard devia-tion) with a gamma distribution to estimate the basic reproduction numberat the early stage of a COVID-19 outbreak, indicating the potential of secondoutbreaks [9]. Maier et al. developed a compartmental model dividing indi-viduals into susceptible, exposed, removed, and quarantined symptomaticallyinfected and showed that the distinctive subexponential increase of confirmedcases in mainland China could be explained as a direct consequence of con-tainment policies that effectively deplete the susceptible population [10].Most of these models are based on assuming the population is homoge-neously mixing, that is, the contacts between people are random and uni-formly distributed throughout the population. However, different individualsmay have varying numbers of acquaintances and contacts in the real world.The important role that heterogenous contact networks play in the trans-mission dynamics of infectious diseases is often overlooked [11]. Models thattake into account contact heterogeneity better represent the actual transmis-sion network through a population and are more likely to capture the trueepidemic dynamics.Disease propagation is closely linked with the structure of social contactnetworks [12]. The ubiquitous diversity in contact patterns and heterogene-ity among individuals depends on differences in social structures, spatialdistances, and behavior [13]. The heterogeneity exists at a wide range ofscales and leads to highly variable transmission dynamics of infectious dis-eases [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, mostnodes can be reached from every other node by a small number of hops orsteps, even if they are not immediate neighbours. This type of network modelallows us to adapt changes to some realistic network structures and examinethe effects of control and intervention countermeasures such as social dis-tancing, self-isolation, and personal protection. The framework and analysiscan be applied to study the transmission dynamics in different regions andmany other infectious diseases.The COVID-19 epidemic in Wuhan ended in April, while the epidemics3n the Greater Toronto Area (GTA, Canada) and the Italian Republic arecontinuing to grow. We fit the parameters of our network model to theconfirmed cases in each of these regions. Although Wuhan, Toronto, andItaly differ in some ways, the way that SARS-CoV-2 is transmitted fromone person to another is quite similar. Individuals may acquire infectionfrom other infectious individuals, even if they do not contact each otherdirectly. The Watts-Strogatz model supplies an ideal tool to study the spreadof epidemics among individuals even if their locations are not considered. Weused the Watts-Strogatz model to generate random networks with the smallworld properties appropriate for infectious disease transmission in these cities[16, 17].The epidemic curves are all fitted very well using the small-world networkstructure models, indicating that the typical small-world property is able tocapture the contact patterns during COVID-19 epidemics. The differencesin these fitted parameters and starting times reflect the differences in theunderlying transmission mechanisms and potential spread in the regions.The model then projected the trends of COVID-19 spread by simulatingepidemics in the Wuhan, Toronto, and Italy networks. Our findings canguide public health authorities to implement effective mitigation strategiesand be prepared for potential future outbreaks.
2. 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 assessthe spread of the disease in heterogeneous populations. We derive the explicitexpression of the epidemic threshold and discuss the final epidemic size forthe network model.
We classify individuals by their average number of contacts in a typi-cal day (time unit for the modeling) represented on the network by theirdegree k ( k = 1 , , · · · , n ). Individuals with degree k are divided into suscep-tible ( S k ), exposed ( E k ), asymptomatically infected ( A k ), symptomaticallyinfected ( I k ), hospitalized ( H k ), recovered ( R k ), and dead ( D k ) states. Our4odel is formulated as follows ds k dt = − βks k P k ′ [1 − (1 − k ′ − k ′ P ( k ′ | k ) i k ′ )(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ )] , de k dt = βks k P k ′ [1 − (1 − k ′ − k ′ P ( k ′ | k ) i k ′ )(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ )] − ǫe k , da k dt = (1 − δ ) ǫe k − γ a a k , di k dt = δǫe k − γi k − µi k − ξi k , dh k dt = ξi k − γ h h k − µh k , dr k dt = γi k + γ a a k + γ h h k , dd k dt = µi k + µh k , (2.1)where s k = S k /N k , e k = E k /N k , a k = A k /N k , i k = I k /N k , h k = H k /N k , r k = R k /N k , and d k = D k /N k represent the fractions of susceptible, ex-posed, asymptomatically infected, symptomatically infected, hospitalized,recovered, and dead individuals with degree k , respectively. Here, N k isthe 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 ) representsthe probability that an edge from a node with degree k connects to a nodewith degree k ′ . For uncorrelated networks, P ( k ′ | k ) = k ′ P ( k ′ ) / h k i [19]. Sincethe node with degree k ′ shares an edge with the node degree k , and only has( k ′ −
1) free edges, a fraction k ′ − k ′ of nodes may acquire the infection.We assume that the transmission rates of symptomatically infected in-dividuals and asymptomatically infected individuals are β and σβ , respec-tively. The factor σ accounts for the different transmission rates betweenasymptomatically infected individuals and symptomatically infected individ-uals. βks k P k ′ k ′ − k ′ P ( k ′ | k ) i k ′ represents the fraction of nodes with degree k infected by symptomatically infected nodes, and σβks k P k ′ k ′ − k ′ P ( k ′ | k ) a k ′ = βks k P k ′ k ′ − k ′ P ( k ′ | k ) σa k ′ represents the fraction of nodes with degree k infected by asymptomatically infected nodes. Here, k ′ − k ′ P ( k ′ | k ) i k ′ representsthe probability that an edge from a degree k node connects to a symptomat-ically infected node with degree k ′ , and k ′ − k ′ P ( k ′ | k ) σa k ′ represents the prob-ability that an edge from a degree k node connects to an asymptomaticallyinfected node with degree k ′ .In Model (2.1), the term (1 − k ′ − k ′ P ( k ′ | k ) i k ′ ) represents the probability ofnot being infected by a symptomatically infected node with degree k ′ , and(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ ) represents the probability of not being infected by anasymptomatically infected node with degree k ′ . Thus, (1 − k ′ − k ′ P ( k ′ | k ) i k ′ )(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ ) is the probability that a node will neither be infected by a5ymptomatically infected nor be infected by an asymptomatically infectedneighbor with degree k ′ , and 1 − (1 − k ′ − k ′ P ( k ′ | k ) i k ′ )(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ ) isthe probability of being infected by a symptomatically infected or an asymp-tomatically infected neighbor with degree k ′ .Therefore, the susceptible individuals are infected at rate βks k X k ′ [1 − (1 − k ′ − k ′ P ( k ′ | k ) i k ′ )(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ )]and enter the exposed state. After incubation period with a mean time of1 /ǫ days, exposed individuals become symptomatically infected and asymp-tomatically infected with probabilities δ and 1 − δ , respectively. Symptomati-cally infected individuals are hospitalized at rate ξ , and die at rate µ . Asymp-tomatically infected individuals, symptomatically infected individuals, andhospitalized individuals recover at rates γ a , γ , and γ h , respectively. Both thehospitalized individuals and symptomatically infected individuals die at rate µ . We derive the epidemic threshold to predict whether the epidemic willspread or die out and derive final epidemic size to quantify the total numberof infected individuals.
To estimate the transmission potential of the epidemic, we derive the im-portant epidemic threshold, R , defined as the average number of secondarycases produced by an infected individual in a completely susceptible popula-tion [20]. There exists a disease-free equilibrium,( s , · · · , s n , e , · · · , e n , a , · · · , a n , i , · · · , i n , h , · · · , h n , r , · · · , r n , d , · · · , d n ) T = (1 , · · · , , , · · · , , , · · · , , , · · · , , , · · · , , , · · · , , , · · · , T =: E . We compute R following the next generation matrix approach presentedby van den Driessche and Watmough [21]. For simplicity, we only considerthe compartments related to infection, namely, e k , a k and i k , and rewrite theequations as the difference between vectors F k and V k following the notationsin [21] [ de k dt , da k dt , di k dt ] T = F k − V k , F k = ( F ik ) = βks k P k ′ [1 − (1 − k ′ − k ′ P ( k ′ | k ) i k ′ )(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ )]00 , V k = ( V ik ) = ǫe k − (1 − δ ) ǫe k + γ a a k − δǫe k + γi k + µi k + ξi k . 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 , , k = 1 , · · · , n for Model (2.1).The Jacobian matrix F is F = (cid:20) ∂ F ik ∂z j (cid:21) E = n × n σβF ′ βF ′ n × n n × n n × n n × n n × n n × n , where z = ( z j ) = ( e , · · · , e n , a , · · · , a n , i , · · · , i n ) and F ′ = P (2 | · · · n − n P ( n | P (2 | · · · n − n P ( n | n P (2 | n ) · · · ( n − P ( n | n ) = 1 h k i P (2) · · · ( n − P ( n )0 2 P (2) · · · n − P ( n )... ... ... ...0 nP (2) · · · n ( n − P ( n ) . The matrices V and V − are V = (cid:20) ∂ V i ∂z j (cid:21) E = ǫ I n n × n n × n − (1 − δ ) ǫ I n γ a I n n × n − δǫ I n n × n ( γ + µ + ξ ) I n , where I n is the n × n identity matrix, and V − = ǫ I n n × n n × n − δγ a I n γ a I n n × nδγ + µ + ξ I n n × n γ + µ + ξ I n . F V − = β (cid:16) − δγ a σ + δγ + µ + ξ (cid:17) F ′ σβγ a F ′ βγ + µ + ξ F ′ n × n n × n n × n n × n n × n n × n . (2.2)Since the rank of matrix F ′ is 1, the spectral radius of F ′ is its trace, i.e., ρ ( F ′ ) = Tr( F ′ ) = 1 h k i X k ( k − kP ( k ) = h k i − h k ih k i . It follows from (2.2) that the basic reproduction number R becomes R = ρ ( F V − ) = β (cid:18) − δγ a σ + δγ + µ + ξ (cid:19) ( h k i − h k i ) h k i , where β − δγ a σ and β δγ + µ + ξ represent the average numbers of secondary casesproduced by an asymptomatically infected individual and a symptomaticallyinfected individual in a homogeneously mixed population, respectively. Theterm h k i−h k ih k i represents the average excess degree of nodes in the network[22]. We shall derive the final size following the approach in [23]. The nonlinearterm in the first and second equations of Model (2.1) can be rewritten as βks k X k ′ [1 − (1 − k ′ − k ′ P ( k ′ | k ) i k ′ )(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ )]= βks k X k ′ [1 − k ′ − k ′ P ( k ′ | k ) i k ′ + k ′ − k ′ P ( k ′ | k ) σa k ′ − ( k ′ − k ′ P ( k ′ | k )) i k ′ σa k ′ ]= βks k X k ′ [ k ′ − k ′ P ( k ′ | k ) i k ′ + k ′ − k ′ P ( k ′ | k ) σa k ′ − ( k ′ − k ′ P ( k ′ | k )) i k ′ σa k ′ ] . When i k ′ ≪ a k ′ ≪ i k ′ a k ′ ≈
0. Hence, βks k X k ′ [1 − (1 − k ′ − k ′ P ( k ′ | k ) i k ′ )(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ )]8 βks k X k ′ [ k ′ − k ′ P ( k ′ | k )( i k ′ + σa k ′ )] . Hence, Model (2.1) can be simplified as ds k dt = − βks k P k ′ k ′ − k ′ P ( k ′ | k ) ( i k ′ + σa k ′ ) . de k dt = βks k P k ′ k ′ − k ′ P ( k ′ | k ) ( i k ′ + σa k ′ ) − ǫe k , da k dt = (1 − δ ) ǫe k − γ a a k , di k dt = δǫe k − γi k − µi k − ξi k , dh k dt = ξi k − γ h h k − µh k , dr k dt = γi k + γ a a k + γ h h k , dd k dt = µi k + µh k . (2.3)We first claim that lim t → + ∞ e k ( t ) = 0. Otherwise, there exist T, η > e k ( t ) ≥ η > ∀ t > T since e k ( t ) ≥
0. Since ddt ( s k + e k ) = − ǫe k ≤− ǫη , s k ( t ) + e k ( t ) ≤ − ǫηt + s ( T ) + e ( T ). Therefore, lim t →∞ ( s k ( t ) + e k ( t )) = −∞ , leading to a contradiction. Similarly, we can show thatlim t → + ∞ e k ( t ) = lim t → + ∞ a k ( t ) = lim t → + ∞ i k ( t ) = lim t → + ∞ h k ( t ) = 0 , ∀ k. For a homogeneous network where all nodes have identical degree k ,Model (2.3) can be reduced to the following model dsdt = − βks h k − k i ( i + σa ) , dedt = βks h k − k i ( i + σa ) − ǫe, dadt = (1 − δ ) ǫe − γ a a, didt = δǫe − ( γ + µ + ξ ) i, dhdt = ξi − γ h h − µh, drdt = γi + γ a a + γ h h, dddt = µi + µh, (2.4)with the initial values s (0) = s , e (0) = e , a (0) = a , i (0) = i , h (0) = h , r (0) = 0 and d (0) = 0.By Model (2.4) and a direct calculation, we have − ( i + σa ) = 1 γ + µ + ξ (cid:18) dsdt + dedt + dadt + didt (cid:19) + (cid:18) σγ a − γ + µ + ξ (cid:19) (cid:18) (1 − δ ) (cid:18) dsdt + dedt (cid:19) + dadt (cid:19) . (2.5)9y the first equation in Model (2.4), we further have ln s (+ ∞ ) s = − βk h k − k i Z + ∞ ( i + σa ) dt, (2.6)where s (+ ∞ ) = lim t →∞ s ( t ). To determine the final size of susceptible indi-viduals, s (+ ∞ ), we set f ( x ) = s exp (cid:26) βk h k − k i (cid:20) x − y − a − i γ + µ + ξ + (cid:18) σγ a − γ + µ + ξ (cid:19) ((1 − δ ) ( x − y ) − a ) (cid:21)(cid:27) , where y = s + e . By (2.5), (2.6) and the definition of f ( x ), we have s (+ ∞ ) = f ( s (+ ∞ )) . It is clear that f ( x ) is a positive, increasing, strictly convex function, and f ( s ) < s . Thus, f has a unique fixed point s + in the interval (0 , s ), whichcan be calculated numerically by using the iteration method and s + = lim m → + ∞ f m ( s ) , where f m denotes composition of f for m times. Then, the final size of sus-ceptible individuals for a homogeneous network, s (+ ∞ ), can be determinedby s + .We now derive the final size for heterogeneous networks. Integrating thefirst equation in Model (2.3) from 0 to t , we have ln s k ( t ) s k (0) = − βk X k ′ k ′ − k ′ P ( k ′ | k ) Z t ( i k ′ + σa k ′ ) du. (2.7)By summing and integrating the equations in Model (2.3), Z t i k du = − δγ + µ + ξ ( y k ( t ) − y k (0)) − i k ( t ) − i k (0) γ + µ + ξ , (2.8)and Z t a k du = − − δγ a ( y k ( t ) − y k (0)) − γ a ( a k ( t ) − a k (0)) , (2.9)10here y k ( t ) = s k ( t ) + e k ( t ). We set g k ( t ) = (cid:18) δγ + µ + ξ + σ − δγ a (cid:19) y k ( t ) + σγ a a k ( t ) + i k ( t ) γ + µ + ξ . By Equations (2.7),(2.8) and (2.9), we have ln s k (+ ∞ ) s k (0) = βk n X j =1 j − j P ( j | k ) ( g j (+ ∞ ) − g j (0))= βk n X j =1 j − j P ( j | k ) (cid:18)(cid:18) δγ + µ + ξ + σ − δγ a (cid:19) s j (+ ∞ ) − g j (0) (cid:19) , where g k (0) ≥ , ∀ k . Therefore, for all k = 1 , ..., n , the final size of susceptibleindividuals satisfies s k (+ ∞ ) = s k (0) exp ( βk n X j =1 j − j P ( j | k ) ·· (cid:18) ( δγ + µ + ξ + σ − δγ a ) ( s j (+ ∞ ) − s j (0)) − w j (0) (cid:19)(cid:27) , where w k (0) = (cid:18) δγ + µ + ξ + σ − δγ a (cid:19) e k (0) + σγ a a k (0) + i k (0) γ + µ + ξ ≥ . We define a map G : R n → R n , x = ( x j ) G ( x ) = ( G ( x ) , ..., G n ( x )) T by G i ( x ) = s i (0) exp ( βi n X j =1 j − j P ( j | i ) ·· (cid:18) ( δγ + µ + ξ + σ − δγ a )( x j − s j (0)) − w j (0) (cid:19)(cid:27) . To analyze the properties of G ( x ), we shall introduce some notations. For Y = ( Y , . . . , Y n ) T , Z = ( Z , . . . , Z n ) T ∈ R n , we denote Y ≤ Z (resp. Y ≪ Z ) if Y l ≤ Z l (resp. Y l < Z l ) , ∀ l = 1 , . . . , n. (2.10)11oreover, we shall claim Y < Z if Y ≤ Z and Y = Z . The abovedefinition defines a partial order in R n . For later use, we could extend thispartial order to n × n matrices as follows. For any n × n matrices A, B , wehave A ≤ B if Ax ≤ Bx, ∀ ≤ x ∈ R n . When 0 ≪ s (0) = [ s (0) , . . . , s n (0)] T and 0 ≤ w (0) = [ w (0) , . . . , w n (0)] T ,by the definition of G ( x ) and partial order defined in (2.10), we have0 ≪ G (0) ≤ G ( s (0)) ≤ s (0) , Since each G i ( x ) is a increasing function, we have0 ≪ G (0) ≤ · · · ≤ G m (0) ≤ G m ( s (0)) ≤ · · · ≤ G ( s (0)) ≤ s (0) , where G m is the composition function of G for m times. By the monotonecriterion, we obtain0 ≪ s := lim m → + ∞ G m (0) ≤ s := lim m → + ∞ G m ( s (0)) ≤ s (0) . Due to the continuity of G , G ( s ) = s and G ( s ) = s . Therefore, we havethe following property [23]. Lemma 2.1.
All the fixed points of G in the interval [0 , s (0)] are containedin [ s, s ] . Due to the continuous differentiability of G , ∂G i ( x ) ∂x j = βi j − j P ( j | i ) (cid:18) δγ + µ + ξ + σ − δγ a (cid:19) G i ( x ) (2.11)for any x ∈ R n and 1 ≤ i, j ≤ n . Moreover, we shall simply write (2.11) interms of the matrix form by DG ( x ) := ∂G ( x ) ∂x = diag ( G ( x ) , · · · , G n ( x )) B, where B = [ b ij ] and b ij = βi j − j (cid:16) δγ + µ + ξ + σ − δγ a (cid:17) P ( j | i ).By the monotony of G , DG is also monotonous, i.e., DG ( x ) ≤ DG ( y ) forany 0 ≤ x ≤ y ≤ s (0). By utilizing the properties of w ( x ) and G ( x ), we canobtain the following theorem. 12 heorem 2.1. Assume that the network is connected, we have (1) w (0) = 0 if and only if G ( s (0)) = s (0) ; (2) when w (0) > , G has a unique fixed point s ++ satisfying ≪ s ++
3. Parameter estimation and model based forecasting
We parameterized the model with reported data on COVID-19 cases andpresented forecasts of the epidemic trends for the three areas.
We simulated the spread of COVID-19 in Wuhan, Toronto, and Italy onthe Watts-Strogatz network with degree k min = 1 and k max = 10.The study period for Wuhan starts from January 11, 2020, after theconfirmed cases were reported, the public becomes aware of the infectionand most people are trying to avoid gathering. The study period starts fromJanuary 26 for Toronto and from January 31 for Italy. In Toronto and Italy,usually people do not gather, especially after lockdown on Wuhan city, theawareness of avoiding exposure to the virus is increasing. Most people stayhome during the study period, and the family sizes in Wuhan, Toronto, andItaly on average are all around 3. Therefore, the range of the node degreesis assumed to be between 1 and 10.The Watts-Strogatz model starts with a ring of N vertices in which eachvertex is connected to its 2 m nearest neighbors ( m vertices clockwise and m counterclockwise). Each edge is connected to a clockwise neighbor withprobability p and preserved with probability 1 − p [19], where the degreedistribution is P ( k ) = min( k − m,m ) X n =0 (cid:18) mn (cid:19) (1 − p ) n p m − n ( pm ) k − m − n ( k − m − n )! e − pm . p →
1, the expression reduces to a Poisson distribution as follows P ( k ) = m k − m ( k − m )! e − m . 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 Table 6, respec-tively. We parameterized the model using the MCMC approach [24] byMATLAB R2016a according to the number of newly confirmed cases andthe cumulative number of cases reported by the Health Commission of HubeiProvince [25] and WHO [26].The rate at which the fraction of the cumulative number of cases changesis dc k /dt = ξi k , where c k ( t ) represents the fraction of the cumulative numberof infected individuals with degree k . The number of newly infected can beexpressed as P k = [ c k ( t ) − c k − ( t )] N k , where P k represents the number of new cases with degree k , and N k representsthe total number of individuals with degree k . We run the MCMC simulationfor 20000 iterations to fit the value of P k .Zhou et al. showed that the median time from illness onset (i.e., beforeadmission) to discharge was 22 days (IQR 18-25), whereas the median timeto death was 18.5 days with IQR between 15 and 22 days [27]. We assumean exponential distribution for the time to recovery for asymptomaticallyinfected individuals, symptomatically infected individuals, and hospitalizedindividuals. This results in the recovery rates γ a = γ = γ h = 1 /
22 perday, and the mortality rate, µ is 1 / . ǫ = 1 /
7. Qiuet al. reported that around 30% −
60% of people infected with COVID-19are asymptomatic or only have mild symptoms, and their transmissibility islower, but still significant [28]. Thus, we assume that the probability that aninfected individual is asymptomatic is 1 − δ = 0 .
6, and σ = 1 for simulations.We divided the Wuhan epidemic into four phases according to the re-ported data [3]. The first phase is before lockdown on Jan 23, 2020. Thesecond phase is between Jan 24, 2020 and Feb 1, 2020 when the hospitalswere short of beds. The third phase is between Feb 2, 2020 and Feb 6, 2020when the Thunder God Mountain Hospital (TGMH) and Fire God Mountain14ospital (FGMH) were put into use. The fourth phase began when door-to-door screening was implemented on Feb 7, 2020 and TGMH, FGMH, andMobile 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 cityannounced the emergence and schools and universities in Toronto were closedon Mar 18.The study period for Italy was divided into two phases. The early epi-demic phase was between Jan 31, 2020 and Mar 8, 2020 when the infectionwas spreading through the northern provinces. The second period begins onMar 9, 2020 when the national lockdown started. The parameters and initial conditions of simulations for Wuhan on theWS network are shown in Table 1. The probability of transmission throughadequate contact is estimated by MCMC. The 5000 realizations of the basicreproduction numbers derived for Wuhan using the parameter values listedin Table 1 are shown in Table 2.From Jan 11 to Mar 31, we estimate that the mean reproduction numberon the WS network decreases from 3 .
41 in the first phase to 5 . × − inthe fourth phase. The epidemic on the WS network is shown in Figure 1.Up to Jan 23, 2020 when Wuhan lockdown started, the estimated epidemicsize is 3 . × . During the second stage, after the lockdown of Wuhan andbefore the TGMH and FGMH were put into use, the predicted final size is2 . × . Thus, the lockdown of Wuhan reduced the expected final size by45 . . × . Hence, the city lockdown and the usage of TGMHand FGMH reduced the final size by 97 . .
70% due to the increase of healthcarecapacity.The variability of the numbers of confirmed new cases is consistent withthe variability of the reproduction numbers listed in Table 2. In the first twophases, the epidemic spread rapidly with larger reproduction numbers thatare larger than 1, and the numbers of infected cases increase. In the last twophases, the spread of disease is controlled, and the reproduction numbersare smaller than 1. In the third phase, because a large number of cases areconfirmed by door-to-door screening and expanded healthcare capacity, the15umulative number of confirmed cases increased. On the other hand, theepidemic will die out because the reproduction number is less than one. Inthe fourth phase, the spread of the disease has been under control with thereproduction number being less than one. Hence, the number of new casesdecreases.The parameters and initial conditions of simulations for the GTA areshown in Table 3. The 5000 realizations of the basic reproduction numbersderived for Toronto using the parameter values listed in Table 3 are shown inTable 2. The reproduction numbers are much smaller due to social distancingpolicy, school closure, as well as behavior changes. The summary of thesimulations is shown in Table 4 and Table 5. Simulation results are shown inFigure 2. The peak size is 60.19 (95%CI: 47.42-72.97), the peak time is Apr2 (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 inTable 6. The 5000 realizations of the basic reproduction numbers derivedfor 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 thatin the first phase due to the awareness of the severity of the epidemic. Thesummary of the simulation results is shown in Table 7 and Table 8. Figure3 shows that the peak number of new cases is 5492 (95%CI: 5277-5708) onMar 26 (95%CI: Mar 24-Mar 27), and the final size is 2 . × (95%CI:2 . × − . × ).
4. The impact of mitigation strategies
The close contacts identified by contact tracing will be quarantined dueto exposure to COVID-19 to see if they become sick. To evaluate the impactof mitigation strategies on the spread of COVID-19, Model (2.1) is rewrittenas follows 16 ds k dt = − βks k P k ′ [1 − (1 − k ′ − k ′ P ( k ′ | k ) i k ′ )(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ )] − qs k + λsq k , dsq k dt = qs k − λsq k , de k dt = βks k P k ′ [1 − (1 − k ′ − k ′ P ( k ′ | k ) i k ′ )(1 − k ′ − k ′ P ( k ′ | k ) σa k ′ )] − ǫe k , da k dt = (1 − δ ) ǫe k − γ a a k , di k dt = δǫe k − γi k − µi k − ξi k , dh k dt = ξi k − γ h h k − µh k , dr k dt = γi k + γ a a k + γ h h k , dd k dt = µi k + µh k , (4.1)where sq k = SQ k /N k represents the fraction of quarantined individuals withdegree k . The parameter q represents the rate at which susceptible individu-als are quarantined, and λ represents the rate at which the quarantined anduninfected close contacts transfer to the susceptible compartment again. Inthe simulations, we let λ = 1 /
14 to approximate a mean time of 14 days inthe exposed state.For Wuhan, the cumulative number of infected individuals after lockdownand 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 capacityare effective in controlling the numbers of confirmed cases.For Toronto, the number of newly infected individuals and the cumula-tive number of infected individuals produced on the WS network after imple-menting additional containment strategies besides school closure are shownin Figure 5. We simulated the scenarios of implementing various containmentstrategies for Toronto. Simulation results showed that personal protection,reducing the node degrees of symptomatically infected individuals, and quar-antine of close contacts are effective in reducing the peak epidemic size andfinal epidemic size. Reducing the transmission rate β , by x % also reduces R by x %. When β is reduced by 20% by personal protection or social distanc-ing, the peak occurs one day earlier, and the final epidemic size is reduced byaround 18%. When β is reduced by 40%, the peak occurs two days earlier,and the final epidemic size is reduced by around 33 . q = 1 / . q = 1 /
4, the peak appears five days earlier, and the finalepidemic size is reduced by 58 . . . . . . .
55% when thenode degrees of symptomatically infected individuals are reduced by 1, 2,and 3, respectively.For Italy, the number of newly infected individuals and the cumulativenumber of infected individuals simulated on the WS network after imple-menting hypothetical containment strategies are shown in Figure 6. Variousscenarios of implementing mitigation strategies showed that the peak epi-demic size and final epidemic size in Italy are greatly reduced by personalprotection, social distancing, behavior change of symptomatically infectedindividuals, and quarantine. The simulations show that the peak would havearrived earlier if the containment had been intensified.When the probability of contact transmission coefficient β , is reduced by20% by personal protection or social distancing, the peak occurs one dayearlier, and the final epidemic size is reduced by 21 . β is reducedby 41 . q = 1 /
8, the peak occurs six days earlier, and the final epidemic sizeis reduced by 52 . q = 1 /
4, the peak occurs eight days earlier,and the final epidemic size is reduced by 67 . . . . . . .
51% when the node degrees of symptomatically infected individualsare reduced by 1, 2, and 3, respectively.
5. Summary and Discussions
Modelling the dynamics of COVID-19 epidemics and assessment of miti-gation strategies could be instrumental to public health agencies for surveil-lance and healthcare planning. For the models to be reliable, the simulatedepidemic must account for the stochastic and heterogeneous contact amongindividuals. Hence, we developed a network model that captured the con-tact heterogeneity among individuals. We applied the model to analyze thetransmission potential, and mitigation strategies for curbing the spread ofCOVID-19 epidemics in the cities of Wuhan, China and Toronto, Canada,and in the Italian Republic. The epidemic threshold derived from our net-work model can be used to predict the risks of spreading scenarios. We also18rovided an explicit expression of the final epidemic size, which facilitatesestimating the scale of an outbreak for any region of interest. Our resultsprovide insights in defining a mathematical framework for the analysis andcontainment of epidemic transmission in the real world.The flexible network model framework can simulate a wide range of mit-igation strategies can be examined by the flexible model framework. It canbe extended to quantify the effectiveness of personal protection, social dis-tancing, reducing the node degree of infected individuals, and quarantine onthe dynamics of epidemics in different regions. When the mitigation strategyis intensified, the model predicts that the number of new cases peaks earlierand the final epidemic size is greatly reduced.The social contact network structure and parameter values determine thetransmission and epidemic course of such an emerging infectious disease. Wechoose the Watts-Strogatz to approximate real social networks, when theexact contact tracing data is unavailable. We assumed that the range ofthe node degree is between one and ten for each network in the absence ofreal contact tracing data, that is, on average each day an infected personwould have between one and ten contacts where they could transmit theinfection to another person. In the real world, the range of the degree willdepend on the distribution of the household sizes of the region and time beingstudied. Moreover, the network structure can be altered by behavior changeof individuals during epidemics. When this happens, the network structurecan be adapted in our model to predict the impact of these changes on theepidemic threshold, epidemic peak value, peak time, stopping time, and finalsize of infected population.The epidemics for the three places under study were fitted very well byour model with a small confidence interval. Hence, the forecasts by themodel can be reliable. We did not provide the stopping time since too manyuncertainties may affect the duration of the epidemics. As shown in thesimulations, the transmission dynamics for four phases in Wuhan are quitedifferent due to the variability on the intensity of interventions, the avail-ability of healthcare facilities, as well as the utilization of personal protectiveequipment (PPE). The dynamics in the first phase is quite different fromthat in the second phase for Toronto. The same phenomenon is observed inItaly.At the early stage, almost no interventions were implemented, and thepublic was not aware of or did not pay much attention to the severity ofthe highly contagious disease. With the increase of the number of reported19onfirmed cases and with the aid of social media, the public becomes aware ofthe severe consequence and has increased the level of personal protection andhave avoided gathering, so that the reproduction number decreases and theestimated epidemic size declines by reducing the node degree of the network.Similarly, after applying the mitigation measures in Italy on March 8 andclosing all schools in Toronto on March 18, the epidemics tend to be undercontrol.Hence, social distancing, self-isolation, quarantine, the utilization of PPE,and other measures of avoiding exposure to the virus can greatly reduce thesize of infection during the COVID-19 outbreak. Therefore, it is essentialto raise the awareness of these countermeasures to avoid contact betweenindividuals. The possibility of recurrent outbreaks of the disease cannot beoverstated. Even if the number of new cases is declining, it is still necessaryto continue taking protective measures to prevent the occurrence of futureoutbreaks. The social media should warn the public not to relax their vigi-lance against the contagion of such a highly infectious disease.
Acknowledgements
LX is funded by Fundamental Research Funds for the Central Universitiesof China. WS is funded by the National Science Foundation for Young Schol-ars of Heilongjiang Province QC2018004, and Fundamental Research Fundsfor the Central Universities of China (HEUCF181106, 3072019CF2411). JGEFis supported by multidisciplinary grant SIP-IPN 20196759. HZ is supportedby Canadian Institutes of Health Research (CIHR), Canadian COVID-19Math Modelling Task Force, and York Research Chair program of York Uni-versity.
Conflict of Interests Statement
We declare that there is no conflict of interest associated with this work.
References [1] R. Anderson, R. May, Infectious Diseases of Humans: Dynamics andControl, Oxford, Oxford University Press, 1991.[2] A. L. Ziff, R. M. Ziff, Fractal kinetics of COVID-19 pandemic, medRxiv, https://doi.org/10.1101/2020.02.16.20023820 , (2020).20 able 1:
The parameter values and initial condition of four-phase simulationsfor Wuhan.
Phase 1 is between Jan 11 and Jan 23, Phase 2 is between Jan 23 and Feb1, Phase 3 is between Feb 1 and Feb 12, and Phase 4 is between Feb 12 and Mar 31.
Parameter Meanvalue Std 95% CI References β (Phase 1) 0 . . × − [0 . . β (Phase 2) 0 . . × − [6 . × − , 0 . β (Phase 3) 2 . × − . × − [0, 7 . × − ] MCMC β (Phase 4) 7 . × − . × − [0, 2 . × − ] MCMC ξ (Phase 1) 0 . . . . ξ (Phase 2) 0 . . . ξ (Phase 3) 0 . . . . ξ (Phase 4) 0 . . . . P k N k (0) 11081000 − − [29] P k S k (0) 11080770 − − Calculated P k E k (0) 200 .
26 18 .
88 [163 .
25, 237 .
26] MCMC P k A k (0) 17 .
81 4 .
13 [9 .
71, 25 .
90] MCMC P k I k (0) 11 .
78 6 .
06 [0, 23 .
66] MCMC P k H k (0) 41 − − [25] P k R k (0) 0 − − Estimated P k D k (0) 0 − − Estimated
Table 2:
Basic reproduction numbers computed by MCMC on the WS network.
Location Period Meanvalue Standardderivation 95% CIWuhan 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 . × − . × − [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]21 able 3: The parameter values and initial condition of simulations for Toronto.
Phase 1 is from Jan 26 to Mar 18, and Phase 2 is from Mar 18 to Mar 29.
Parameter Meanvalue Std 95% CI References β (Phase 1) 7 . × − . × − [5 . × − , .
01] MCMC β (Phase 2) 1 . × − . × − [0 , . × − ] MCMC ξ (Phase 1) 0.1421 0.0734 [0, 0.2858] MCMC ξ (Phase 2) 0.1140 0.0709 [0, 0.2530] MCMC P k N k (0) 5928000 - - [30] P k S k (0) 5927990 - - Calculated P k E k (0) 4.12 4.83 [0, 13.58] MCMC P k A k (0) 2.43 2.22 [0, 6.78] MCMC P k I k (0) 3.22 2.89 [0, 8.88] MCMC P k H k (0) 1 - - [26] P k R k (0) 0 - - Estimated P k D k (0) 0 - - Estimated Table 4:
The peak number of new cases, peak time, and final epidemic size aftercontainment strategies are implemented in Toronto.
Scenarios Peak size (95%CI) Peak time (95%CI) Final size (95%CI) β . β . β q = 0 60.19 (47.42, 72.97) Apr 2 (Mar 29, Apr 7) 2712 (1603, 3820) q = 1 / q = 1 / Table 5:
The peak number of new cases, peak time, and final epidemic size forToronto when varying node degrees of symptomatically infected individuals.
Degree Peak size (95%CI) Peak time (95%CI) Final size (95%CI) k k − k − k − able 6: The parameter values and initial condition of simulations for Italy.
Phase 1 is between Jan 31 and Mar 8, and Phase 2 is between Mar 8 and Mar 26.
Parameter Mean value Standardderivation 95% CI References β (Phase 1) 0.0179 1 . × − [0.0156, 0.0203]] MCMC β (Phase 2) 4 . × − . × − [0 , . × − ] MCMC ξ (Phase 1) 0.0996 0.0458 [9 . × − , . ξ (Phase 2) 0.1312 0.0250 [0.0823, 0.1801] MCMC P k N k (0) 59430000 - - [31] P k S k (0) 59429892 - - Calculated P k E k (0) 69.01 56.92 [0, 180.58] MCMC P k A k (0) 22.60 18.24 [0, 58.35] MCMC P k I k (0) 32.11 24.23 [0, 79.59] MCMC P k H k (0) 2 - - [26] P k R k (0) 0 - - Estimated P k D k (0) 0 - - Estimated Table 7:
The peak number of new cases, peak time, and final size after contain-ment strategies are implemented in Italy.
Scenarios Peak size(95%CI) Peak time(95%CI) Final size (95%CI) β . × (2 . × , . × )0 . β . × (1 . × , . × )0 . β . × (1 . × , . × ) q = 0 5492 (5277, 5708) Mar 26 (24, 27) 2 . × (2 . × , . × ) q = 1 / . × (1 . × , . × ) q = 1 / . × (7 . × , . × ) Table 8:
The peak size, peak time and final size for Italy when varying nodedegrees of symptomatically infected individuals.
Degree Peak size(95%CI) Peak time(95%CI) Final size (95%CI) k . × (2 . × , . × ) k − . × (1 . × , . × ) k − . × (1 . × , . × ) k − . × (1 . × , . × )23 igure 1: Fitting the number of reported new cases and the cumulative number of reportedcases 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 thecumulative number of reported cases on the Watts-Strogatz network. T o t a l c o n f i r m e d n e w c a s e s A T o t a l c o n f i r m e d c a s e s × B MCHTGMH and FGMHcity lockdown MCHcity lockdownTGMH and FGMH igure 2: Fitting the number of reported new cases and the cumulative number of reportedcases for Toronto on the Watts-Strogatz network. (A) Fitting the number of reported newcases on the Watts-Strogatz network. (B) Fitting the cumulative number of reported caseson the Watts-Strogatz network. T o t a l c o n f i r m e d n e w c a s e s A T o t a l c o n f i r m e d c a s e s B containmentcontainment igure 3: Fitting the number of reported new cases and the cumulative number of reportedcases for Italy on the Watts-Strogatz network. (A) Fitting the number of reported newcases on the Watts-Strogatz network. (B) Fitting the cumulative number of reported caseson the Watts-Strogatz network. T o t a l c o n f i r m e d n e w c a s e s A T o t a l c o n f i r m e d c a s e s × B containmentcontainment igure 4: The impact of the variability on the healthcare capacity on the spread of theepidemic in Wuhan on the Watts-Strogatz network.
20 40 60 80Time (days)024 T o t a l c o n f i r m e d c a s e s × A Phase 1 without any measures 20 40 60 80 100 120Time (days)024 T o t a l c o n f i r m e d c a s e s × B Phase 2 city lockdown20 40 60 80 100 120Time (days)051015 T o t a l c o n f i r m e d c a s e s × C Phase 3 city lockdownTGMH and FGMH 20 40 60 80Time (days)0246 T o t a l c o n f i r m e d c a s e s × D Phase 4 city lockdownTGMH and FGMHMCH [3] H. Wang, Z. Wang, Y. Dong, et al., Phase-adjusted estimation of thenumber of coronavirus disease 2019 cases in Wuhan, China, Cell discov-ery 6 (1) (2020) 1–8.[4] A. J. Kucharski, T. W. Russell, C. Diamond, et al., Early dynamicsof transmission and control of COVID-19: a mathematical modellingstudy, The lancet infectious diseases 20 (5) (2020) 553–558.[5] M. U. Kraemer, C. H. Yang, B. Gutierrez, et al., The effect of humanmobility and control measures on the COVID-19 epidemic in China,Science 368 (6490) (2020) 493–497.[6] M. Chinazzi, J. T. Davis, M. Ajelli, et al., The effect of travel restric-tions on the spread of the 2019 novel coronavirus (COVID-19) outbreak,Science 368 (6489) (2020) 395–400.[7] N. Ferguson, D. Laydon, G. Nedjati Gilani, et al., Report 9: Impact ofnon-pharmaceutical interventions (NPIS) to reduce COVID-19 mortal-ity and healthcare demand,
DOI:https://doi.org/10.25561/77482 ,(2020). 27 igure 5: The impact of mitigation strategies on the spread of COVID-19 epidemic inToronto on the Watts-Strogatz network. In this figure and the following figure, the dashedlines represent 95% confidence intervals. In (A) and (B), the red, purple, and green linesrepresent that the transmission rates are unchanged, reduced by 20%, and reduced by40%, respectively. In (C) and (D), the red, purple, and green lines represent the rate ofquarantine, q = 0, 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 individualsare reduced by 0, 1, 2, and 3, respectively. (A) The number of newly infected individualsafter reducing the transmission rates by personal protection and social distancing. (B) Thecumulative number of infected individuals after reducing the transmission rates by personalprotection and social distancing. (C) The number of newly infected individuals after closecontacts are quarantined. (D) The cumulative number of infected individuals after closecontacts are quarantined. (E) The number of newly infected individuals after the nodedegrees of symptomatically infected individuals are reduced. (F) The cumulative numberof infected individuals after the node degrees of symptomatically infected individuals arereduced. T o t a l c o n f i r m e d n e w c a s e s A containment 1/26 3/6 4/15 5/25 7/4Time (days)020004000 T o t a l c o n f i r m e d c a s e s B containment1/26 3/6 4/15 5/25 7/4Time (days)020406080 T o t a l c o n f i r m e d n e w c a s e s C containment 1/26 3/6 4/15 5/25 7/4Time (days)020004000 T o t a l c o n f i r m e d c a s e s D containment1/26 3/6 4/15 5/25 7/4Time (days)020406080 T o t a l c o n f i r m e d n e w c a s e s E containment 1/26 3/6 4/15 5/25 7/4Time (days)020004000 T o t a l c o n f i r m e d c a s e s F containment igure 6: The impact of mitigation strategies on the spread of COVID-19 epidemic inItaly on the Watts-Strogatz network. (A) The number of newly infected individuals afterreducing the transmission rates by personal protection and social distancing. (B) Thecumulative number of infected individuals after reducing the transmission rates by personalprotection and social distancing. (C) The number of newly infected individuals after closecontacts are quarantined. (D) The cumulative number of infected individuals after closecontacts are quarantined. (E) The number of newly infected individuals after the nodedegrees of symptomatically infected individuals are reduced. (F) The cumulative numberof infected individuals after the node degrees of symptomatically infected individuals arereduced. T o t a l c o n f i r m e d n e w c a s e s A containment 1/31 3/11 4/20 5/30 7/9Time (days)0123 T o t a l c o n f i r m e d c a s e s × B containment1/31 3/11 4/20 5/30 7/9Time (days)0200040006000 T o t a l c o n f i r m e d n e w c a s e s C containment 1/31 3/11 4/20 5/30 7/9Time (days)0123 T o t a l c o n f i r m e d c a s e s × D containment1/31 3/11 4/20 5/30 7/9Time (days)0200040006000 T o t a l c o n f i r m e d n e w c a s e s E containment 1/31 3/11 4/20 5/30 7/9Time (days)0123 T o t a l c o n f i r m e d c a s e s × F containment R in models for infectiousdiseases in heterogeneous populations, Journal of mathematical biology28 (4) (1990) 365–382.[21] P. Van den Driessche, J. Watmough, Further notes on the basic repro-duction number, in: Mathematical epidemiology, Springer, 2008, pp.159–178.[22] Y. Wang, J. Cao, A. Alofi, A.M. Abdullah, A. Elaiw, Revisiting node-based sir models in complex networks with degree correlations, PhysicaA: Statistical Mechanics and its Applications 437 (2015) 75–88.[23] Y. Wang, J. Cao, G. Huang, Further dynamic analysis for a networksexually transmitted disease model with birth and death, Applied Math-ematics and Computation 363 (2019) 124635.[24] H. Haario, M. Laine, A. Mira, E. Saksman, DRAM: efficient adaptiveMCMC, Statistics and computing 16 (4) (2006) 339–354.[25] Health Commission of Hubei Province, http://wjw.hubei.gov.cn/bmdt/ztzl/fkxxgzbdgrfyyq/ , accessed Apr 7, 2020 (2020).[26] World Health Organization, , accessedApr 7, 2020 (2020).[27] F. Zhou, T. Yu, R. Du, et al., Clinical course and risk factors for mortal-ity of adult inpatients with COVID-19 in Wuhan, China: a retrospectivecohort study, The Lancet 395 (10229) (2020) 1054–1062.[28] J. Qiu, Covert coronavirus infections could be seeding new outbreaks, , Nature(2020).[29] Hubei Provincial Bureau of Statistics, Available from: http://tjj.hubei.gov.cn/tjsj/sjkscx/tjnj/qstjnj/ , accessed Feb 17, 2020(2020). 3130] Toronto Population 2020, https://worldpopulationreview.com/world-cities/toronto-population/ , accessed Mar 17, 2020 (2020).[31] World Health Organization,