key: cord-0755098-emxz3qyf
authors: Zhuang, Qiao; Wang, Jin
title: A spatial epidemic model with a moving boundary
date: 2021-09-03
journal: Infect Dis Model
DOI: 10.1016/j.idm.2021.08.005
sha: 1758880498081f9d456791f030fab2a98290d0b9
doc_id: 755098
cord_uid: emxz3qyf

We present a new mathematical model to investigate the spatial spread of an infectious disease. The model consists of a nonlinear PDE system with an unknown velocity field, defined on an epidemic domain that changes with time. The moving boundary of the domain represents the wavefront of the epidemic. We conduct an equilibrium analysis to the simplified models represented by ODE systems. We also perform a numerical study on the original PDE system for a range of scenarios, including one under a realistic epidemic setting.

Mathematical modeling is an important theoretical tool in epidemiology. Traditional mathematical epidemiology dates back to Daniel Bernoulli who created a mathematical model for smallpox in 1760 (Hethcote, 2000; Thieme, 2003) . The Kermack-McKendrick and Reed-Frost epidemic models proposed in 1920s made the foundation for compartmental models (Brauer & Castillo-Chavez, 2001; Murray, 2002) . These classical mathematical models for infectious diseases divide the total population into several classes, such as S (susceptible), I (infectious), and R (recovered), each referred to as a compartment, and ordinary differential equations (ODEs) are typically employed to describe the patterns of movement and progression between the compartments.

In order to incorporate the spatial dynamics of epidemics, models based on ODEs have been extended to meta-population frameworks such as multi-group and multi-patch models (Cosner et al., 2009; Hanski, 1999; Levins, 1969) . In such a study, the entire system is divided into a number of parts based on different locations, and the disease dynamics on each location are described by an ODE sub-system. The communications between different locations can be incorporated into each sub-system through a Lagrangian approach (for multi-group models) (Hasibeder & Dye, 1988; Rodriguez & Torres-Sorando, 2001; Ruan, Wang, & Levin, 2006) or an Eulerian approach (for multi-patch models) (Allen, Bolker, Lou, & Nevai, 2007; Arino & van den Driessche, 2003; Gaff & Gross, 2007; Hsieh, van den Driessche, & Wang, 2007) .

Meanwhile, models based on partial differential equations (PDEs) have been introduced to describe the spatial movement of human hosts and the dispersal of pathogens (Allen, Bolker, Lou, & Nevai, 2008; Bertuzzo, Casagrandi, Gatto, Rodriguez-Iturbe, & Rinaldo, 2010; Cosner, 1991, 2003; Magal, Webb, & Wu, 2019; Thieme, 2009; Wang, Gao, & Wang, 2015; Wang & Zhao, 2012; Wu, 2008; Yang and Wang, 2020a) . Most of these PDE systems employ reaction-diffusion equations, representing the random movement of human hosts as a diffusion process. The diffusion rates in these studies are typically fixed as constants. Traveling wave solutions for such PDE models are extensively studied, including the analysis of the thresholds for the existence of a traveling wave, the stability of the traveling wave, and the minimum speed and asymptotic speed (which are usually shown to be equal) of wave propagation (Capasso, 1993; Diekmann & Heesterbeek, 2000; Grenfell and Dobson, 1995; Kopell & Howard, 1981; Kot, 2001; Li, Li, & Hethcote, 2009; Metz and van den Bosch, 1996; Murray, 2002; Rass & Radcliffe, 2003) .

