key: cord-0592930-2nu6t2ek
authors: Li, Tianyi
title: Simulating the Spread of Epidemics in China on the Multi-layer Transportation Network: Beyond the Coronavirus in Wuhan
date: 2020-02-27
journal: nan
DOI: nan
sha: 962582ae057d9ef92d2af93303e25ffcd9e7826f
doc_id: 592930
cord_uid: 2nu6t2ek

Based on the SEIR model and the modeling of urban transportation networks, a general-purpose simulator for the spread of epidemics in Chinese cities is built. The Chinese public transportation system between over 340 prefectural-level cities is modeled as a multi-layer bi-partite network, with layers representing different means of transportation (airlines, railways, sail routes and buses), and nodes divided into two categories (central cities, peripheral cities). At each city, an open-system SEIR model tracks the local spread of the disease, with population in- and out-flow exchanging with the overlying transportation network. The model accounts for (1) different transmissivities of the epidemic on different transportation media, (2) the transit of inbound flow at cities, (3) cross-infection on public transportation vehicles due to path overlap, and the realistic considerations that (4) the infected population are not entering public transportation and (5) the recovered population are not subject to repeated infections. The model could be used to simulate the city-level spread in China (and potentially other countries) of an arbitrary epidemic, characterized by its basic reproduction number, incubation period, infection period and zoonotic force, originated from any Chinese prefectural-level city(s), during the period before effective government interventions are implemented. Flowmaps are input into the system to trigger inter-city dynamics, assuming different flow strength, determined from empirical observation, within/between the bi-partite divisions of nodes. The model is used to simulate the 2019 Coronavirus epidemic in Wuhan; it shows that the framework is robust and reliable, and simulated results match public city-level datasets to an extraordinary extent.

