Modeling epidemic flow with fluid dynamics

1.   Introduction

Predicting the spatial spread of infectious diseases has been one of the most important tasks for epidemic management and public health administration. A large number of mathematical, statistical and computational models have been developed for this purpose. No doubt that these modeling studies have significantly advanced our understanding in regard to the complex spatiotemporal dynamics of infectious diseases. However, the outcomes of these model predictions have not met our expectation. In particular, the COVID-19 pandemic and its onset, progression and persistence underscore the gap between the complexity of epidemic spreading and our current knowledge in epidemiology ^[1,2,3]. A myriad of mathematical and computational models have been published for COVID-19 ^{[4,5,6,7,8,9,10,11,12,13]}, which have certainly improved our understanding of the transmission and spread of the disease. On the other hand, as pointed out in ^[14], most of the epidemic models based on early COVID-19 data have failed to correctly predict the pandemic progression, often by an order of magnitude.

A common mathematical modeling approach for the spatial dynamics of infectious diseases is based on reaction-diffusion type partial differential equations (PDEs) ^{[15,16,17,18,19,20,21,22,23,24]}. Although these PDE models are analytically important, it is usually a challenge to apply them for practical problems. Particularly, the diffusion coefficients in such models are difficult to calibrate, even if we simply assume that they are constants (independent of time and space). Meanwhile, multi-patch or multi-group models based on ordinary differential equations (ODEs) have also been developed to represent the spatial variations ^{[25,26,27,28,29,30,31]}. These models, however, generally involve a large number of parameters due to the meta-population structure. Consequently, estimating the model parameters and determining their identifiability and sensitivity could be highly nontrivial. With these difficulties, currently available mathematical models can only offer limited insights into the spatial dynamics of infectious diseases regarding when, where and how a real epidemic would spread.

In the present paper, we aim to develop a novel mathematical modeling framework for the transmission and spread of infectious diseases based on principles of fluid dynamics. The essential idea is to model the process of epidemic spreading as a fluid motion, which we refer to as "epidemic flow". This is partly motivated by prior studies on traffic flow where fluid dynamics are incorporated to model the movement of vehicles ^{[32,33,34,35]}. Specifically, the infected/infectious host population is treated as a fluid, while each infected individual is treated as a fluid element. Our study is primarily macroscopic. The details of an individual element are not of our concern, but collectively these elements form a continuum that will represent the macroscopic behavior of the epidemic. In particular, the displacement of the fluid motion would provide useful information on the spatial distribution of the disease, and the velocity of the fluid motion would reveal how fast the epidemic is spreading out. Consequently, through such epidemic flow simulation, we will be able to predict important characteristics associated with the spatial dynamics of the epidemic, which would provide useful guidelines for the design of effective disease control strategies.

Our motivations of modeling epidemic flow are two-fold: 1) There is a long and highly productive history in the modeling and analysis of fluid flow, and mathematical equations (mostly partial differential equations) are well developed and extensively studied in fluid mechanics ^[36,37,38]. 2) Computational fluid dynamics (CFD) has been a very active research area for many decades, and numerous accurate and efficient computational methods are available for solving fluid equations ^[39]. In particular, we will use high-order weighted essentially non-oscillatory (WENO) methods ^[40] to compute our epidemic flow model in the present paper. Hence, we can take advantage of the theory, analytical techniques and computational algorithms in fluid dynamics to improve our understanding in epidemiology, especially in the spatial spread of infectious diseases.

The remainder of this paper proceeds as follows. In Section 2, we present the details of our model formulation for the spatial dynamics of infectious diseases. In Section 3, we describe the numerical treatment of our CFD-based epidemic flow model. In Section 4, we verify the accuracy of our computational methods and, subsequently, apply our model to the COVID-19 outbreak in Wuhan, China as a demonstration of our modeling and simulation framework. In Section 5, we conclude the paper with some discussion.

2.   Model formulation

Our fluid dynamics-based epidemic model, described below, is mainly concerned with the transmission and spread of infectious diseases at the macroscopic level. Our fundamental assumption is that the spatial spread of an epidemic can be represented as an inviscid fluid flow and studied by the theory and computational methods from fluid dynamics.

Let $X$ denote the spatial position. In particular, if we represent a region of our interest in the two-dimensional (2D) space, then $X = (x, y)$ where $x$ and $y$ are the standard horizontal and vertical spatial coordinates. Let $s(t, X)$, $i(t, X)$ and $r(t, X)$ denote the densities of the susceptible, infected (and infectious), and recovered individuals, respectively, at time $t$ and location $X$. The total population density is then given by