These modeling investigations have certainly advanced our knowledge in regard to the complex spatiotemporal dynamics of infectious diseases. However, many difficulties and challenges remain in the study of the spatial spread of infectious diseases based on current models (Bertozzi, Franco, Mohler, Short, & Sledge, 2020; Brauer, 2017; Brauer, Castillo-Chavez, & Feng, 2019; Chowell, Sattenspiel, Bansal, & Viboud, 2016; Riley, Eames, Isham, Mollison, & Trapman, 2015) . For example, a meta-population model generally involves a great number of parameters to represent the host mobility patterns, the crosstransmission rates, and the communication between different locations. Given the large amount of parameters, determining their identifiability and sensitivity and estimating their values are highly nontrivial for such meta-population ODE systems. Thus far, very few PDE-based epidemic models have been applied to practical disease outbreaks, partly due to the difficulty in calibrating the model parameters with realistic data, even if the parameters are simply assumed to be constants. Meanwhile, almost all these PDE systems are defined on fixed spatial domains with prescribed convection and/or diffusion rates, where the model settings may not represent the practical situations. In particular, the traveling-wave solution may not be a good approximation to the propagation of a realistic epidemic wave.

In this paper, we propose a new model based on a system of nonlinear PDEs to investigate the spatial spread of an epidemic, as an effort to partially address the aforementioned challenges. The most significant feature of this model is that it involves an epidemic domain with a moving boundary, which depicts the epidemic wave front. The epidemic domain is changing with respect to time due to the movement of the human hosts (both infected and uninfected), and the expansion of the domain in space represents the spread of the infection. By explicitly incorporating the evolution of the epidemic domain and the (unknown) speed of the epidemic spatial spread into the model, we hope to gain important insight into the spatial dynamics of the epidemic. In particular, we aim to provide explicit information for practical questions such as where the epidemic is going, how far the epidemic is spreading, and how soon the epidemic is reaching a certain location. Our focus in this work is a disease outbreak that starts at one or more "point sources" or small areas, referred to as the onset locations, and our main goal is to study how fast such a disease would reach a higher proportion of the host population and spread to a larger spatial domain.

Our PDE model is closely related to the familiar SEIR (Susceptible-Exposed-Infectious-Recovered) compartmental model. Based on this standard framework, we emphasize the spatial dynamics by introducing a velocity field, which has to be determined as part of the solution, and assuming that human hosts in each compartment undergo a convection process with the velocity field. Consequently, the expansion of the epidemic domain, which contains all the infected hosts at any time, characterizes the spatial spread of the epidemic. Our model differs from existing reaction-convection or reaction-diffusion type epidemic models in that our epidemic domain and velocity field are not prescribed as constants or known functions. Instead, they need to be computed in the solution process and they both change with time, capable of catching the realistic features of the epidemic spread.

The remainder of this paper is organized as follows. In Section 2, we present the detailed formulation of our spatial epidemic model that involves a moving boundary. In Section 3, we consider some special scenarios of the model that lead to simplified dynamical systems based on ODEs, and conduct a careful equilibrium analysis of the spatially homogeneous stationary solutions. In Section 4, we describe a numerical procedure to efficiently compute the original PDE system, and conduct several numerical tests for model validation and application. Particularly, we perform a proof-of-concept study by applying our model to simulate the spatial spread of the coronavirus disease in Wuhan, China, the first epicenter for COVID-19. Finally, we conclude the paper with some discussion in Section 5.

We aim to model the spatial spread of a disease outbreak starting from a single spot, with a focus on the movement of the epidemic wavefront. The proposed model is primarily based on the conservation laws. Consider a quantity q that is transported in a velocity field V. Let r be the density of the quantity q per unit volume. The general continuity equation (Attard, 2012) states that vr vt þ V,ðrVÞ ¼ s;

(2.1) where rV represents the flux of q and s is the generation of q per unit volume per unit time.

We represent the region of our interest, where an epidemic takes place, in the two-dimensional (2D) space. We focus our attention on an urban area which has a relatively high and uniform population density. Let X ¼ (x, y) represent the spatial location, where x and y are the standard horizontal and vertical spatial coordinates. Let O(t) denote the epidemic domain which contains all the individuals that have been infected by time t. The domain O(t) is expanding in the 2D space as time progresses, due to the movement of the human hosts. We let S(t, X), E(t, X), I(t, X) and R(t, X) denote the densities of the Q. Zhuang, J. Wang Infectious Disease Modelling 6 (2021) 1046e1060 susceptible, exposed, infected and recovered individuals, respectively, within the epidemic domain O(t) at the spatial location X and time t. To describe the human movement, we introduce a velocity field Vðt; XÞ ¼ ðuðt; XÞ; vðt; XÞÞ with horizontal speed u(t, X) and vertical speed v(t, X), and assume that all the human hosts undergo a convective process characterized by the velocity field V(t, X). In addition, we let Zðt; sÞ ¼ À Z x ðt; sÞ; Z y ðt; sÞ Á denote the outer boundary of the epidemic domain at time t, parameterized by the length s, where Z x (t, s) and Z y (t, s) are the horizontal and vertical positions, respectively, of the boundary.

