Application of fluid dynamics in modeling the spatial spread of infectious diseases with low mortality rate: A study using MUSCL scheme

2.1 Model assumptions

Here, we outline the assumptions made for this model.

1. The model is primarily concerned with the macroscopic behavior of the epidemic, meaning it looks at the large-scale dynamics rather than the detailed characteristics and heterogeneities of individual hosts.

2. The velocity fields of the fluids are macroscopic and reflect the collective movement of the population densities rather than the motion of individual person.

3. The total population density is assumed to be constant over short time periods due to the low death rate from the infection. This allows the model to ignore changes in total population size in the spatial domain.

4. These fluids interact with each other, where susceptible individuals “convert” into infected individuals through contact ( β 1 , β 2 ), modeled as fluid interaction, continuously changing the densities of the populations over time and space.

5. The susceptible population density initially occupies most of the space, while the infected population density is localized to small areas where the outbreak begins.

6. It is assumed that the disease spreads from high-prevalence areas to low-prevalence areas, analogous to fluid flow moving from regions of high density to low density.

Let us assume we have a control volume where the quantity represented by the density ρ can be generated or consumed as well as flow in and out of the volume with velocity field q ¯ . If there is an additional source/sink term σ (representing generation or consumption per unit volume), the total rate of change of ρ within this volume can be described by the conservation principle [2,4,7]:

where V is the control volume, ℓ is the surface enclosing the control volume V , and n is the outward-pointing unit normal vector on the surface ℓ . Applying the divergence theorem, equation (2.2) becomes:

Equation (2.3) is usually the general form for continuity equation [3], with ρ q ¯ representing the mass flux. Now, applying the one-dimensional version of equation (2.3) to the susceptible, infected, and treated population densities, respectively, together with their velocities, we have:

where u , v , w denote the spatial velocities of the susceptible, infected, and treated fluid motion, respectively. These velocities does not depict the contact rate, rather, take into account how fast the epidemic spreads out spatially. With the following initial conditions: S ( 0 , x ) > 0 , I ( 0 , x ) ≥ 0 , T ( 0 , x ) ≥ 0 , and the Neumann boundary conditions: S ( t , 0 ) = 0 , I ( t , 0 ) = 0 , T ( t , 0 ) = 0 to ∂ S ∂ x ( t , L ) = 0 , ∂ I ∂ x ( t , L ) = 0 , ∂ T ∂ x ( t , L ) = 0 having that t ∈ [ 0 , ∞ ) and L ∈ Ω . This assumes that no individuals enter or leave the boundaries of the domain, meaning that the flux of the susceptible, infected, and treated individuals is zero at the boundaries.

On the other hand, the parameter β 1 denotes the contact rate when infected individuals comes in contact with susceptible individuals, β 2 denotes the contact rate when the infected individuals who are undergoing treatments come in contact with the susceptible persons (note that β 2 < β 1 ), τ represents the ratio of infected individuals who receives treatment, γ denotes the ratio of treated individuals who obtains recovered, while μ represents the disease mortality rate (which is assumed low).

We drew ideas from the work of Cheng and Wang [13,14], which expressed the mathematical form of the ideal gas equation [33,62],

and incorporated it into the Euler’s equation of motion (having that pressure force is the only body force), and characterized the pressure of the ideal gas in terms of the state variables (S,I,T) [1,3,44,45,54]. Note that ℛ is the ideal gas constant, T is the temperature (which is assumed constant), and d is the density. Thus, we can regard the susceptible, infected, and treated individuals as independent “particles” whose movements are unaffected by other particles. On a macroscopic level, the movement of each population density can be likened to the flow of an ideal gas. Additionally, the distribution of fluid density influences the motion of each fluid, as demonstrated by Cheng and Wang [14]. Consequently, the spread of disease will tend to move from areas of high prevalence to those with lower frequency. At a constant temperature and for a given specific gas constant, the one-dimensional velocity field that describes the epidemic flow can then be expressed as [9,16,20,42,64]. The motion of the susceptible, infected, and treated population densities is dictated by the Euler equations from fluid dynamics, where their movement is driven by density gradients, similar to how fluid pressure operates. The velocity field is therefore expressed as