Here, $s(t, X)$, $i(t, X)$ and $r(t, X)$ are commonly understood as number densities; i.e., they represent the numbers of susceptible, infected, and recovered individuals per unit area. Since the focus of our study is the macroscopic behavior of the epidemic, we do not attempt to resolve the detailed characteristics and heterogeneities of individual hosts. If we assume that each individual host has approximately the same mass (or, equivalently, we use an average mass to represent each individual), then $s(t, X)$, $i(t, X)$ and $r(t, X)$ may also represent the mass densities of the susceptible, infected, and recovered hosts, subject to a constant multiplication. We will use both interpretations in an interchangeable manner in this work.

Inspired by prior modeling studies on traffic flow, where the movement of vehicles is treated as an inviscid flow described by the Euler equation ^{[32,33,34,35]}, we will model the infected population, with density $i(t, X)$, as an inviscid fluid and the spatial spread of the epidemic as an inviscid fluid flow. The susceptible population, with density $s(t, X)$, is treated as a medium (or, background), through which the fluid moves and interacts. At the initial time $t = 0$, the medium (i.e., the susceptible population) may occupy the entire spatial domain, whereas the fluid (i.e., the infected population) may only cover one or more small areas representing the onset location(s) of the disease outbreak. The early period of an epidemic is typically characterized by an ascending phase. Through the contact between infected and susceptible individuals, more susceptible individuals become infected and the epidemic spreads out to larger areas. This process is modeled as the interaction between the fluid and the medium, through which the medium is continuously transformed into the fluid. Consequently, the density of the fluid and the area occupied by the fluid increase during the ascending phase of the outbreak and keep changing throughout the entire epidemic.

The motion of the fluid depicts the spatial spread of the infection. We let $V(t, X)$ denote the fluid velocity at time $t$ and location $X$ in reference to the medium. This velocity does not necessarily reflect the physical motion of individual hosts; instead, it characterizes the speed and direction of the epidemic movement, which is shaped by the collective behavior and interaction among the human hosts. The velocity field is mathematically represented as a vector; particularly, $V(t, X) = \big(u(t, X), \, v(t, X) \big)$ in the 2D space with horizontal speed $u(t, X)$ and vertical speed $v(t, X)$. Initially the velocity field may be restricted to one or more small areas where the fluid is present. With the interaction and mixing between the susceptible and infected hosts, the velocity field expands and the fluid spreads out.

In general, if $d$ is the density of a physical quantity $q$ in a velocity field $V$, then the continuity equation ^[41,42] states that

where $d V$ represents the flux of $q$ and $\sigma$ is the generation of $q$ per unit volume per unit time. In the present study, we let $\beta$ denote the disease transmission rate, $\gamma$ denote the recovery rate, and $\omega$ denote the disease induced death rate. Applying the continuity Eq (2.2) separately to the susceptible, infected and recovered populations, we obtain

and

Equations (2.3)–(2.5) represent an extension of the basic SIR (susceptible-infected-recovered) model with spatial movement for the infected population relative to the medium. Since the duration of an epidemic is typically short, we ignore the vital dynamics here. Equation (2.5) is usually not needed in the analysis and computation of the model, since the other two equations do not depend on the variable $r$.

Meanwhile, the velocity field $V(t, X)$ is described by the Navier-Stokes equation ^[36,37], the most fundamental equation in fluid dynamics. Based on our inviscid flow assumption, the Navier-Stokes equation is reduced to the Euler equation

where $p$ is the pressure, which is an internal force of the fluid. An equation of state ^[41,43] is needed to close Eqs (2.3)–(2.6) and provide a characterization of the fluid pressure $p$ in terms of other state variables. Here we assume

for a constant $c$. Equation (2.7) implies that the fluid motion is driven by the difference of the fluid density, corresponding to the empirical observation that the spread of an epidemic is often driven by the gradient of the infection level; i.e., it typically spreads from high-prevalence areas to low-prevalence areas.

We make an additional remark for Eq (2.7) by comparing this equation with the ideal gas law ^[39,44]. If a fluid with density $d$ and pressure $P$ can be treated as an ideal gas, then we have

where $R_s$ is the specific gas constant and $T$ is the absolute temperature. If we assume that $T$ is also a constant (i.e., the temperature does not change), then Eq (2.8) is in a form similar to that of Eq (2.7). In this analogue, we may think of infected individuals as moving "particles" that are not subject to interparticle interactions, so that the motion of the infected population may be qualitatively analogous to an ideal gas flow.

A linear stability analysis for this model is presented in the Appendix. In what follows, we will focus our attention on the numerical calculation of the model.

3.   Numerical treatment

Multiplying the vector $V$ to both sides of Eq (2.4) and adding the result to Eq (2.6), we obtain

In 2D, Eqs (2.4) and (3.1) become

We put system (3.2) in a conservation form

where