For some infectious diseases, COVID-19 in particular, individuals in the exposed class do not show symptoms but are capable of transmitting the disease (Chan et al., 2020) . We assume in this study that the exposed and infected individuals are both infectious. In other words, the E and I compartments in our model represent the densities of the asymptomatic infectious and symptomatic infectious individuals, respectively. Susceptible individuals contract the disease through contacting exposed and infected individuals, with transmission rates b E and b I , respectively. Exposed individuals become infected (with developed symptoms) after an incubation period of length a À1 . Infected individuals have a disease-induced death rate w, which is assumed to be small, and a recovery rate g.

Applying the continuity equation (2.1) to S, E, I and R, we obtain for 0 s l(t). In order to derive the equations for u(t, X) and v(t, X), we need to introduce additional assumptions about the velocity field. Instead of doing that for the general setting, we will consider a simplified case below.

In this work, we focus our attention on a special case of the 2D setting, where we assume that the human movement is radially symmetrical with respect to the outbreak onset location, so that the infection would spread to all directions with the same speed; i.e., the spatial spread of the epidemic is isotropic. This scenario is illustrated in Fig. 1 . Based on the radial symmetry, the epidemic domain O(t) can be represented by a circular disk. Let L(t) be the radius of the disk O(t) at time t. Then L(t) completely determines the moving boundary of the domain O(t) which depicts the wavefront of the epidemic. Consequently, the densities S, E, I, R and the velocity V all depend on the radial coordinate r and time t, where t ! 0 and 0 r L(t). Since the duration of an epidemic is often short (typically ranging from a few weeks to several months), and since the disease induced mortality is typically very low for most infections in the present age, we assume that the total population density is approximately a constant during the epidemic period. We further normalize and re-scale this total population density to 1; i.e.,

so that S, E, I and R also represent, respectively, the percentages of the susceptible, exposed, infected and recovered components in the total population density.

With the radial symmetry, system (2.2) is reduced to

Adding up the equations in (2.5) and using equation (2.4), we obtain the following equation for the velocity V which characterizes the epidemic flow in our setting: 1 r v vr ðrVÞ ¼ Gðt; rÞ À wIðt; rÞ; 0 < r L; t ! 0:

(2.6)

A trivial case corresponds to V ¼ 0 for 0 < r L and t ! 0, which would reduce the PDE system (2.5) to a normal SEIR model based on ODEs. In what follows we assume that G(t, r) À wI(t, r) s 0, which would hold for a small w > 0 in particular, to ensure a nontrivial velocity field. The kinematic condition applied to the moving boundary L(t) yields dL dt ¼ Vðt; LÞ:

(2.7)

Since L(t) represents the spatial spread of the epidemic, we naturally expect that V(t, L) ! 0, which indicates that the expansion rate of the epidemic domain O(t) (i.e., the domain that contains all the individuals who have been infected) is nonnegative. This is also a reflection of the fact that the number of the cumulative infected cases is always non-decreasing. Finally, at the center r ¼ 0, the radial symmetry implies that Vðt; 0Þ ¼ 0:

(2.8) Equations (2.5)e(2.8) constitute a spatial epidemic model with a moving boundary, which is our focus in this study. We will discuss related initial and boundary conditions in Section 4.

Even though the model (2.5) has only one space variable r, it is a strongly coupled nonlinear PDE system with a moving boundary and an unknown velocity field, which makes mathematical analysis difficult. Nevertheless, we may gain useful insight by examining some simplified models that reduce the original PDE system to ODE systems.