The 2019-nCov epidemic, originated from Wuhan, China [14, 34, 39, 40] has incurred heavy casualties and tremendous economic loss. With the first case confirmed as early as Dec. 8th, 2019 [20] and effective measures implemented in a national scale in China around 50 days later, the new Coronavirus has claimed over 70,000 cases and more than 2700 deaths in China (as on Feb. 26th). The epidemic took place right before the Chinese New Year, and the massive population flow across the entire country aggravated its vicious spread. Not long after its burst in the Hubei Province in early January, the disease started to propagate worldwide, invading almost all countries in East Asia and soon reached Europe, America and Australia [16, 27] , due to the fact that Wuhan is a top-10 Chinese metropolitan and technically the single transportation hub in the central region of China, with a population as large as 11 millions. Facing the severe threat of the epidemic, since late January, the Chinese government has taken strong measures to quench the disease, mobilizing and coordinating available forces in a unprecedented scale, and thanks to the collective efforts of Chinese citizens, the epidemic has been largely put under control in most Chinese provinces by mid-February. Nevertheless, although the situation in China is picking up, recent reports indicate that the virus has begun to trigger vicious dynamics in other countries, and the epidemic will inevitably leave profound effects on the global economy. * tianyil@mit.edu Right after the burst of 2019-nCov, studies set out to model the epidemic and use simulation results to nowcast and forecast its intensity [e.g., 21, 22, 32] . The 2019 Coronavirus is compared with the SARS in 2003, and people discovered that this time the virus is more contagious but less fatal. Arguably, the incubation period of 2019-nCov is believed to be longer than SARS, and the early symptoms are less salient [17] ; to a large extent, it is for these two reasons that this epidemic was not paid enough attention to in its early period, and an exponential growth after the sufficient breeding soon took off during the Chinese New Year. Wu et al. [35] is the first formally published simulation model for 2019-nCov. Because this disease is originated only in Wuhan and is carried to other cities, Wu et al. [35] assembles an SEIR model for Wuhan with population in-and outflows, treating infections outside Wuhan as imported cases. The model is calibrated with the time series data of the confirmed cases in a list of global cities, and a posterior estimate of the basic reproduction number R 0 = 2.68 is carried out. Simulation models similar to the flavor of Wu et al. [35] are of great importance for policy decisionmaking, whose findings may significantly help guide the rescue and response for the epidemic. Nevertheless, a few problems arise with the model in Wu et al. [35] , hampering it to reach sufficient resolution. One big concern is that the model does not differentiate means of transportation during the population flow. The public transportation system is composed of airlines, railways, sail routes and buses (highways). It is very important to note that especially in China where public transportation is exploited to a great extent, different means of transportation leads to different contact rates, and thus different transmissivities of the disease during travel. The spread of epidemic is substantially easier on trains or buses than on airplanes, while for (private) car travels, cross-infection is essentially negligible since the transportation is end-to-end. The differentiation of means of transportation is a must for high-resolution epidemic models regarding inter-city population flow since people spend non-trivial time on the route and may engage in various activities. Moreover, there exists crossinfection during travel due to path overlaps, which could be accounted for if means of transportation are not separated. Since railway, sails and bus travels are often not end-to-end, the spread of disease will likely take place along the way, among people taking the same vehicle yet having different destinations. Again, such an infection scenario is not negligible in China where the transportation system is crowded. Note that similar cross-infection concerns may arise in transportation service locations such as train stations and airports; yet concerns in these places are not as severe as during the travel where people are more densely located and have little choice of isolating from their neighbors who have unknown points of departure and destinations. Second, the model initiates an SEIR compartment model only in Wuhan but not in other cities. This only holds true if the disease is not sufficiently contagious, or after the quarantine procedures are successfully implemented in all other cities, and this model could not capture the local evolution dynamics of the imposed diseases, which is non-trivial at least for the other cities in Hubei Province. Ideally, the same SEIR dynamics could be initiated in every city, contributing to the spread of the disease in a nationwide landscape. Third, the model neglects an important aspect of the population flow in China, the non-trivial occurrence of transit events in public transportation. Although highly developed, the Chinese transportation system could not (even get close to) realize cheap end-to-end travels between the over 340 prefectural-level cities, as the primarily concerned level of resolution in the epidemic spread. Yet Chinese citizens are also much more likely than citizens of other countries to travel a long distance to a major city (e.g., to make a living) due to an imbalanced development across the nation, which gives rise to the non-trivial role of transfer during inter-city public transportation. Without taking care of the transit issue, the direct-import model may possess a fundamental system error. Last, the model in Wu et al. [35] allows population out-flow of Wuhan from all four SEIR compartments. This may contain another system error as it is more appropriate to assume that the infected population are not traveling but instead stay local. Also, since the model did not trace the source of the population inflow, it has to put the recovered population into the susceptible compartment and essentially assumes repeated infections, which are often not the case for virus-triggered diseases. The above problems could be resolved by assembling a network of the transportation system on top of the local evolution of epidemics and formulate open-system compartments, which will bring the model resolution to the next stage. On this network, nodes represent population districts (communities/cities/countries) and edges describe transportation availability between nodes. An identical compartment model (e.g., SEIR model) is initiated at each node, which generates its own dynamics of the epidemic under the in-and out-flow of population that it exchanges from the overlying transportation network. Among various methodologies in epidemic research [28] , this modeling approach has been widely taken by previous studies, and simulators are built to study the spread of epidemics as well as help design corresponding control policies after the burst of SARS [18] , H5N1 [13] , H1N1 [4] , on a national [2, 10, 11] or global [18] scale. In these studies, the transportation system under concern is often considered as aggregated or single-layer [e.g., the airline network, 7] with few exceptions [e.g., 2], and people make simplifying assumptions to determine the flow matrix, which is often considered as stochastic [e.g., 18 ]. Based on existing works, one notes that in a finer resolution, the public transportation system could be further modeled as a multi-layer network [3, 8, 9] , with each layer representing a specific means of transportation [1, 5, 6, 12, 19, 33, 41] . On this network, nodes are maintained at different layers, in which case the network is sometimes termed as "multiplex" networks [25, 26, 30] . As in the single-layer representation, transportation takes place along network edges [23] , captured by a set of flowmaps that record the flow between each pair of connected nodes on each layer. One expects that the multi-layer representation of the transportation system outstands the aggregated single-layer representation since by differentiating means of transportation, the model could account for different diffusion properties in the spread of epidemics; from the perspective of policy analysis, such an increase of model resolution might provide valuable insights. In this study, upon a multi-layer network model for the Chinese inter-city public transportation system, a simulator for the spread of epidemics in Chinese cities is built. On each of the over 340 prefectural-level cities, an identical SEIR model similar to the type in [35] is assembled to characterize the local dynamics of the disease. Per the above discussion, the flow model on the network accounts for a set of important realistic concerns, including the transit of inbound flow, the cross-infection on public transportation media due to path overlap, the deactivation of outflow of the infected population, and the unlikelihood of repeated infections. Inspired by the realworld situation of Chinese administrative districts, the model also adopts a bi-partite structure that partitions nodes into central cities and peripheral cities; this division determines the flow strength between connected nodes. Essentially, this multi-layer network model could serve as a general-purpose simulator for the evolution of epidemics in China (and potential other countries as well), which extends beyond the spread of 2019-nCov in Wuhan. Given the epidemiological parameters (basic reproduction number, incubation period, infection period, zoonotic force), geological sources (single location or multiple locations) and occurrence time (which determines the windowing of the seasonal flowmap) of a specific epidemic, the model could simulate the spread of the disease on all prefectural-level cities in China. With edge connectivities and the flowmaps being easily updated according to potential changes in public transportation infrastructures and in the distribution of population, this modeling framework could be of significant use for future policy decision-making on the control and emergency response for epidemics in China.