If we introduce $u_1 = i$, $u_2 = iu$ and $u_3 = iv$, and substitute Eq (2.7) for $p$, we may rewrite Eq (3.4) as

The Jacobian matrices for $E$ and $F$ are given by

and

It can be easily verified that for any real numbers $a$ and $b$, the matrix $\, a E'(U) + b F'(U) \,$ has three real eigenvalues:

and a complete set of eigenvectors. Hence, system (3.3) is hyperbolic.

Since our main focus is on the spatial spread of the fluid, we seek to accurately resolve the spatial dynamics by using a high-order numerical method. To that end, we will apply the fifth-order weighted essentially non-oscillatory (WENO) technique ^[40] from CFD to compute this model. In what follows, we present the details of our numerical treatment. We consider a rectangle $[x_l, \, x_r] \times [y_l, \, y_r]$ as our computational domain. We divide the domain into uniform meshes in the $x$ and $y$ directions, respectively, by $x_i = x_l+(i-0.5)h_x$ and $y_j = y_l+(j-0.5)h_y$, for $i = 1, \, 2, \cdots, \, N_x$ and $j = 1, \, 2, \cdots, \, N_y$, where

We denote the numerical approximation of the variable $U$ at the grid node $(i, j)$ by ${U}_{ij}$. We then formulate the following finite difference scheme to discretize the governing system (3.3):

where $i = 1, \, 2, \cdots, \, N_x$ and $j = 1, \, 2, \cdots, \, N_y$. In this equation, the numerical fluxes $\hat{{E}}_{i \pm 1/2, \, j}$ and $\hat{{F}}_{i, \, j \pm 1/2}$ in the $x$ and $y$ directions, respectively, are computed by the fifth-order WENO method. The essential idea is to convert the conserved variables and circulation into a characteristic space, split the flow fluxes, perform WENO reconstruction, and then transfer the numerical circulation in the characteristic space back to the original space. Below we describe the computational procedure for $\hat{{E}}_{i+1/2, j}$. The treatment of the other numerical fluxes is similar.

$1)$ Calculate an average state variable ${E}_{i+1/2, j}\,$:

Then conduct eigendecomposition for the corresponding Jacobian matrix: ${E}'({U}_{i+1/2, j}) = {R}_{i+1/2, j}{\Lambda}_{i+1/2, j}{R}^{-1}_{i+1/2, j}$.

$2)$ For $i_0 = i-2, \, i-1, \, i, \, i+1, \, i+2, \, i+3$, the original variables and fluxes are transformed into a feature space:

$3)$ Conduct the flow flux splitting. Here we apply the Lax Friedrichs splitting method:

where $\tilde{E}_{l, i_0, j}$ denotes the $l$-th component of $\tilde{{E}}_{i_0, j}$, $\; \alpha_l = \max \limits_{i-2\le i_1\le i+3}|\lambda_l({E}'({U}_{i_1+1/2, j}))|$, and $\lambda_l \big({E}'({U}_{i+1/2, j}) \big)$ denotes the $l$-th eigenvalue of the matrix ${E}'({U}_{i+1/2, j})$.

$4)$ Based on $\tilde{{E}}^{\pm}_{i_0, j}$ calculated in Step 3, $\; \hat{\tilde{{E}}}^{\pm}_{i+1/2, j}$ is reconstructed by the fifth-order WENO scheme. For more details of the WENO reconstruction, the readers are referred to ^[45].

$5)$ Obtain the synthetic numerical flux and convert it to the original space:

In the same way we compute $\hat{{E}}_{i-1/2, j}$, $\hat{{F}}_{i, j+1/2}$ and $\hat{{F}}_{i, j-1/2}$. Then we can evaluate

For the temporal approximation of both $U(t, X)$ and $s(t, X)$, we use the third-order TVD Runge-Kutta scheme ^[46] to advance the numerical solution from the time step $n$ to time step $n+1$:

As a summary of our numerical treatment for the model, a time marching process based on the third-order TVD Runge-Kutta method, described by Eq (3.16), is applied to advance the solution in time. At each time step, the fifth-order WENO scheme represented by Eq (3.15) is applied for the spatial discretization.

4.   Numerical results

4.1.   Accuracy test

To ensure that we accurately resolve the spatial dynamics, we first verify the spatial accuracy of our numerical method by using a simple example. In general, an exact solution to Eq (3.3) is difficult to obtain. We can, however, pick a known function $U^*(t, X) = \big[ i^*(t, X), \, i^*(t, X)u^*(t, X), \, i^*(t, X)v^*(t, X) \big]^T$, and replace Eq (3.3) by

where

Equation (4.1) can be re-written as

Note that, once $i^*(t, X)$ is known, Eq (2.3) can be integrated to find