We first consider the scenario when all the time derivatives in system (2.5) vanish; i.e., all the population densities reach a steady state: S ¼ S(r), E ¼ E(r), I ¼ I(r), and R ¼ R(r). Also assume G ¼ G(r). Then system (2.5) becomes 1 r v vr 

Clearly, the velocity field V is also independent of time; i.e., V ¼ V(r), with the boundary condition V(0) ¼ 0 based on equation (2.8). Consequently, equation (3.2) can be integrated to obtain

Equations (3.1), (3.2) and (3.4) form ODE systems which can be easily solved numerically. If we substitute equation (3.2) into system (3.1), we obtain V dS dr

(3.5) Essential properties of system (3.5) are determined by its equilibria, which represent the spatially homogeneous solutions of the system. Standard equilibrium analysis can be conducted to investigate such properties. Before we analyze system (3.5), however, let us present another important scenario of the PDE system (2.5), which is concerned with uniform spatial distribution of the population densities so that they are only dependent on time: 

(3.10)

Comparing systems (3.5) and (3.10), we observe that the two systems share the same (constant) equilibria, which are referred to as the spatially homogeneous stationary solutions. These solutions represent the spatially homogeneous states for system (3.5), and the temporally stationary states for system (3.10). In particular, there is always a unique disease-free equilibrium for the two systems:

ðS; E; I; RÞ ¼ ð1; 0; 0; 0Þ:

(3.11)

We may gain more detailed information on the constant equilibrium solutions of system (3.5) or (3.10) on a special case when G is a constant. Using the next-generation matrix technique (van den Driessche & Watmough, 2002), we can derive the basic reproduction number of (3.10) as

(3.12)

where the two parts represent the contributions from the exposed individuals and infected individuals, respectively, in shaping the transmission risk of the disease. Based on the classical stability theory for epidemic models (van den Driessche & Watmough, 2002), the disease-free equilibrium is stable when R 0 < 1, and is unstable when R 0 > 1. Now we investigate the existence of a non-trivial equilibrium (S*, E*, I*, R*), also referred to as an endemic equilibrium. Setting the right-hand side of each equation in (3.10) to 0 and substituting the expression of B(t), we obtain

aE À ðw þ gÞI À IðG À wIÞ ¼ 0; gI À RðG À wIÞ ¼ 0:

(3.13)

The third equation in (3.13) yields

(3.14)

Adding the first two equations in (3.13) and substituting equation (3.14), we obtain S ¼ aFðIÞ À G wI À G À FðIÞdGðIÞ: We proceed to analyze the roots of equation (3.16), which will determine the equilibrium points of system (3.10). It is straightforward to observe that F(0) ¼ 0 and P(0) ¼ 0. Thus, the two curves F(I) and P(I) have an intersection at I ¼ 0, which corresponds to the disease-free equilibrium. Meanwhile, we have

(3.18)

By differentiating the function P, we have

(3.20)

þ2wF 0 ðIÞ þ wIF 00 ðIÞ:

(3.22) By our assumption, we are considering a small disease induced mortality rate; i.e., w ≪ 1. When w ¼ 0, it can be easily verified that F 0 (I) > 0, F 00 (I) ¼ 0, G 0 (I) < 0, G 00 (I) ¼ 0, and, consequently, P 00 ðIÞ ¼ 2 aþG ½b E F 0 ðIÞ þ b I G 0 ðIÞ < 0. By continuity, when w is sufficiently small, we also have P 00 (I) < 0, so that the curve represented by P(I) is concave. On the other hand, the curve F(I) is always convex since F 00 (I) ¼ 2w/a ! 0 for any w ! 0.

We already know that P(0) ¼ F(0) ¼ 0. Moreover, through direct algebraic manipulation, we can verify that when w < G, we have P(1) < F(1); i.e., when I ¼ 1, the curve P(I) is below the curve F(I). Therefore, the two curves have an intersection at some I* 2 (0, 1) if and only if P 0 (0) > F 0 (0). Due to the different concavities of the two curves, such an intersection must be unique if it exists. Using equations (3.18) and (3.21), we can easily obtain that P 0 ð0Þ > F 0 ð0Þ if and only if R 0 > 1;