Consider the transportation network G. At each node i ∈ G, an SEIR model [e.g., 24, 31] is assembled, where the population P i is divided into 4 compartments (Susceptible, Exposed, Infected, Recovered):

Adopting the notations in [35] , in the base model, the dynamics are governed by:

where R 0 , D E , D I are the basic reproduction number, the incubation period and the infection period, respectively; z = z(t) is the zoonotic force, which is only nonzero for a certain period at the source node of the epidemics. Note that the venue from E i directly leading to R i is cut off, which exists in the general SEIR model. These dynamics work for a closed population which could not account for the spread of diseases due to the population flow between cities (nodes). The flow is maintained on the overlying transportation network G which is cast as multi-layer, with each layer representing a specific means of transportation, including airline (A), railway (R), sail route (S) and bus (B). We denote the set of layers as Q = {A, R, S, B}. These layers share the same set of nodes V = {i}, the 347 prefectural-level cities in China, while having different sets of edges L q for q ∈ Q according to the availability of transportation infrastructures. Edges are undirected since the transportation between two cities is bi-way. Note that we do not have a separate layer for the (private) car travel and instead assume that cars share the routes in the bus layer (see following sessions for further discussion). International population flow are not included in the current model, which are mainly carried out by airlines. Therefore we

Two properties differentiate means of transportation, besides layers' specific edge connectivities determined by the availability of infrastructures. First, for each q, there is a specific transfer rate T R q i ∈ [0, 1] for every city i, which represents the proportion of inflow to city i that will be in transit and leave for the next destination, thus not entering the population stock P i . Such transfer rates will be small for small cities and large for cities that are regional transportation centers. Ideally, the transfer rate will also be higher for certain transportation means such as railways and lower for others, say airlines (direct flights are primary among domestic air travels); yet for simplicity in the current model we assume that it is the same for all means of transportation at a certain city, i.e., T R q i = T R i . Second, on different transportation media the likelihood that an exposed person will infect others during travel is not the same, which is represented by an in-travel basic production number R q T > 1. Expectedly, on airplanes, such a number will be lower than that on trains or buses where people talk more often, and in cars, the number is essentially 0; it is safe to assume that R q T is homogeneous on different edges of layer q, i.e., the value only depends on the means of transportation. Features of the multi-layer network are illustrated in Figure 1 . The in-travel transmissivity R T characterizes the crossinfection of the disease during public transportation. Importantly, cross-infection has spillovers due to path over-laps on public transport media, which deals with the topologies of network layers: a patient traveling from city j to city i by railway will likely spread the virus to everyone in the train during the travel between j and i, including those who might get on the train at a different city k. This effect is pervasive in the real world, especially in China where the public transportation is crowded, and could be exempted only for an end-to-end transportation (e.g., airlines, cars). To model this path-dependent feature, we make two assumptions here ( Figure 2 ). First, it is assumes that only those people that get off the train at the same destination i as the patient will be infected. Implicitly, this is assuming that passengers sharing the same destination are more likely to stay close during the travel and the getting-off, which does not deviate from realistic situations. In a complete manner, the spillover will likely take place between any two routes that share a finite part; yet this will open a new dimension in the computation of the spillover and incur redundant complexity. Second, we assume that the strength of the spillover effect is proportional to the ratio of the (shortest path) distance between city k and i and the distance between j and i, for a node k that lies between j and i in the path; for a node k that lies beyond j in the path (i.e., earlier in the path), the k → i distance is thresholded by the j → i distance.