Obviously, the function $U^*$ is an exact solution of Eq (4.3). For the accuracy test, our numerical methods described before can be applied to the left-hand side of Eq (4.3), while the right-hand side of Eq (4.3) can be computed by using the known function $U^*$.

Specifically, we set our computational domain as $[-\pi, \, \pi]\times[-\pi, \, \pi]$ and pick the following analytical solution:

with a periodic boundary condition. We set $\beta = \gamma = \omega = c = 1$, and choose a time step size

to ensure numerical stability. Using the fifth-order WENO method described in the previous section, we compute the numerical solution to Eq (4.3) up to $t_{end} = 0.5$, and compare with the analytical solution in (4.5) to obtain the numerical error. Standard grid refinement tests are performed to check the order of accuracy for the numerical method. Table 1 displays the $L^2$ error and $L^{\infty}$ error for $i(t, x, y)$ with different mesh sizes. We clearly observe that, as the grid is refined, fully fifth-order spatial accuracy is achieved. Similar results hold for $u(t, x, y)$ and $v(t, x, y)$ and are not presented here.

4.2.   COVID-19 simulation

We now demonstrate a real-world application of our modeling framework by considering the COVID-19 outbreak in Wuhan, China in 2020. It is known that the incubation period of COVID-19 varies between 2 to 14 days, during which time there may not be any symptoms. Those infected individuals, before they are confirmed and isolated, may move around which subsequently leads to the spread of the disease.

For simplicity, we assume that the total population density $n(t, X)$ in Wuhan city is uniform in space and time; i.e., $n(t, X) \equiv {\rm Const}$, since the duration of the epidemic in Wuhan is short and the disease-induced mortality rate is low. Without loss of generality, we normalize this constant total density to $1$; i.e.,

for any time $t$ and space location $X$. Consequently, $s(t, X)$, $i(t, X)$ and $r(t, X)$ can be regarded as the percentages of the susceptible, infected and recovered components, respectively, in the total population density.

Wuhan has an urban area of 1528 square kilometers ^[47]. We represent the city as a square of the dimension 40 km $\times$ 40 km. The spatial domain in our model is thus defined as

There are several different ways to handle the boundary conditions for hyperbolic flows ^[39,48]. For the simulation results presented in this section, a simple extrapolation for the variables is applied on the four boundary edges of the square domain. The Huanan Seafood Market is commonly believed to be a primary source and the main onset location of the epidemic in Wuhan ^[49]. The market is within the Jianghan District, a chief administrative and commercial district located in the central region of Wuhan city. The market has a reported floor space of about 50,000 square meters. Considering the presence of multi-storey buildings, we assume that the lot area covered by the market is 40,000 square meters and, for simplicity, represented by a square of 0.2 km $\times$ 0.2 km in the center of the domain $\Omega$:

We start our model simulation from January 24, 2020, when the coronavirus disease outbreak in Wuhan began to attract global attention ^[7,49]. We set the initial conditions as

with a relatively high infection density $\theta = 0.05$ at the initial time, and

with a low infection density $\alpha = 10^{-4}$ at the initial time. The fluid representing the infected population is assumed to be stationary initially; i.e.,

The fluid motion is started and driven by the difference of the pressure, defined in Eq (2.7), between the main onset location $\Gamma$ and the outside region.

We make an additional comment to justify our initial model setting on January 24, 2020. Although the majority of the initially reported cases in Wuhan were related to the market (the small square $\Gamma$ in our setting), most of these infected people would not stay inside the market for long. We assume that they were already dispersed and uniformly distributed to the entire city by January 24, 2020. The initial distribution level of the infection outside the market, though, was very low, corresponding to $\alpha = 10^{-4}$ in our setting. Given that the total population in Wuhan at that time was about 9 million ^[4], this would lead to an infection number of around 900 in our initial setting based on the value of $\alpha$. On the other hand, the reported cases in Wuhan were widely believed to be underestimates of the outbreak, as there were many undocumented cases. For example, a Science paper ^[7] estimated that about 86% of the total infections were undocumented before January 23, 2020, and about 35% were undocumented afterwards. We note that the number of reported cases on January 24, 2020 was 572 ^[4,49]. If we adjust this number by the undocumented rate of 35%, then the actual cases on January 24, 2020 would likely be $572/(1 - 0.35) = 880$, which is very close to 900 in our model setting.

Our model based on system (3.2) has four parameters: the transmission rate $\beta$, the recovery rate $\gamma$, the disease-induced death rate $\omega$, and the pressure coefficient $c$. We assume that all the four parameters are constants. The first three of these are standard epidemiological parameters and their values can be found from prior modeling studies for the Wuhan epidemic ^[4,50], with an appropriate scaling by the total population size: $\beta = 0.347$ per day, $\gamma = 1/15$ per day, and $\omega = 0.01$ per day. With these parameter values, we can easily verify that the basic reproduction number $\mathcal{R}_0 > 1$, where $\mathcal{R}_0$ is defined in Eq (A.11) in the Appendix.