Equation (2.6) represents Euler’s equations of motion for an inviscid (non-viscous) fluid, where u , v , and w are the velocity components corresponding to S , I , and T , respectively. The constants k s , k i , and k t suggest that temperature remains constant throughout the epidemic cycle. It is important to note that equation (2.4) works in tandem with (2.6); the former describes the dynamics of the disease model, while the latter projects the velocity field associated with the compartments S ( t , x ) , I ( t , x ) , and T ( t , x ) , exhibiting fluid-like behavior. Readers interested in how equation (2.6) was derived can refer to [13].

3.1 Spatial treatment

In this section, we will explore the one-dimensional spatial analysis of the model. Let us focus on equations (2.4) and (2.6). By applying the velocity vectors to both the left-hand and right-hand sides of equation (2.4) and then combining the outcomes with equation (2.6), we arrive at the following result:

Then, expressing the terms of equations (2.4) and (3.1) spatially, we have

In its conservative form, equation (3.2) becomes

where

For simplicity, we define m 1 = S , m 2 = S u , m 3 = I , m 4 = I v , m 5 = T , and m 6 = T w . We then rewrite

Computing the Jacobian ( J ) matrix of the flux function H ( U ) , we have

Then, re-substituting m i for i = 1 , … , 6 in equation (3.6), we have

It can be easily verified that the eigenvalues of matrix (3.7) are

Since the eigenvalues are real and distinct, together with the set of eigenvectors, system (3.7) is hyperbolic. Hence, the hyperbolicity indicates the existence of solution and well-posedness of the model equations (2.4)–(2.6) [14].

3.2 Linear analysis

It is often challenging or even impossible to solve nonlinear PDEs analytically. In order to make complex nonlinear PDE problems more manageable for analysis and numerical computation, a usual mathematical strategy is to linearize the PDE. Take for example, if U 0 is an equilibrium state, and U = U 0 + δ (where δ represents small perturbations), then a nonlinear function F ( U ) can be approximated by F ( U ) ≈ F ( U 0 ) + ∂ F ∂ U U 0 δ , where ∂ F ∂ U U 0 is the Jacobian matrix of F evaluated at U 0 . A typical example can be found in [32,58].

Given that the total population density N ( t , x ) is re-scaled to “1” [67], then the equilibrium state otherwise known as the steady state, for equation (3.3) becomes:

As such, it can be ascertained from equation (3.8) that at equilibrium state, S = S 0 = 1 , I = I 0 = 0 , T = T 0 = 0 , and u = v = w = 0 for all x ∈ Ω and t ∈ [ 0 , ∞ ) . This is commonly known as the disease-free equilibrium (DFE) state [14]. To linearize equation (3.3), we begin by evaluating the Jacobian of H ( U ) at the equilibrium state. Thus, we have:

The source term can be approximated as

where ð = U − U 0 denotes a minute perturbation and Q ( U ) ≈ Q ( U 0 ) + ∂ Q ∂ U U 0 = J Q ( U ) ( U 0 ) . Since Q ( U 0 ) = 0 , the Jacobian of Q ( U ) computed from equation (3.5) becomes

Evaluating equation (3.11) at U 0 , we have

Therefore, the linearization of equation (3.3) becomes

where J H = ∂ H ∂ U , J Q = ∂ Q ∂ U , and for convenience, ð = U − U 0 → U . To further proceed, we make an artificial guess of the solution (otherwise called ansatz)

where c is the wave speed in the spatial direction, j is an imaginary unit with the property j 2 = − 1 . Then, differentiating the initial guess (3.14) and substituting the results in equation (3.13), the following result was obtained:

But e η t + j c x ≠ 0 , and in matrix form, η = η ℐ (where ℐ represents the identity matrix). Therefore, (3.11) becomes

For a nontrivial solution of U ¯ , η is expected to be an eigenvalue of the resulting matrix ∂ Q ∂ U ( U 0 ) − ∂ H ∂ U ( U 0 ) represented as

with eigenvalues