Base on the established model setting, we formulate the flow equations. On the transportation network, we track the in-and out-flow of the exposed (E), susceptible (S) and recovered (R) population of each city. For each means of transportation q, the flowmap between nodes is characterized by a matrix F q = {f q j,i }; the matrix is not guaranteed to be symmetric since the flow f q i,j from city i to j is not equal to the flow f q j,i from city j to i. Using values on the |Q| = 4 flowmaps, at time t, the inflow of the exposed population to city i summarize over all cities on the transportation network and over all means of transportation:

with the fraction of people in transit deducted. Here

is the adjusted exposed population flow from city j to i by means q, taking care of the spillover effect in crossinfection, in which d q i,j represents the shortest path distance between i and j on layer q. And we have

and

which are the time-stamped proportion of the exposed and recovered population among the total outflow population of city i, respectively. Therefore, (1 − µ i − η i ) is the proportion of the susceptible population among the total outflow of city i. The exposed population among the total outbound of i is contributed by both the exposed originated from city i (∆E out i , see below) and the exposed among the transit of inbound (second term in the numerator of µ i , proportion to ∆E in i by a term of the transfer rate), and is dynamically updated at each step; we denote the proportion by µ i (t). Similarly, we track the stock (via proportion η i (t)) of the recovered population among the total outflow of a city, which is also updated at each step. This tracking is an important caveat in this model since we consider that a recovered person will unlikely to be infected again. Therefore, by this means, only the susceptible population (1 − µ i − η i ) among the total outflow along each edge will be subject to cross-infection during travel.

The outflow population from city i's population P i is the total outbound flow minus the transferred inbound flow, which is contributed to by the S, E, R compartments (not I), different from the case in Wu et al. [35] , where it assumes that all four compartments contribute to the outbound flow of Wuhan. Proportionally, the outflow of the exposed is:

Note that the flowmaps and the transfer rates should ensure that at each city, the inbound transfer flow should

Illustration of the spillover of cross-infection during travel due to path overlap. If there are exposed cases in the population flow from city j to city i, then cross-infection takes place on all routes k −→ i that shares a part with the route j −→ i. The spillover f q j,i µj(R q T − 1) of exposed cases is allocated proportionally to the susceptible population f k,i (1 − µ k − η k ) among the outflow of city k, weighted by the shared length of the two route. Routes that share a finite part with j −→ i but not terminate at i are not considered in the spillover. always be smaller than the outbound flow, i.e.,:

For the recovered population, the inflow is tracking all the recovered people (through η j ) upon arrival:

and the outflow is proportional to ∆E out i :

(10) For the susceptible population, according to the flow balance, ∆S in i is the total un-transferred inflow subtracting the recovered and the exposed inflow:

and

In the end, we arrive at the open-system SEIR model accounting for the transportation flow:

with µ i and η i updated according to (5) and (6) . Note that in the model we don't assume a balance of the inbound and outbound flow at each city; there is no such equilibrium for the population flow in real practice, and each city's transient population is in constant dynamics. However, the inequality (8) regarding the inbound transit and the outbound should always hold at each city.

The model incorporates 4 layers of the public transportation system (airline, train, sail, bus), which corresponds to the situation in China but may be applicable to other countries as well. Based on empirical observations, simplifying assumptions on the topologies of the transportation network in China are to be made. Most importantly, the multi-layer network could be further represented by a bi-partite structure (i.e., factor graph [38] ), with nodes V divided into two categories: central cities V c and peripheral cities V p (Figure 3 ). Upon this division, we assume different flow strength on edges between V c ↔ V p , V c ↔ V c and V p ↔ V p when constructing the flowmaps F . We are also able to assign a homogeneous value of the transfer rate T R to each category of nodes. These treatments bring down the model's parameter space to a feasible region (2 for T R, 4×3 = 12 for F ). This division of cities' roles is valid in China, since Chinese cities are often categorized on a level base [ Note 1 ]. In the current study we apply a two-bin partition and regard all provincial capital cities and level 0-2 cities as central cities V c (54 counts, Table 1 ), with the rest 293 cities belong to V p ; expectedly, in a finer scale, cities could be categorized into more bins, and the flowmap will then have a multi-partite structure similar to the celebrated stochastic block model [15] . When real datasets could be acquired, these parameters may be determined in a more objective manner. Ideally, the flow f q j,i from node j to node i by means q, is supposed to be cast as a time series that incorporates seasonal features but approximately remains invariant on the annual basis. Then for a specific starting time of an epidemic, a cursor is placed on the annual curve and a window of the flowmap is then cut out from this point on for the usage in simulation. Moreover, in cases where only the aggregated flow of all transportation layers between two cities F i,j = q {f q j,i } is available, one may apply a multinomial logit model to determine the flow for each layer, as in [37] . For the transfer rate, if the city's (average/annual) transient population P T i is available, then T R could be calculated by comparing the aggregated outflow and the transient population:

However, this detailed formulation might be subject to error when data quality is not guaranteed, since one has to ensure that the outbound is always greater than the transient population.

As in many countries, different layers of the public transportation system in China play different roles in serving the demand of population flow, and therefore edge connectivities are different for each layer [Note 2] . For airlines, we obtain public data on the airline schedules for 90 cities; for railways, there are 11 primary railroads in China which occupy the major railway passenger flow, connecting ∼ 120 cities into the network. For sails, there are roughly 19 passenger lines along the east coast of China and the Yangtze river; yet sails take a very small part (< 2%) in Chinese public transportation. For buses, schedules are difficult to collect, but according to empirical observations, we construct edges on the bus layer between cities within a mid-range geological distance (150km), and between a provincial capital city to other cities in the province. The four layers (G A , G R , G S , G B ) of the Chinese inter-city public transportation network are shown in Figure 4 . Based on the determined edge connectivities, we calculate nodes' centrality measures on the network, which could be used to distinguish central cities from peripheral cities. Three measures are computed for each node i ∈ V : node activity A i [26] , which is the count of layers in which node i has at least one edge; betweenness B i :

and coupling strength C i [23] , defined as:

where σ s,t in both cases indicates the number of shortest paths between nodes s and t, σ i s,t are the number of such shortest paths that traverse node i, and σ coupled s,t is the number of shortest paths that use edges on more than one layer. These three measures characterize the centrality of nodes on a multilayer network from different perspectives: A i indicates the presence of i on different layers, B i indicates its significance in bridging network paths, and C i indicates the extent of its embedding into the multiplexity of the network. Note that B i is calculated on the aggregated network that is the superposition of multiple layers. Results of these centrality measures for the 54 central cities are shown in Table 1 , with their average values compared to those of peripheral cities ( Figure  5 ). Many central cities are active on 3 or 4 layers, with an average A i of 2.87, compared with 1.65 for peripheral cities. The average betweenness B i of central cities are 35 times larger than that of peripheral cities, and ranking #1 − 35 largest betweenness are all on central cities. For coupling strength C i , central cities and peripheral cities are almost indistinguishable. This is consistent with expectation, since on one hand, central cities have multiple means of transportation that likely increases its coupling strength, while on the other hand, central cities often have well-developed airlines that are always the shortest paths, which reduces C i ; these two competing forces decide that the coupling strength of central cities is not going to be clearly larger than that of peripheral cities. Overall, these results suggest that our hard division of V c and V p according to public information is reasonable on the basis of real edge connectivities.

We use the model to simulate the spread of the Wuhan Coronavirus in Chinese cities under assumed parameters determined from empirical considerations. For epidemiological parameters, we let R 0 = 2.68, D E = 6 days, D I = 2.4 days, according to Wu et al. [35] ; we starts the simulation at Dec. 8th, 2019, when the first case was reported, with a nonzero zoonotic force z = 27 cases/day (the total number of confirmed cases in December) before Jan. 1st the date of the closure of the Wuhan Huanan seafood market, which is the claimed source of the virus. We tested the values suggested by Wu et al. [35] (z = 86 cases/day extending 30 days) and results show a misfit (see below) three times larger than using our values for z; therefore the values in Wu et al. [35] are not used. For parameters of the transportation network, the in-travel basic production number R T is assumed to be:

, smallest for airlines, largest for buses. Flowmaps are determined based on empirical constraints (all units are people/day). Airlines are the major transportation means between central cities, and we determine f A cc/cp/pp = 1000/500/0, assuming no airline transport between peripheral cities, four flights between central cities per day, and half the amount between a central city and a peripheral city. For railways the flow is determined asf R cc/cp/pp = 2000/200/500, and f i,j =f i,j /d i,j , i.e., passenger flows are inversely proportional to the shortest path distance between cities, under the carriage capacity 2000 people/train. For sail routes, f S cc/cp/pp = 100/100/100, which occupy a small fraction of the total flow. For buses, we let f B cc/cp/pp = 0/3000/1000, assuming no bus travel between central cities; f A pp = 0 and f B cc = 0 reinforce the bi-partite structure of the model. Buses are heavily used in provincial transport, for travelers in peripheral cities to go to either the local central city or nearby peripheral cities. We let the numbers (3000 and 1000) be three/two times the normal flow (50 people/bus, 20 or 10 buses between central-peripheral and peripheral-peripheral) to account for the private transportation by cars that nevertheless contributes to the entire population flow, which is mostly provincial rather than regional (to a nearby province). Implicitly, as car travels are aggregated into bus travels, the R B T represents the effective value watered down from a real transmissivity on buses. Note that the assumed constant flowmaps do not consider seasonal effects, which might be significant in certain cases, exactly like for the Wuhan epidemic which took place right before the Chi-nese New Year when massive transports are carried out. Although not obtained from a full-parameter inversion, the values for R T and F that we use are nevertheless well tested by extensive forward simulation runs and are confirmed as quasi-locally-optimal values. The current runtime of the simulator is approximately 55 seconds per simulation time step, estimated on a 2018 Macbook Pro; for a simulation range of 50 time steps, one forward run takes rough 0.75 hours. Given the large parameter space even for the simplified model (2 for T R, 12 for F , 4 for R T and 5 epidemiological parameters), a fullparameter inversion is not feasible on personal computing devices. Per the above discussion, with other parameters determined through empirical considerations, T R c and T R p are set as open parameters and we initiate a partial inversion. At each run, the model is simulated for 53 time steps (days), from the date of the first case (Dec. 8) minus an incubation period (6.4 days), to the day before the Chinese New Year's Eve (Jan. 23). This stopping time is reasonable since most population flow in China during winter travels took place before this date, right after which the government took urgent measures and asked the entire domestic population to self-quarantine and abandon inter-city public transportation. After Jan. 23, confirmed cases of the epidemic are collected from each prefectural-level city and are reported to the public; we use this dataset for inversion [Note 3] . We observe that the fraction of confirmed cases in each city ρ among national headcounts remains roughly invariant after around Feb.15th ( Figure 6 ) due to the reduction of daily new headcounts as the effect of nationwide measures, with Wuhan having ∼ 60% of all confirmed cases and all cities in the Hubei Province having > 80% cases. Therefore, instead of matching the simulated time series with the entire data series, for the inversion we minimize the misfit between the simulated and realistic ρ obtained at the end of the two time series, summing over each city except Wuhan in prevention of double-counting because i ρ i = 1. This choice of misfit is also consistent with the fact that we don't have a reliable estimate of the zoonotic force, and also that the transportation network setting focuses more on the relative strength of the epidemic in each city. Since the simulation stops at Jan.23 but the confirmed cases are reported gradually afterwards, when calculating the ρ of the simulation result we regard the exposed (E) as additional infected (I) cases as they will eventually enter the stock, a valid treatment given that the incubation period is unlikely to be as long as 30 days (from Jan. 23 to Feb. 23). Although real datasets for flowmaps are not applicable and a full-parameter inversion is not conducted, quite surprisingly, results of the coarse partial inversion nevertheless demonstrated great fitting performance of our model. The best-fit T R c = 14.95%, T R p = 0.36%, i.e., on average around 15% inbound flow at central cities are in transit, and the number is near 0 for peripheral cities, which is largely consistent with the empirical observation. We show the best-fit ρ of Chinese provincial dis-Jan. 24 Feb . Figure 7) . It suggests that our results (red dots) recover the data (bars) to a satisfying extent: the severe situation of the disease in Hehan, Hunan, Shandong and Jiangxi is correctly revealed, while the spread in provinces such as Chongqing, Anhui, Zhengjiang, Guangdong and Heilongjiang is underestimated. This result might be consistent with the public news that after the burst of the epidemic there are uncommonly large population flow directed into these provinces from Wuhan, a situation not captured by our simplified flowmaps. Entire best-fit time series are also shown (right, Figure 7) , on five example cities: Wuhan (origin of the epidemic), Beijing (capital of China), Huanggang (peripheral city in Hubei), Harbin (central city outside Wubei) and Kiamusze (peripheral city outside Wubei). The data (solid lines) are S-shape as the epidemic is gradually under control after Chinese New Year, while the simulation time series (dash lines) are demonstrating exponential growth, as one expects from the SEIR model. Nevertheless, applying a uniform time-shift (∆t = 20) assuming implicitly that all cities took measures at the same time, which was roughly the case in China, simulation results match the initial part of the real S-curve by a large margin. However, the matching of the two curves should not be over-interpreted as the model is always able to generate exponential growth and the value of the time-shift is not warranted. The absolute values of the confirmed cases, instead of the fraction ρ, are also fit reasonably well (bars in Figure 7 right), suggesting that our choice of model parameters makes certain sense. Overall, given the coarse treatment during the parameter determination and data-processing process, these partial inversion results are believed to be acceptable, suggesting that the model is promising in generating reference dynamics for the spread of epidemics in China. 