(3.23)

where R 0 is defined in equation (3.12). Hence, we conclude that when G is a constant and w is small, there is a unique positive equilibrium (S*, E*, I*, R*) for system (3.10) if and only if R 0 > 1. This equilibrium corresponds to a non-trivial, spatially homogeneous stationary solution for the original PDE system (2.5).

To examine the stability of the endemic equilibrium (S*, E*, I*, R*), we look at the Jacobian matrix of system (3.10). Since the first three equations in (3.10) do not depend on the variable R, we only need to check these three equations in terms of S, E and I. When w ¼ 0, the Jacobian matrix at the endemic equilibrium is given by

Using the Routh-Hurwitz stability criterion and the algebraic relationship in (3.13), we can directly verify that all the three eigenvalues of the characteristic polynomial associated with the matrix J* have negative real parts. This property also holds when w is sufficiently small due to the continuous dependence on the parameter. Hence, the positive equilibrium of system (3.10) is locally asymptotically stable for small w > 0. Therefore, our equilibrium analysis concludes that for the steady-state system (3.5) or the spatially invariant system (3.10), with constant influx G and small disease-induced death rate w, there is a forward transcrtical bifurcation taking place: when R 0 < 1, solutions locally converge to the disease-free equilibrium (1, 0, 0, 0); when R 0 > 1, solutions locally converge to the endemic equilibrium (S*, E*, I*, R*). These two equilibria represent two spatially homogeneous steady states for the PDE system (2.5), one indicating the elimination of the disease and the other indicating the persistence of the disease.

Finally, suppose that when R 0 > 1 the endemic equilibrium is reached at t ¼ t*, then the function B ¼ G À wIdB* becomes a constant. Consequently, equation (3.8) yields Q. Zhuang, J. Wang Infectious Disease Modelling 6 (2021) 1046e1060

(3.25) and equation (3.9) yields

(3.26) which shows that the radius L(t) would increase exponentially with a constant rate 1 2 B * for t ! t*. The same result can be obtained from equations (3.3) and (3.4). Note that equation (3.4) leads to dL=dt ¼ 1 2 B * L.

For the original PDE system (2.5), we will focus on a numerical study. One challenge in dealing with the proposed model is that it involves an unknown moving boundary L(t) which has to be determined as part of the solution procedure. To overcome this difficulty, we make the following coordinate transformation (t, r) / (t, z):

Under this transformation, the original variable domain 0 r L(t) is mapped to a fixed domain 0 z 1. It is straightforward to calculate the partial derivatives under the transformed coordinates:

Hence, the equations in system (2.5) are mapped to the following equations A numerical procedure can be naturally developed with standard methods for PDEs (Thomas, 1995) . Using the initial conditions (see equation (4.10)), equation (4.4) can be solved to obtain the initial profile of the velocity V. At each time step t nþ1 , equation (4.5) is solved to update L(t nþ1 ). Then S, E, I and R are computed from equation (4.3). Note that Vðt; zÞ À zL 0 ðtÞ L ¼ 0 at z ¼ 0 and z ¼ 1;

(4.7)

based on equations (4.5) and (4.6). Hence, no boundary conditions are needed to handle system (4.3). Instead, at z ¼ 0 and z ¼ 1, an ODE system is solved to update S, E, I and R at the two boundary points. Next, equation (4.4) is computed to update the velocity V. This would complete one computational cycle and the procedure repeats at the next time step. Using the mapped equations and the afore-mentioned numerical procedure, we now present several numerical tests to validate our models and analyses presented in previous sections. We first consider a special form of G(t, z) that represents a constant influx for the susceptible population:

Note that since t ¼ t, we will use t to refer to the time variable in the remainder of this section. We have conducted some mathematical analysis in Section 3 on the simplified models based on ODEs. Now we run numerical simulation to the original PDE model to test the different scenarios associated with the spatially homogeneous stationary states. To start the numerical simulation, we set up a small circular disk as the initial epidemic domain and assume uniform distribution of human hosts within this initial domain. We then examine different values of the basic reproduction number R 0 , defined in equation (3.12) for the spatially homogeneous ODE system (3.10), by varying the constant influx G. Specifically, we fix all the other model parameters, and pick different values of G so that R 0 < 1 and R 0 > 1 can be achieved. For each case, we numerically simulate the original PDE model based on the formulation in equations (4.3)-(4.6). Fig. 2 shows a typical scenario for R 0 < 1, where we set G ¼ 1 person per meter per day which leads to R 0 z0:45. We observe that the density of the susceptible individuals quickly approach a steady state at S ¼ 1, while the densities of the exposed, infected and recovered individuals all approach 0. Thus, the solutions converge to the disease-free equilibrium (S, E, I, R) ¼ (1, 0, 0, 0), and this pattern is uniform for all the spatial positions. For illustration, Fig. 2 only displays the result at z ¼ 0.5, but the results are the same for all 0 z 1 in the spatial domain. Fig. 3 shows a typical scenario for R 0 > 1, where we set G ¼ 0.1 person per meter per day which results in R 0 z2:39. We clearly see that after 50 days, all the solution curves level off at a positive steady state, referred to as the endemic equilibrium, where (S, E, I, R) z (0.41, 0.25, 0.20, 0.14). Again we present the result at z ¼ 0.5, but the pattern is identical for all the spatial points. These results are consistent with our analytical predictions in Section 3. In addition, Fig. 4 displays the evolution of the radius L with respect to time t, under the same parameter setting as in Fig. 3 . The exponential growth of L(t) after t > 50, when the endemic equilibrium of the state variables is reached, can be clearly observed. This behavior agrees with our analytical result in equation (3.26). Next, we perform a case study with more practical relevance, to illustrate the potential of this model for real-world applications. At present, COVID-19 remains a global pandemic. The hundreds of millions of infections caused by COVID-19, the worldwide spread of this disease, and its huge impact on public health and the world economy, are all unprecedented. Here we simulate a scenario that could partially represent the COVID-19 outbreak in Wuhan, China, the first epicenter for the coronavirus disease.

To implement the model, we consider the following form of the influx G with both temporal and spatial dependence: Gðt; zÞ ¼ cz ½1 À Sðt; zÞ; t ! 0; 0 z 1: (4.9) By our definition of the epidemic domain, there is no infection outside the domain O(t) at time t. Equation (4.9) assumes that all the individuals outside the epidemic disk O(t) are susceptible with a (normalized) uniform density of 1, and that the increment of the susceptible density is linearly distributed along the radial direction.

The following initial setting is considered. We set January 23, 2020 as day 0, when the Chinese government ordered to lock down Wuhan city. We assume that at t ¼ 0 the infection takes place at a small neighborhood of the origin z ¼ 0, represented by a small circular disk with a radius e > 0. We also assume that a small number of infections are present inside this domain, represented by a (normalized) density q > 0. We further assume that there are no exposed and recovered individuals initially and the human individuals within this domain are uniformly distributed at t ¼ 0. Thus the initial conditions for these variables are given by

Lð0Þ ¼ e; Sð0; zÞ ¼ 1 À q; Eð0; zÞ ¼ 0; Ið0; zÞ ¼ q; Rð0; zÞ ¼ 0:

(4.10) Values of the model parameters are provided in Table 1 . These parameter values are based on the reported data and previous modeling studies for in Wuhan. Note that, however, the two transmission rates b E and b I have been normalized by the total population size in Wuhan City; it is estimated that about 9 million people remained in the city by January 23, 2020 (The government of Wuhan h). The Huanan Seafood Market, a live animal and seafood wholesale market in Wuhan, was widely believed as a primary source and onset location of the epidemic Yang and Wang, 2020b) . The initial epidemic domain in our model would represent the area occupied by this seafood market. The market has a reported floor space of about 50, 000 square meters, though the actual area covered by the market would be less than that considering the presence of multi-storey buildings in the market. Hence, we set the radius of the initial epidemic disk as e z 100 m. It is reported that 55% of the first 425 confirmed cases were linked to this marketplace . Within our initial epidemic domain, we set the percentage of infection as q z 5%.