The eigenvalues η 1 and η 2 represent the wave-front of the susceptible fluid, with a propagation speed c in the direction defined by the unit vector ± ψ , where ψ = k s k s . Similarly, the eigenvalues η 3 and η 4 describe the wave-front of the infected fluid, also moving at speed c in the same direction ± ψ . Assuming that β 1 − ( τ + μ ) > 0 , the eigenvalues η 3 and η 4 will have positive real parts, while η 5 and η 6 will have negative real parts. This indicates that the DFE is a saddle point and thus unstable [14]. Trajectories may diverge along directions corresponding to eigenvalues with positive real parts, while they may converge along those associated with negative eigenvalues. In the case of disease spread, this might translate to a scenario where:

1. A minor perturbation of the illness in the direction of the positive eigenvalue, such as the introduction of a few new infections, may generate an epidemic or outbreak as the perturbation increases.

2. In contrast, perturbations in the directions corresponding to negative eigenvalues will eventually fade out, resulting in the disease’s final elimination.

To reinforce this perspective, we examine the reproduction number of the corresponding ordinary differential equation (ODE) version of the model [19,63]:

R 0 > 1 means the disease persists in the population, and vice versa. Note that with regard to equation (3.8), S 0 = 1 at DFE.

4.1 Time updating

The third-order Runge-Kutta scheme [8,40,66] for time integration

where 1 ≤ i ≤ Z x , n = 0 , 1 , 2 , … , and L ( U ) denotes the spatial discretization operator.

4.2 Numerical scheme

To spatially resolve the resulting nonlinear PDE (3.3) with MUSCL scheme, we begin by discretizing the domain: we divided the spatial domain into grid points with cell center x i and cell interface x i + 1 2 , where i denotes the indexes of the grid point, thereby letting Δ x denotes the spatial grid spacing, and Δ t denote the time step. At each grid point x i , we let U i 0 = U ( x i , t = 0 ) represent our initialized solution.

Accuracy of the numerical method for spatial discretization of our PDE is the utmost priority. To achieve this, we reconstruct the variable U and compute the left U i + 1 2 L and right U i + 1 2 R states at the interfaces of the cell x i + 1 2 . To ensure the monotonicity of our solution, we used the van Leer slope limiter [56,65], where the slope ratio r i for each cell i and the van Leer limiter ϕ ( r ) are given by

and

Then, considering a given cell i , reconstructing the left and right state at the boundary i + 1 2 , we have

where Δ U i = U i + 1 − U i . Next based on the reconstructed states, we calculate the flux F i + 1 2 at each cell interface using the Harten-Lax-van Leer with contact (HLLC) Riemann solver [5,51]

where H i + 1 2 L = H U i + 1 2 L and H i + 1 2 R = H U i + 1 2 R represent the HLLC of the left and right flux functions. The intermediate state flux H * of the HLLC (4.5) is then calculated with the wave speeds S L , S R , and S M , defined by left speed, right speed, and middle (intermediate) speed, respectively, i.e.,

and

where the left and right velocities were chosen based on the maximum and minimum velocity signals in the system. For simplicity, we chose a unit difference in each of the velocities u R , v R , w R , and u L , v L , w L . For the middle wave speed, we use Rankine-Hugoniot condition [23,28,35]

Let “ a ” denote the sound speed of the system. Therefore, we define a S , a I , and a T to be the sound speed in each compartment, respectively. For the purpose of the system, we approximate each sound speed to be: a S = k s ∂ S ∂ x , a I = k i ∂ I ∂ x , and a T = k t ∂ T ∂ x . Different cases were considered, where ( a S , a I , a T ) = ( 1 , 0 , 0 ) at DFE, and ( a S , a I , a T ) = ( 1 , ε , ε ) at any other equilibrium point, where 0 < ε < 1 . Implementing ( a S , a I , a T ) into equations (4.6) and (4.7), the new left and right speed becomes

and

Process (4.3)–(4.8) and (4.1) are repeated until discretization is exhausted.

Then, for the time step Δ t , solution is updated at cell i by:

Therefore, the steps in applying the numerical method involves discretizing the spatial domain and initializing the solution; use the van Leer limiter to rebuild the variables at cell interfaces; use the HLLC Riemann solver to compute the numerical fluxes at cell interfaces; and use the third-order Runge-Kutta scheme to integrate the solution in time. The given PDE can be successfully solved using the MUSCL scheme with van Leer limiter, HLLC Riemann solver, and third-order Runge-Kutta time integration by following these steps.