For this study, an SEIR model is used as the baseline epidemic model; therefore, as discussed earlier, this simulator is only reliable to generate dynamics for the period where no effective government intervention has been implemented. After actions are taken, the SEIR compartments will be invalidated, and more elaborated compartment models should be adopted to account for government's measures such as quarantine and the reportage of suspected cases [36] , as well as the shutdown of transportation at certain places. Another extension to the current model framework is to relax the assumption in the spillover of cross-infection, and instead allow that cross-infection occurs between all routes sharing a finite part in the path. This will increase one search depth in the computation and will be feasible on massive-scale clusters; a super computing device will also facilitate a full-parameter inversion for the Wuhan coronavirus, ideally with real datasets assembled for flowmaps, whose results will undoubtedly uncover further information for the study of this on-going epidemic. Overall, constructed on the multi-layer network flow model with flexible inputs of edge connectivities, flowmaps and arbitrary system parameters, this general-purpose simulator for the city-level spread of epidemics in China has an adaptive nature that could be tuned for specific usage, which might be helpful for policy analysis in emergence response, and earlywarnings of future events.

A multilayer perspective for the analysis of urban transportation systems

Multiscale mobility networks and the spatial spreading of infectious diseases

The structure and dynamics of multilayer networks

The hidden geometry of complex, network-driven contagion phenomena, science

Modeling the multi-layer nature of the European Air Transport Network: Resilience and passengers re-scheduling under random failures

When human networks collide: the degree distributions of hyper-networks

The role of the airline transportation network in the prediction and predictability of global epidemics

The physics of spreading processes in multilayer networks

Mathematical formulation of multilayer networks

Strategies for containing an emerging influenza pandemic in southeast asia

Strategies for mitigating an influenza pandemic

Anatomy and efficiency of urban multimodal mobility

Mitigation strategies for pandemic influenza in the United States

Return of the Coronavirus: 2019-nCoV

Stochastic blockmodels: first steps

First case of 2019 novel coronavirus in the United States

Clinical features of patients infected with 2019 novel coronavirus in

Forecast and control of epidemics in a globalized world

Layered complex networks

Early transmission dynamics in Wuhan, China, of novel coronavirus?infected pneumonia

Transmission dynamics of 2019 novel coronavirus

Early transmissibility assessment of a novel coronavirus in Wuhan

Transport on coupled spatial networks

Networks: an introduction

Growing multiplex networks

Measuring and modeling correlations in multiplex networks

Importation and human-to-human transmission of a novel coronavirus in Vietnam

Large-scale spatial-transmission models of infectious disease

Epidemic spreading on interconnected networks

Congestion induced by the structure of multiplex networks

Business dynamics: systems thinking and modeling for a complex world

An updated estimation of the risk of transmission of the novel coronavirus

Decomposing multilayer transportation networks using complex network analysis: a case study for the Greek aviation network

A novel coronavirus outbreak of global health concern

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study

A compartmental model for the analysis of SARS transmission patterns and outbreak control measures in China

Solving a discrete multimodal transportation network design problem

Phase transitions and optimal algorithms for semi-supervised classifications on graphs: from belief propagation to graph convolution network

A pneumonia outbreak associated with a new coronavirus of probable bat origin

A novel coronavirus from patients with pneumonia in China

Topological relation of layered complex networks

15 level-1 cities, 30 level-2 cities, 70 level-3 cities, 90 level-4 cities and 128 level-5 cities

Time series used in the study are from