The pressure constant $c$, on the other hand, is a new parameter that we introduce into our model based on fluid dynamics. We thus seek to find an appropriate value of $c$ such that we can obtain a reasonable match between our simulation result and the real epidemic data. To that end, we compare the numerical result and the reported data using the cumulative cases. Based on our model, the cumulative cases at time $t$ can be computed by

where $N$ is the total population size of Wuhan city (approximately 9 million people), $A(\Omega)$ is the area of the domain $\Omega$ (about 1,600 km$^2$), and $N/A(\Omega)$ is the (unnormalized) total population density. The first integral in Eq (4.12) measures the number of active infections, the second integral measures the number of recovered individuals, and the third integral calculates the number of disease-induced deaths as of time $t$. We then use our model to fit the number of reported cumulative cases in Wuhan from January 24, 2020 to February 10, 2020, for a period of 17 days in total that corresponds to the ascending phase of this epidemic. We use a spatial grid of $N_x \times N_y = 400 \times 400$ to discretize the domain $\Omega$. Through the simulation, we find that when $c = 5.5$, our model output based on Eq (4.12) best fits the reported data. Figures 1 displays the data fitting result with $c = 5.5$.

With the fitted parameter $c = 5.5$, Figure 2 shows the contour plots of $i(t, x, y)$ at several different times: day 1, day 5, day 7 and day 10. We observe that at $t = 1.0$ (i.e., day 1), the epidemic wave has already spread out from the center of the domain, which has the highest concentration of infection at the initial time. The circle formed by the wavefront has a radius around 5 km, and the largest value of $i$ (about $3 \times 10^{-4}$) occurs on the outer wavefront. Due to the symmetry of our domain and initial conditions, the distribution of $i(t, x, y)$ exhibits a radially symmetric pattern. A linear analysis is provided in the Appendix for such a scenario with a radial symmetry. In particular, it is shown that the wavefront generated by the infected fluid propagates in both the positive and negative radial directions. The epidemic wave continues propagating and, at $t = 5.0$, the outer wavefront moving along the positive radial direction has already been very close to the boundary of the domain, whereas the inner wavefront moving along the negative radial direction has apparently reached the center of the domain. The largest value of $i$ becomes about $3.5 \times 10^{-4}$, and this highest infection density occurs in several places that include a center region and four small regions near the boundary, indicating an impact of the nonlinear dynamics on wave propagation. At $t = 7.0$, the largest value of $i$ becomes about $6 \times 10^{-4}$ which is concentrated in a center region. The value of $i$ decreases from the center to the boundary (with the lowest value around $1.5 \times 10^{-4}$), forming a set of consecutive layers that indicate a pattern of mixing between infected and susceptible individuals. At $t = 10.0$, the infection density reaches as large as $1.0 \times 10^{-3}$ and the highest prevalence again occurs in a center area, with a similar (but more regular) layered pattern as that at $t = 7.0$, showing a higher degree of population mixing and interaction as the number of infected individuals increases.

Before we proceed, we provide another piece of evidence to confirm the accuracy of our simulation results. Figure 2 is generated from a spatial grid of $N_x \times N_y = 400 \times 400$. We have also tested the simulation with higher spatial resolution such as $800 \times 800$ and $1200 \times 1200$, and we obtain almost identical results. In particular, Figure 3 displays the contour plots for $i(t, x, y)$ at $t = 1.0$ with two fine grids ($800 \times 800$ and $1200 \times 1200$), where no difference is noticed in comparison with Figure 2(a).

Now we plot the contours of the velocity components $u(t, x, y)$ and $v(t, x, y)$ in Figures 4 and 5, respectively. We see that, at $t = 1.0$, the velocity field is concentrated at a relatively small circular region centered around the origin with a radius about 5 km, corresponding to the infection profile displayed in Figure 2(a). The velocity field has significantly expanded at $t = 5.0$, consistent with the propagation of the wavefront shown in Figure 2(b). At $t = 7.0$, the velocity field has spread out to the entire domain. The highest value of the horizontal velocity $u(t, x, y)$ occurs near the left and right boundaries of the domain, while the highest value of the vertical velocity $v(t, x, y)$ occurs near the top and bottom boundaries of the domain, representing the spread of the outer wavefront. At $t = 10.0$, both the horizontal and vertical velocity components have developed a clear stripe pattern, where the bands of high values are near the boundary and the bands of low values are near the center of the domain.