Wuhan City has an urban area of 1528 square kilometers (Wikipedia: Wuhan. Availab). When L(t) ¼ 22 km, the area of the epidemic domain O(t) would approximately equal that of the entire city of Wuhan. Our simulation result for the time evolution of the radius L(t) is shown in Fig. 5 , where we observe that L(t) reaches the threshold value L* ¼ 22 km at t z 25.6 days. Since day 0 in our model corresponds to January 23, 2020, our result indicates that by February 17, 2020, the epidemic had spread to the entire city of Wuhan. For this epidemic (perhaps also for most other epidemics), the specific day and time that the infection had spread to the entire city were unknown. Nevertheless, based on the reported COVID-19 data (Yang & Wang, 2021 ; The government of Wuhan h), the number of cumulative cases in Wuhan was exponentially increasing for the first 25e28 days (i.e., from January 23 to February 16e19, 2020), and then the rate of increase slowed down afterwards. Meanwhile, according to a study on the spread of COVID-19 to places surrounding Wuhan City (You, Wu, Yang, Liu, & Liu, 2020) , a number of cities, towns and villages near Wuhan started reporting a sharp increase for the number of cases during the period of February 10e19, 2020. These patterns of reported data seems to be consistent with our numerical estimate of the growth of the epidemic domain, which is approximated as a regular circular disk in this study. Fig. 6 shows our simulation results for the number of cumulative cases, in comparison to the daily reported data. Based on our model, the number of cumulative cases at time t is computed as where the first integral represents the number of individuals who are currently infected, the second integral represents the number of individuals who have recovered, and the third integral represents the number of individuals who have died from the infection, at time t. Compared to the reported data, our simulation results generally underestimate the number of cases. Note that we have not employed any data fitting method in this simulation. Another possible reason could be our (simplistic) assumption that the epidemic starts and spreads from a single location. In reality, there were initial cases reported in a number of different places in the Wuhan region (Du, Javan, Nugent, Cowling, & Meyers, 2020; Li et al., 2020; You et al., 2020) , all of which would contribute to the spread of the COVID-19 epidemic. Fig. 7 shows a three-dimensional plot for the velocity field V(t, z) versus time and space. We observe that V generally increases with t and z. Meanwhile, Fig. 8 shows a three-dimensional plot for the density of infected individuals I versus time and space. We observe that the infection level remains significant and increases with time near the center of the epidemic domain (z ¼ 0) which, by our assumption, is the onset site for the outbreak. For the majority of the spatiotemporal domain, however, the infection density (or, equivalently, the percentage of infected individuals among the whole population) appears low, due to the large population size (about 9 million) in the city of Wuhan.

Predicting the spatial spread of an infectious disease has been one of the most important tasks in epidemic management and public health administration. Although a large number of mathematical, statistical and computational models have been developed for this purpose, their applications to realistic epidemic scenarios have not met our expectation. In particular, the outbreak and worldwide spread and persistence of COVID-19 underscore the gap between the complexity of epidemic spreading and our current quantitative understanding in epidemiology.

As a pilot effort, we have presented a new PDE-based epidemic modeling framework that involves a moving boundary, with a focus on the spatial spread of the infectious disease. We have introduced a time-dependent epidemic domain in our model, where its moving boundary represents the wave front of the epidemic that is determined by the speed of the human movement. By modeling the collective human movement as a convection process, our study is primarily macroscopic, treating the human population as one or more homogeneous classes, without paying special attention to the individual heterogeneities. In this sense, our approach is built on the classical compartmental, population-level models of mathematical epidemiology. Our major innovation, however, is to introduce a moving-boundary problem into the epidemic modeling, where both the epidemic domain and the velocity field are unknown and have to be determined as part of the solution. Our proposed model explicitly accounts for the interplay between disease transmission, host movement, and epidemic spread, in order to gain deeper insight into the realistic spread of an epidemic in space and its impact on the epidemic progression and evolution. On the other hand, we have kept our mathematical model sufficiently simple, involving a minimal number of equations and parameters, so that it can be possibly applied to the simulation of realistic epidemics. In this work, we have illustrated a potential application of our model by incorporating the COVID-19 data from Wuhan City in China. Our modeling results in this proof of concept provide useful information on the spatial spread of COVID-19 in Wuhan during the first period of the pandemic.