In Figure 1, the two plots were analyzed separately at τ = 0.1 and τ = 0.315 . In both instances, it was assumed that the entire population was susceptible (blue curve) prior to the onset of disease transmission. The graphical plots (Figure 1) illustrate the numerical solutions for the population densities of Susceptible (S), Infected (I), and Treated (T) groups under different treatment response scenarios. The figure compares two cases, where the treatment parameters are set at τ = 0.1 and τ = 0.315 . These plots show how the populations change spatially across the domain. In both treatment scenarios, the densities of S, I, and T vary along the spatial domain x , with noticeable differences in infection levels and treatment effectiveness. Notably, as τ increases from 0.1 to 0.315, the infected population density decreases, suggesting that higher treatment rates are more successful in reducing infections. In contrast, the density of the treated population rises, indicating a direct relationship between the treatment rate and the number of individuals receiving treatment.

Figure 1

Population densities over space at different treatment responses ( τ = 0.1 and τ = 0.315 ), respectively.

The surface plot (Figure 2) reveals that with a high treatment rate, the infection (infected population density) fades away over time. Conversely, the treated population density plot shows a positive correlation between the treatment rate and the decline in the infected population density. This surface plot offers a visual depiction of how the susceptible, infected, and treated population densities change over time and space. It mathematically demonstrates how these population densities are affected by spatial movement (advection), interaction terms (such as infection and treatment), and natural processes (recovery and death). This visualization makes it easier to assess the outbreak’s spread, the impact of treatments in controlling the disease, and whether the system reaches a steady state or remains dynamic.

Figure 2

Surface plot of the population density over time and space at τ = 0.315 .

The graphical simulation (Figure 3) of the epidemic disease (2.4)–(2.6) with the corresponding initial conditions, is being analyzed. To this effect, we assigned the parameter values k s = k i = k t = 1 and utilized the same parameter data contained in [49] to simulate. Consequently, we now move forward with the spatial simulation.

Figure 3

(a)–(c) The spatio-temporal dynamics of susceptible, infected, and treated population density of model (2.4) with the parameters: β 1 = 0.5 , β 2 = 0.23 , τ = 0.12 , γ = 0.42 , and μ = 0.0005 . (d)–(f) The velocity wave-front of S , I , and T , while showing the homogeneous distribution of the population level in each compartment. The plots above centers on the endemic equilibrium of the model E * = ( S * , I * , T ) with the basic reproduction number R 0 = 4.6940759123943545 .

The graphical analysis in Figure 3(a)–(c) presents the surface plot for the model equations (2.4)–(2.6), where the system’s transmission rate ( β 1 ) is notably high. In contrast, plots Figure 3(d)–(f) illustrate the flow velocity for each population density, depicting the spatial spread and progression of the disease across different locations ( x ) for compartments S , I , and T . Notably, the treated population density reaches its peak at a specific spatial point in Figure 3(a)–(c). Within a confined space with a constant density population, the spatial dynamics exhibit sharp variations in the endemic state, particularly in the infected and treated density populations Figure 3(d)–(f). In summary, the analysis shows that higher infection rates lead to sustained and widespread infection levels within communities.

The graphical representation in Figure 4(a)--(c) illustrates the relative impact of a low contact rate within the system. This leads to the establishment of a DFE, where the infection is completely eradicated, and the population remains immune to further disease transmission. Lower transmission rates effectively curb the spread of the disease by reducing the likelihood of individuals transitioning from uninfected to infected. As a result, the disease fails to propagate and eventually ceases to affect the population density. In contrast, plots Figure 4(d)–(f) depict the velocity profile and spatial behavior of each compartment across location x . Our findings suggest that successful disease management and eventual eradication depend on maintaining low infection rates among the susceptible density population and keeping key epidemiological factors at consistently low levels. Consequently, the DFE state ( S 0 , 0 , 0 ) is asymptotically stable, with a basic reproduction number R 0 = 0.539419 , which is less than 1.

Figure 4

(a)–(c) The spatio-temporal dynamics of susceptible, infected and treated population densities of model (2.4) with the parameters: β 1 = 0.05 , β 2 = 0.02 , τ = 0.1 , γ = 0.5 , and μ = 0.0005 . (d)–(f) The velocity wave-front of S , I , and T , while showing the homogeneous distribution of the population level in each compartment. The above plots centers on the DFE of the model ( S 0 , I 0 , T 0 ) = ( S 0 , 0 , 0 ) with the basic reproduction number R 0 = 0.4901960784313726 .