Moreover, we present the contour plots of $i(t, x, y)$, $u(t, x, y)$ and $v(t, x, y)$ on the final day of our simulation, $t = 17.0$, in Figure 6. We observe that the distribution of the infection becomes more uniform than that at earlier times. The value of $i$ ranges from $8 \times 10^{-4}$ to $1.1 \times 10^{-3}$. Most parts of the domain have a infection density around $1.0 \times 10^{-3}$, indicating that about $0.1\%$ of the total population become active infections. The velocity components display a more regular stripe pattern than that shown in Figures 4(d) and 5(d), with high values near the boundary and low values near the center. We also note that the initial profile for $i(t, x, y)$ at $t = 0$ is discontinuous, and this discontinuity is propagated with time, due to the nature of the hyperbolic system. This can be clearly seen from Figure 2. As time goes on, however, the discontinuity is gradually smoothed out with the increased level of mixing and interaction between the susceptible and infectious populations, as indicated from Figure 6(a) for $t = 17.0$.

Next, we present some additional simulation results to test our model setting. We have assumed that the main onset area of the outbreak is located in the center of the domain with symmetric initial conditions for all the variables. This setting leads to the propagation of the epidemic wave in a symmetric manner from the center to the boundary. Now we consider a scenario where the main onset location is not in the domain center, with $\Gamma$ in Eq (4.8) replaced by

while other parts of the setting remain the same. We then run the simulation and generate the contour plots of $i(t, x, y)$ at $t = 1.0, \, 5.0, \, 7.0$ and $17.0$ in Figure 7. Comparing to the original results in Figures 2 and 6(a), we observe essentially the same pattern, except that the density distribution in Figure 7 is not radially symmetric anymore. Instead, it looks as though the center of symmetry is shifted from the origin $(0, 0)$ to the new center $(2, 2)$ of the region $\widetilde{\Gamma}$.

For simplicity, we have assumed that the initial densities for $i(t, x, y)$ and $s(t, x, y)$ are uniformly distributed outside the region $\Gamma$, described by Eq (4.10). We now test a scenario where the initial conditions for $i(t, x, y)$ and $s(t, x, y)$ are not homogeneous. Specifically, we assume that initial infection level would take a higher value for a region closer to the domain center (i.e., the main onset location) and a lower value for a region further away from the domain center. We define

We then assign the following initial conditions for the infected density in three different regions: $i(0, x, y) = \delta$ for $(x, y) \in \Gamma_1 - \Gamma$, $\; \, i(0, x, y) = \delta* 95\%$ for $(x, y) \in \Gamma_2 - \Gamma_1$, and $\, i(0, x, y) = \delta* 90\%$ for $(x, y) \in \Gamma - \Gamma_2$, with $\delta = 1.09 \times 10^{-4}$. It can be easily verified that under this non-homogeneous setting, the total number of initially infected individuals still adds up to approximately $900$, the same as our original setting. The initial condition for the susceptible density is subsequently given by $s(0, x, y) = 1 - i(0, x, y)$ for $(x, y) \in \Omega - \Gamma$. The simulation results for the contours of $i(t, x, y)$ at $t = 1.0, \, 5.0, \, 7.0$ and $17.0$ are presented in Figure 8. We again observe a very similar pattern, both qualitatively and quantitatively, as that in Figures 2 and 6(a), indicating a reasonable degree of robustness of our computational results with regard to perturbations of the initial setting.

Finally, we validate our simulation results by implementing and comparing with a different model that can also be used to study the transmission and spread of infectious diseases. To that end, we consider a simple reaction-diffusion system where the susceptible population is treated as the background and where the infected population undergoes a diffusion process relative to the background:

The parameter $d_f$ is the diffusion rate and is assumed to be a constant. To compute this model, a third-order Runge-Kutta method is used for the temporal discretization and a fourth-order centered difference scheme is applied for the spatial discretization. When there are more than one diffusion coefficients in a reaction-diffusion system, it is often a challenge to estimate these diffusion rates, and their values are typically prescribed in an application ^[10,13]. Here we seek to estimate the single diffusion constant $d_f$ in Eq (4.15) through a fitting to Wuhan COVID-19 data. We thus use similar initial conditions and parameter values for $\beta$, $\omega$ and $\gamma$ as before, and estimate the value of $d_f$ by fitting the model to the reported data on cumulative cases. The best fit is found when $d_f = 5.0$, and the fitting result is presented in Figure 9. We clearly observe that the quality of fitting for this reaction-diffusion model is not as good as that for our epidemic flow model (compare Figure 9 with Figure 1). In fact, a simple calculation shows that the $L^2$ error of the fitting for the reaction-diffusion model is about three times larger than that for our model. On the other hand, a standard analysis of system (4.15) yields that the critical speed for the traveling wave, which approximates the speed of the epidemic spread, is given by $2 \sqrt{d_f (\beta - \gamma - \omega)} \approx 1.98$ km per day. Comparing this to the approximated infection wave speed from our model $\sqrt{c} \approx 2.34$ km per day, we see the two estimates are reasonably close to each other. However, our model is also capable of providing other useful information, such as the detailed velocity field associated with the epidemic, which is not available from the reaction-diffusion model.