There are a number of limitations in this pilot study. Our model does not represent the detailed mobility patterns of the human hosts. Practically, the movement of individual hosts could be highly heterogeneous and would be better described by an agent-based model. We have introduced a simple deterministic formulation that human hosts move in the radial direction described by a convective velocity field, which makes the analysis and simulation more tractable than that for an agent-based model. We may extend the current modeling framework by incorporating diffusion terms that represent the random component in the motion of human hosts. On the other hand, the model may be unrealistic by assuming that all individuals at the same location move with the same speed. In particular, severely infected individuals may be hospitalized and isolated and may not be able to freely move. Our focus is the collective behavior of human movement that determines the propagation of the epidemic wave, and the velocity field should be understood as the average speed of the human movement. Moreover, we have not utilized any parameter fitting techniques in our numerical simulation, including the application to the COVID-19 epidemic in Wuhan, which may have impacted the quality of our simulation results in comparison to the reported number of cases. Based on this initial effort of model development, we plan to incorporate data fitting, which is a non-trivial task for any PDE-based epidemic models, in our future efforts to improve the accuracy of our model predictions and to enable more realistic applications. Additionally, we may improve our modeling work by using the more general 2D formulation presented in Section 2, and by incorporating multiple epidemic onset locations into our model setting.

Both authors, Qiao Zhuang and Jin Wang, declare that there is no conflict of interest in this work.

Asymptotic profiles of the steady states for an SIS epidemic patch model

Asymptotic profiles of the steady states for an SIS epidemic reaction-diffusion model. Discrete and Continuous Dynamical Systems

A multi-city epidemic model

The challenges of modeling and forecasting the spread of COVID-19

On spatially explicit models of cholera epidemics

Mathematical models in population biology and epidemiology

The effects of spatial heterogeneity in population dynamics

Mathematical structures of epidemic systems

A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: A study of a family cluster

Mathematical models to characterize early epidemic growth: A review

The effects of human movement on the persistence of vectorborne diseases

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Using the COVID-19 to influenza ratio to estimate early pandemic spread in Wuhan

Modeling tick-borne disease: A metapopulation model

Ecology of infectious diseases in natural populations

Metapopulation ecology

Population dynamics of mosquito-borne disease: Persistence in a completely heterogeneous environment

The mathematics of infectious diseases

Impact of travel between patches for spatial spread of disease

Target pattern and spiral solutions to reaction-diffusion equations with more than one space dimension

Elements of mathematical ecology

Some demographic and genetic consequences of environmental heterogeneity for biological control

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

Periodic traveling waves in SIRS endemic models

On the basic reproduction number of reaction-diffusion epidemic models

Velocities of epidemic spread

Mathematical biology

Spatial deterministic epidemics, mathematical surveys and monographs

Five challenges for spatial epidemic models

Models for infectious diseases in spatially heterogeneous environments

Epidemiological parameter review and comparative dynamics of influenza, respiratory syncytial virus, rhinovirus, human coronvirus, and adenovirus

Mathematics in population biology

Spectral bound and reproduction number for infinite-dimensional population structure and time heterogeneity

Numerical partial differential equations: Finite difference methods

Influence of human behavior on cholera dynamics

Basic reproduction numbers for reaction-diffusion epidemic models

Spatial structure: Partial differential equations models

Basic reproduction numbers for a class of reaction-diffusion epidemic models

Transmission rates and environmental reservoirs for COVID-19 e a modeling study

The authors would like to thank the editor and reviewers for their helpful comments that have improved the original manuscript. This work was partially supported by the National Science Foundation under grant number 1951345.