It can be seen for the STIR model, the transmission rate ( β ) plays a substantial role in determining the endemic and disease-free dynamics of the system. However, subsidizing the transmission rate and enhancing the treatment parameter is paramount in order to eradicate diseases [39,52].

Figure 5 illustrates the impact of varying treatment parameters τ on the treatment population density. Given the assumption of a short epidemic duration, the figure indicates that as treatment efficiency improves and its objective is met at a specific point in time, it becomes prudent to discontinue the treatment strategy within the system.

Figure 5

Spatial fragmentation of T ( t , x ) at different treatment rates τ .

4.3 Sensitivity analysis

To perform the sensitivity analysis on the basic reproduction number R 0 of the epidemiological model, we need to calculate the partial derivatives of R 0 with respect to each parameter ( β 1 , β 2 , τ , γ , μ ) and evaluate how changes in each of these parameters influence R 0 . Sensitivity analysis quantifies the responsiveness of R 0 to changes in these parameters [37] (Figure 6).

Figure 6

partial rank correlation coefficient plot depicting the sensitivity analysis of the system’s infection cycle.

Therefore, we have

In summary, the spatio-temporal patterns of the velocities reveal that the velocity component u is influenced by the spatial gradient of S ( t , x ) , v is driven by the spatial gradient of I ( t , x ) , while w is modulated by the gradient of T ( t , x ) . In Figure 3(a), due to the steep slope of S ( t , x ) (indicating a rapid spatial variation), the term − k s ∂ S ∂ x becomes large, resulting in significant changes in u over space and time. Conversely, in regions where S ( t , x ) is nearly constant ( ∂ S ∂ x ≈ 0 ), the velocity u remains relatively stable, leading to a more uniform velocity field.

In Figure 3(b), the large gradient driven by a high infection rate β causes rapid variations in I ( t , x ) , leading to substantial shifts in v . Areas with sharp increases or decreases in infection levels will show notable velocity changes, whereas regions with a more uniform infection spread will exhibit a smoother velocity pattern.

For Figure 3(c), it is important to consider that T ( t , x ) is influenced by its transition dynamics and the treatment rate of I ( t , x ) . Regions with higher treatment rates ( τ ) will significantly affect the gradient of T ( t , x ) , thereby impacting the velocity w . Specifically, w will experience large fluctuations in areas where T undergoes rapid transitions, such as at the boundaries of an outbreak. In contrast, w will be more stable in regions where T changes gradually.

4.4 Contribution to knowledge

This work fundamentally connects epidemiology with fluid dynamics, providing a framework that adeptly models the complex spatio-temporal dynamics of epidemics. Advances in population flow modeling, together with high-resolution numerical methods, create a predictive capability that exceeds traditional epidemic models in capturing the real-world intricacies of disease spread. This approach represents a substantial leap forward in epidemic modeling, supporting more adaptable and spatially targeted public health responses.

The system of PDEs used in this context provides an improved mathematical framework for understanding the dynamics of spatial epidemic propagation. By incorporating specific epidemiological characteristics and the effects of spatial mobility, the equations are designed to model the interactions between susceptible (S), infected (I), and treated (T) population densities. Numerical simulations of these complex interactions can be accurately and stably performed by using third-order Runge-Kutta time integration, in combination with advanced numerical schemes, such as the MUSCL scheme with van Leer slope limiter and the HLLC Riemann solver. This work has contributed to knowledge in the following ways:

• MUSCL scheme with van Leer slope limiter: The scheme MUSCL minimizes the numerical diffusion while preserving stability by improving the spatial resolution of the numerical solution. By stopping oscillations close to steep gradients, the van Leer slope limiter improves this even further and guarantees physical correctness.

• HLLC Riemann solver: For accurate computation of fluxes across cell interfaces, the HLLC Riemann solver is utilized. The solution’s shocks and discontinuities are captured, which is important for simulating abrupt changes in population density during the course of an epidemic.

• Third-order Runge-Kutta time integration: The accuracy and computing efficiency of this temporal integration technique are balanced. Even when the PDEs contain nonlinear elements, it permits steady integration over time.