5.   Conclusions and discussion

5.1.   Conclusions

We have presented a new modeling framework for infectious diseases using theory and computational methods from fluid dynamics. The focus is to investigate the spatial spread of the infection. Our study is primarily macroscopic and coarse-grained, classifying the human population as three homogeneous classes (i.e., susceptible, infected/infectious and recovered classes) without paying special attention to the individual heterogeneities. In this sense, our approach is closely related to classical mathematical epidemiology based on compartmental, population-level models. Our major innovation, however, is to introduce fluid dynamics into the epidemic modeling, treating the onset, evolution, and spatial progression of the infection as a fluid flow.

This work represents a pilot effort in epidemic flow modeling. We describe the infected (and infectious) population as an inviscid fluid, whose motion characterizes the spatial spread of the epidemic. We emphasize that we do not need to assume the susceptible population is stationary in this modeling framework. Instead of explicitly describing the movement of susceptible people, we regard the susceptible population as a medium and the infected fluid flow as a motion relative to this medium. Although at the individual level, different persons could have very different and random movement, the underlying assumption of our model is that at the macroscopic level, the spread of an epidemic can be approximated by an inviscid flow that is described by the Euler equation. In order to solve the model system, which consists of strongly nonlinear partial differential equations, we have employed advanced computational methods for inviscid flow, including a high-order WENO scheme in particular. Our work represents a computationally intensive approach for modeling infectious diseases, and the results demonstrate a novel and meaningful application of CFD methods in epidemiology.

5.2.   Discussion

Through the simulation to the COVID-19 outbreak in Wuhan, China, with a focus on the ascending phase of the epidemic, our model generates numerical results that match the reported data (specifically, the cumulative cases) with very good accuracy. More importantly, our numerical findings provide new insight and useful information, which are either unavailable or difficult to obtain from currently existing models, regarding the spatial dynamics of an infectious disease, COVID-19 in particular. These include: 1) Our model can generate detailed predictions on the spatial distributions and severity levels of the infection at any time during the epidemic. Such information can help with the management of the epidemic, especially for directing efforts toward areas with high prevalence levels. 2) Our model can output the velocity profile at any location inside the domain of our interest throughout the course of the epidemic. The velocity field measures the speed and direction of the epidemic movement, and such information can help us to determine where and how fast the disease is spreading, which is critical for the design of effective intervention strategies. 3) The model outcome can be utilized to compare the disease risk at different areas throughout the spatial domain, in terms of both the current situation and potential future development, so that public health administrations could strategically scale resources and efforts for different locations.

An additional advantage of our proposed model is reflected in the fitting of data. Standard compartmental models based on ODE systems have been widely used in mathematical epidemiology, and many statistical techniques are available to estimate their parameter values, making their application to epidemic simulation relatively easy. However, because of their lacking of spatial components, such ODE models can only provide limited knowledge on the spatial spread of an epidemic, especially when heterogeneous populations and environments are concerned. Extension from an ODE system to a reaction-diffusion PDE system is a natural and popular way to incorporate spatial dynamics into epidemic modeling, but the diffusion coefficients (which generally take different values for different population groups and spatial locations) may be difficult to calibrate. On the other hand, extension from a standard ODE system to a meta-population (e.g., multi-patch or multi-group) model typically involves a large number of parameters to be estimated and fitted. These challenges hinder the practical applications of currently available ODE and PDE models. In contrast, our proposed CFD-based epidemic model only requires one single parameter (the pressure constant $c$) to be fitted, in addition to a few epidemiological parameters that are well defined and readily evaluated from standard ODE models. This makes model calibration to real data much easier than that of current spatial epidemic models.

We have utilized a relatively simple setting in our COVID-19 simulation to facilitate the implantation of our computational model. Nevertheless, we still observe rich patterns in the spatial distribution and temporal evolution of the infection profile and the velocity field, indicating the high complexity of the spatiotemporal dynamics associated with the epidemic spread. In more realistic scenarios, the onset of a disease outbreak may not take place exactly in the center of the domain, there may be more than one major onset locations, and the distribution of the population density may be heterogeneous. Consequently, the spatial spread of the epidemic may not be symmetric in general. Nevertheless, our modeling framework and computational techniques can be similarly applied in such situations to investigate the detailed spatiotemporal dynamics of the epidemic spread. Additionally, instead of treating the relative motion of the infected population with respect to the susceptible population, we may consider the movement of both populations and treat them as two different fluid flows. This approach might lead to a higher resolution into the population dynamics of both the susceptible and infected people and potentially yield more insight into the spatiotemporal evolution of the epidemic spread. We plan to pursue these research efforts in a future study.

Acknowledgments

ZC was partially supported by the Fundamental Research Funds for the Central Universities of China. JW was partially supported by the National Science Foundation under grant number 1951345. The authors would like to thank the two anonymous reviewers for their comments that have improved the quality of the original manuscript.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix: Linear stability analysis

It is clear that system (3.3) has a trivial solution

which represents a homogeneous steady state of the model with $s = s_0 = 1$, $i = i_0 = 0$ and $u = v = 0$ everywhere in the domain. This is commonly referred to as a disease-free equilibrium. Note that we have rescaled the total population density to $1$. Evaluating the Jacobian matrices $E'(U)$ from Eq (3.6) and $F'(U)$ from Eq (3.7) at $U_0$, we obtain

Meanwhile, from Eq (3.5) we have

with $G(U_0) = 0$ and the Jacobian matrix given by

Linearizing system (3.3) around the constant solution $U_0$, we obtain

With a variable transformation ${U} - {U}_0 \longrightarrow {U}$, we may re-write system (A.5) as

where the variable $U$ in system (A.6) now represents a small perturbation to the equilibrium solution ${U}_0$.

To proceed, we introduce the ansatz

where $j$ is the imaginary unit satisfying $j^2 = -1$, and $k$ and $m$ are wave numbers associated with the horizontal and vertical spatial directions, respectively. Substituting (A.7) into Eq (A.6), we obtain

where ${I}$ denotes the identify matrix. To ensure a nontrivial solution for $\overline{U}$, we know that $\lambda$ must be an eigenvalue of the matrix

The three eigenvalues associated with this matrix are

Clearly, when $\beta -\gamma -\omega > 0$, the eigenvalues $\lambda_1$, $\lambda_2$ and $\lambda_3$ all have positive real parts, indicating that the disease-free equilibrium is unstable. Consequently, a small number of infections introduced into the system, which represent a small perturbation to the steady state ${U}_0$, would grow in size and spread out to other regions. This implies that if we define

then $\mathcal{R}_0 > 1$ would indicate the persistence and spread of the infection. Note that $\mathcal{R}_0$ defined in Eq (A.11) is the same as the basic reproduction number associated with the standard SIR model based on ordinary differential equations; i.e., Eqs (2.3)–(2.5) with the spatial term $\nabla \cdot (iV)$ removed, where $s_0$ has been normalized to $1$.

Meanwhile, Eq (A.10) shows that the wavefront generated by the infected fluid propagates with a speed $\sqrt{c }\,$ along the positive and negative directions of the unit vector $\phi = \big(\frac{k}{\sqrt{k^2 + m^2}}, \, \frac{m}{\sqrt{k^2 + m^2}} \big)$. Thus, the parameter $c$ not only acts as a pressure coefficient, but also represents the speed of propagation for the infection wave under this linear approximation.

Furthermore, if we assume that the dynamical behavior of the model is radially symmetric in the 2D domain, then the analysis presented above may be simplified and more insight may be gained into the linearized solution. For this special case, the infection would spread in all directions with an equal probability and with the same speed. Thus, the variables will only depend on the time, $t$, and the radial coordinate, $\rho$. We let $V(t, \rho)$ denote the radial velocity for the infected fluid, which is now a scalar.

With the coordinate transformation $(t, x, y) \longrightarrow (t, \rho)$ and the assumption of radial symmetry, Eqs (2.4) and (2.6) become

where we have incorporated Eq (2.7). We assume that $\rho$ is relatively large; i.e., we consider a location relatively far from the origin. Under this approximation, we obtain

where, for convenience of analysis, we have introduced $L = {\rm ln}\, i$ to replace the terms associated with $i$. We now consider the linearization of Eq (A.13) around the disease-free equilibrium with $s = s_0 = 1$ and $V = V_0 = 0$, which yields

From Eq (A.14), it is easy to derive

i.e., $L$ satisfied the standard wave equation in the $(t, \, \rho)$ coordinate system. Consequently, the infection wave would propagate in both the positive and negative radial directions at a speed $\sqrt{c}$. Additionally, the solution of the linear system (A.14) can be easily constructed by using the characteristic lines $\rho \pm \sqrt{c} \, t = {\rm Const}$. It is then clear to see that the solution (i.e., a small perturbation about the disease-free equilibrium) will grow with time if $\beta -\gamma -\omega > 0$ and decrease with time if $\beta -\gamma -\omega < 0$. These results are consistent with the information revealed in Eqs (A.10) and (A.11).