Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to COVID-19, mathematical analysis, and numerical study.

The outbreak of COVID-19 in 2020 has led to a surge in interest in the research of the mathematical modeling of epidemics. Many of the introduced models are so-called compartmental models, in which the total quantities characterizing a certain system may be decomposed into two (or more) species that are distributed into two (or more) homogeneous units called compartments. We propose herein a formulation of compartmental models based on partial differential equations (PDEs) based on concepts familiar to continuum mechanics, interpreting such models in terms of fundamental equations of balance and compatibility, joined by a constitutive relation. We believe that such an interpretation may be useful to aid understanding and interdisciplinary collaboration. We then proceed to focus on a compartmental PDE model of COVID-19 within the newly-introduced framework, beginning with a detailed derivation and explanation. We then analyze the model mathematically, presenting several results concerning its stability and sensitivity to different parameters. We conclude with a series of numerical simulations to support our findings.

Introduction

Many phenomena in the physical and social sciences feature a compartmental structure, in which the total quantities characterizing a system of interest may be decomposed into two (or more) species that are distributed into (two or more) homogeneous units called compartments. As the system evolves in time, the relative distribution of species across the compartments changes, as different physical conditions alter the species state in each compartment and induce species migration from one compartment to another. Compartmental models have been used extensively in biological, ecological, and chemical applications. Notable examples include the susceptible-infected-recovered (SIR) models and their variants for epidemic modeling [1-3], the Lotka–Volterra models for predator-prey dynamics [1, 4, 5], pharmacokinetic models used extensively in pharmacology [6], and demographic and migration models found in sociology and demography [7-9].

The majority of compartmental models encountered in the literature consist of systems of ordinary differential equations (ODEs). These models, while simple to formulate, analyze, and solve numerically, are limited in their ability to describe dynamics in both space and time. A common strategy to introduce spatial variation into such ODE models is by defining regional compartments corresponding to different areas in physical space, with coupling terms added to the model equations to account for the movement of species among the different regions [10-13]. This approach was recently employed in [14, 15] to model the spread of COVID-19 among the different administrative regions in Italy. While this approach may be effective for some applications, description of complex spatial dynamics within compartments is difficult and possibly even non-feasible in this framework.

In contrast, compartmental models based on partial differential equations (PDEs) incorporate spatial information more naturally. Specifically, PDE models allow for a space-continuous description of the relevant dynamics, enabling one to describe dynamics in time and space across all scales. This represents a significant advantage over ODE models, whose ability to describe spatial information is limited by the number of spatial compartments one includes. Examples of compartmental models based on PDEs can be found in [16-21]. Likely owing to their apparent increased mathematical complexity and more significant computational burden, compartmental PDE models are less common and, to the authors’ knowledge, a systematic study of compartmental PDE models in a general setting has not been performed.

The present work has two primary goals. First, we aim at formalizing PDE compartmental models in a general framework more familiar to continuum mechanics. Accordingly, we reinterpret such models as fundamental equations of balance and compatibility, with the relationship between the balance and compatibility equations defined by a constitutive relation. We believe that such a framework may be useful to researchers seeking to better understand general compartmental models, and may ultimately help facilitate interdisciplinary collaboration. Our second goal is to improve our understanding of a specific compartmental PDE model, which describes the spatiotemporal spread of COVID-19, and has demonstrated good agreement with the measured data [22], from the physical, mathematical, and numerical points of view. As the reinterpretation of the COVID-19 model within the continuum mechanics framework plays a significant role in this study, the stated goals are complementary.

The current work is organized as follows. We begin by introducing compartmental PDE models within the continuum mechanics framework in Sect. 2. As a preliminary example to aid understanding, we derive a simple two-compartment Lotka–Volterra-type model within this setting. In Sect. 3 we turn our attention to the COVID-19 model discussed in [22], beginning with its derivation within the newly-introduced notational system from continuum mechanics. We analyze the model mathematically and establish its formal sensitivity to diffusivity and its stability in the norm. We also use an ODE variant of the model, which does not incorporate diffusion, to define a basic viral reproduction number , which is extensively used as an epidemiological indicator of infectious disease spread. A brief spectral analysis is also performed on the ODE variant. Then, Sect. 4 presents a series of numerical simulations in 1D and 2D to examine different aspects of the COVID-19 model behavior. In 1D, we seek to observe how changes to the spatial and temporal discretization affect the model’s numerical solution. In 2D, we analyze how changes in diffusion affect the physical behavior of the model. For both the 1D and 2D problems, we evaluate the effectiveness of the ODE-derived as a predictor of model behavior, demonstrating the significance of spatial diffusion on modeled viral reproduction. We conclude by summarizing the presented results and suggesting directions for future research in Sect. 5.

General formulation of compartmental models in a continuum mechanics framework

We consider a system which may be decomposed into N distinct species: . Each is a function describing the spatiotemporal distribution of the given species with spatial variable and time variable t. It is often the case that has a natural interpretation: for example, the may represent well-defined subgroups of a given population, with their sum then yielding the total population. However, this does not always hold. For instance, the may describe the populations of different animal species, rendering their summation physically meaningless without additional normalization. It is always the case, however, that the are the fundamental quantities of interest describing the system dynamics, and change in response to some or all of the other species in the model.

We arrange the in a vector in such that . Rather than using the more traditional notation found in mathematical and biological references, we opt here for a general notational convention more common to continuum mechanics. Hence, over a spatial domain and a time interval our equations read:plus appropriate initial and boundary conditions. In the system above, Eq. (1) represents a force balance in terms of an internal force , which is thermodynamically conjugate to , and an externally applied force . Physically, we may interpret as describing the changes in the extensive properties of a given species. Eq. (2) represents the compatibility equation in terms of species and specie gradient . Physically, we may interpret as the specie gradient in space. Then, the relationship between the balance and compatibility equations is defined by the constitutive relation Eq. (3).

The role of the externally applied forces defined in Eq. (4) is fundamental in compartmental models, and warrants some additional discussion. As these forces depend on the unknown variable , often in a nontrivial way, their description as ‘externally applied’ may initially seem inconsistent. To understand why such an interpretation is well-motivated, we recall that Eqs. (1)–(4) describe N different species and their relative distribution in time and space. A species may therefore act on a different species , such that for this represents an external force. These intra-species interactions are described by in Eqs. (1) and (4). We note additionally that may only depend algebraically on . Often, these terms are referred to in the literature as ‘reaction terms’ [17]. We consider Eqs. (1)–(2) to be fundamental and fixed; i.e., any PDE compartmental model will share these equations. The relations in Eqs. (3)–(4) thus define the specific behavior of a given model.

Example: two-compartment Lotka–Volterra-type model

We illustrate the continuum mechanics framework presented above by first considering a two-compartment Lotka–Volterra model, also known as the predator-prey equations. This model describes the interaction between two animal populations, predator and prey, in time and space [5]. Let where represents a population of rabbits (prey) and a population of wolves.1 For ease of dimensional analysis, we denote characteristic population, length and time scales as P, L and T. Our model assumptions are: As the compatibility equation Eq. (2) describes the change in a population resulting from its movement in space, with our constitutive relation Eq. (3) we therefore seek to describe the natural tendency for a given population to move. This tendency to move (or resist movement) can be seen as internal forces that regulate the rate at which movement occurs. Specifically, the source of such forces in the current setting may be the level of exertion required for a member of a population to move a certain distance. Therefore, we consider the following definition for the constitutive relation:2where and are scalar “diffusion” parameters with units . The line above and is to indicate that these are constant, scalar quantities, a convention we will use throughout the present work. The constitutive relation Eqs. (5)–(6) can also be seen as arising from the limit of a probabilistic random walk [23]. That and are scalars (and not tensors, as may be the case in general) results from assumption 1, which implies that movement exhibits no directional preference [23].

The movement of both rabbits and wolves exhibit no spatial preference and are independent of each other;

The food supply for the rabbits is sufficiently plentiful such that it does not depend on rabbit population (in biological terminology, we say there is no intraspecific competition [5]) ;

The wolves have no sources of food other than rabbits;

The mortality rate of the wolves, as well as the non-predation mortality rate of the rabbits, does not depend on population size.

We now define the external forces . Assumption 2 implies that the reproduction rate of rabbits grows with population size without any limiting factor, as their food supply is unconstrained. In mathematical terms, this is expressed as:where is the reproduction rate of the rabbits and has units .

Assumption 3, however, implies that the reproduction rate of wolves is naturally limited by the size of the rabbit population. Accordingly:with the reproduction rate of the wolves a function of the rabbit population r. We consider the simplest possible case and postulate is a linear in r:with . Note that has units , reflecting its dependence on the local rabbit population. We naturally expect, in turn, that the number of rabbits eaten by wolves increases with the number of wolves. Then,where is the predation rate and depends on w. We again assume this function to be linear in w, giving:where has units .

Assumption 4 simply states that the mortality of wolves and the mortality of rabbits has no dependence on the population size of either species. Mathematically, we may write this as :where the mortality rates and are both nonnegative and with units .

From Eqs. (7)–(13), we may define as:The relations in Eqs. (5) and (14) are sufficient to define the model in terms of Eqs. (1)–(4). Written in a notation more common to mathematical biology, the model reads:

Spatiotemporal model of COVID-19 infection spread

We now discuss the COVID-19 model proposed by Viguerie et al. in [22]. We consider a population p of individuals divided into compartments corresponding to disease status, modeling the movement in space and time of the subpopulation in each compartment. Specifically, these compartments are the susceptible population s, the exposed population e, the infected population i, the recovered population r, and the deceased population d. Note that d refers only to deaths due to COVID-19. We denote the living population pool as . Due to the names of the compartments used, this model may be called a susceptible-exposed-infected-recovered-deceased (SEIRD) model. We therefore formulate the problem in terms of the vector containing the different compartments.

Model derivation and explanation

Following the example of the Lotka–Volterra-type model shown in Sect. 2.1, we begin by making several model assumptions: As in the Lotka–Volterra model, the compatibility equation describes the changes in a population due to movement, and the constitutive relation will describe the natural extent to which a given population moves. Assumption 1 above implies that such movement is proportional to the living population size n, while assumption 2 sets the movement of the deceased population to zero. Therefore, the constitutive relation for this model is given by:Note that the , , and have units L T P, in contrast to the units in Eqs. (5)–(6), which were L T. This reflects the population-dependent movement rate implied by our assumption 1. Another way to interpret Eqs. (18)–(19) above is as a heterogeneous diffusion process, where the amount of diffusion is proportional to population size.

Movement is proportional to population size; i.e., more movement occurs within heavily populated regions;

No movement occurs among the deceased population;

There is a latency period between exposure and the development of symptoms;

The probability of contagion increases with population size;

Some portion of exposed persons never develop symptoms, and move directly from the exposed compartment to the recovered compartment (asymptomatic cases);

Both asymptomatic and symptomatic patients are capable of spreading the disease;

All living persons are capable of reproduction (the population is not age-structured);

The non-COVID-19 mortality rate is independent of the population compartment;

New births are susceptible to the virus.

Having quantified the internal forces with the constitutive relation, we now focus on the external forces. Assumption 3 implies that all persons who come into contact with the virus first move to the exposed compartment e from the susceptible compartment s:However, assumption 6 implies that this contact could come from both patients showing symptoms (infected population i) or patients not showing symptoms (exposed population e), and from assumption 4 we conclude that such contact must depend on population size n. Therefore,Color-coding has been introduced for ease of understanding, to clearly demonstrate that any addition to compartment e must be accompanied by an equal subtraction from compartment s. We further assume the functions and to be linear in e and i respectively:Parameters and are called the contact rates (units ) and correspond to the likelihood of contagion resulting from contact with an asymptomatic or symptomatic person, respectively. We now define the function A(n) as:One can naturally see that for , the term , increasing with population as desired. is referred to as the Allee parameter (units P) and has to be carefully selected [5].

From assumption 3 we know that some portion of the exposed population e will become symptomatic after a latency period, and hence move to the infected compartment i:where is a parameter corresponding to the latency (or incubation period) with units . However, from assumption 5, we also know that some portion of the exposed population e will never develop symptoms, moving directly to the recovered compartment r. These are called asymptomatic cases. Therefore:where is the asymptomatic recovery rate with units . In Eqs. (27)–(30), we see again that subtraction from one compartment is coupled with an equal addition to another.

Some portion of infected patients will recover, leading to movement into the recovered compartment r:while others will die, moving into the the deceased compartment d:Parameters and are the symptomatic recovery rate and disease mortality rate respectively, both with units . Finally, assumptions 7 and 9 imply that:such that new births enter into the susceptible compartment s, with the birth rate defined by the parameter (units ). Lastly, assumption 8 states that the deaths that are not due to COVID-19 have no compartmental dependence, implying:with representing the general mortality rate, with units . Similar terms appear in the exposed, infected, and recovered compartments as well. The terms in Eqs. (35) and (36) are not color-coded because they are not accompanied by a corresponding term of opposite sign in a different compartment.

Finally, Eqs. (25)–(36) allow us to define the external forces for this model asWe note that the signs in are reversed when compared to Eqs. (25)–(36), as we have now placed these terms on the external force term of the left hand side of the equilibrium equation (Eq. (1)).

Additionally, the standard formulation of the COVID-19 model in mathematical biology would be :

Mathematical analysis

In this section, we examine four results: the sensitivity equations for the diffusion, the nature of the equilibria of the non-diffusive (space-independent) system, the growth/decay behavior of the total population n and the resulting stability in the norm of Eqs. (39)–(43), and a derivation of the basic viral reproduction number for an ODE variant of Eqs. (39)–(43).

Sensitivity equations for the diffusion

A fundamental difference between a PDE model such as Eqs. (39)–(43) and an ODE model is the presence of diffusion. Understanding the nature in which the model solution depends on diffusion is therefore of critical importance. To quantify the dependence of the solution on the diffusion parameters , we compute the sensitivity equations, to determine the quantitiesWe proceed by applying standard arguments of perturbation analysis to Eq. (1) using the constitutive relation introduced in Eqs. (18)–(19) and the external forces defined as in Eq. (38). For we find the equations:where , , , , the third-order tensor readsand .

For the sake of notation, set . Recalling that , we have that for . From now on, for simplicity, we set .

Notice that the matrix has rows from 3 through 5 null since the entries of are constant. Then, the entries in the rows readfor (in the columns 1,2,3,4), while column 5 is null. This leads to the submatrixwhile all the other entries of the matrix are 0.

Equations (45) are equipped with homogeneous initial and boundary conditions of the same type of the conditions for . The resulting solution then describes the sensitivity of a given point in time and space to vary with changes in a given diffusion coefficient.

Equilibria of the non-diffusive system

An analysis of the equilibria of the non-diffusive (space-independent) system provides guidelines on what to expect for the asymptotic behavior of the solution in time also for our PDE system. The equilibria of the space-independent case are obtained by solving the nonlinear algebraic systemIt is promptly computed that this system has the following solutions:here, are constants depending on the initial conditions. Notice that this solution is not incompatible with the occurrence of n at the denominator of some entries of , as the singularity is eliminated by the multiplication of the term by the terms se or si.

for the equilibrium reads ;

for the equilibrium reads ;

for the equilibrium reads ;

A local stability analysis around those equilibria can be rapidly achieved, as for the eigenvalues of are explicitly computed asThese eigenvalues are all real, so we do not expect an oscillatory behavior of the solution in time. If we have a positive eigenvalue, that means that in absence of diffusion the solution asymptotically diverges in time (as one may expect). For , the equilibrium for the first four components is stable.

For , we may expect an asymptotic behavior with the susceptible, recovered and deceased converging to a steady nontrivial equilibrium, while exposed and infected tend to the depletion of the epidemic.

Even though this analysis is conducted on a linearized problem, extensive numerical simulations with different values of the parameters show that the results are realistic. In Fig. 1 we report some illustrative results (obtained by a python in-house solver using the NumPy library), probing the asymptotic behavior under different values of the parameters. We tested the case (with to have a stable equilibrium) in panel (a), the case in panels (b-d) with different values of the parameters and in panel (e). The specific parameter values for each simulation are also reported in Fig. 1. These values do not represent necessarily realistic scenarios; they are just used to probe the asymptotic behavior of the solution under different conditions. We found the expected equilibria. The comparison among panels (b), (c) and (d) pinpoints the importance of the depensation (Allee) parameter A in the pattern of the evolution of the different populations. Depensation impact is particularly evident when n becomes small. The accurate identification of the parameters is clearly a major issue for the practical use of these models (see e.g. [24]).

Fig. 1

Non-diffusive SEIRD model for different values of the parameters (specified on the right-bottom panel). For all the simulations the initial conditions are Persons. The evolution of the solution is consistent with the predictions of our asymptotic analysis

Growth/decay of the total population and stability

In this section, we examine the behavior of Eqs. (39)–(42). One easily observes from the lack of diffusion in Eq. (43) and final column of in Eq. (38) that d does not influence the dynamics of the system. Adding Eqs. (39)–(42) together, we observe cancellation of all the colored terms except , leaving:We let , , , , , , , and rewrite Eq. (50), yielding:We now multiply Eq. (51) by a test function w and integrate over , giving:Applying the divergence theorem, we obtainFrom the assumption of zero-flux boundary conditions, the boundary term on the right-hand side of Eq. (53) vanishes. Note that is the outward-pointing normal vector on and is not to be confused with n. We now let globally, causing the third term on the right-hand side of Eq. (53) to vanish as well, leaving:We define the total population N(t) and total infected population I(t) as:respectively. Then, Eq. (54) can be rewritten as the ODE:whose solution is:which describes the growth/decay behavior of the total population N. This also amounts to an stability result for the system, assuming (as is the case in the present applications).

Remark

If one assumes that the diffusivities are constant and all equal, Eq. (50) reduces further to:The above suggests that one may interpret the global behavior of the system as a nonlinear continuity equation for n transported over the convective field . It can also be interpreted as a reaction–diffusion equation.

Determination of

The basic viral reproduction number serves an important role in the discussion of SIR-type models. In a wholly susceptible population, describes the average number of additional infections caused by each infected individual. Naturally, implies growth of the epidemic, whereas implies decay in infectious spread [5]. The concept of is well-defined for ODE models. However, its extension to a PDE model is unclear, owing to the influence of diffusion. We derive for the ODE version of the PDE model given by Eqs. (39)–(43) and will evaluate its efficacy with numerical tests in Sect. 4.

The ODE version of the COVID-19 model reads:Here, we denote the time derivatives with dots, as we now consider the derivative of a function of a single variable, rather than partial derivatives as done previously in this work. For simplicity, we are not considering non-COVID19 deaths, new births, and the Allee term (hence, ; although their inclusion is not a problem for the analysis shown here.

We proceed using the next-generation matrix procedure outlined in [25]. This approach considers all compartments regarded as ‘diseased’ in a given model. ‘Diseased’ in this context means groups capable of transmitting the infection to others. The terms in the model corresponding to new diseased cases are grouped into a matrix , while the terms describing the movement of existing diseased cases into different compartments are grouped into a matrix . The basic reproduction number is then obtained as the spectral radius of . The justification for this is based on the Perron-Frobenius theorem and is not straightforward. The interested reader is referred to [25].

In our model, there are two compartments that we consider ‘diseased’: the exposed and the infected compartments. Thus, we consider the equation:As stated above, is the matrix containing the new appearances of diseased patients into any compartment, and contains terms which transfer already diseased individuals from one compartment to another. In this case, we stress that movement from e to i is due to the matrix , as an exposed patient moving to the infected category is not considered a new entry into the ’diseased’ category and hence does not participate in . Thus, we define and as:A simple computation shows:which in turn yields:Hence,Applied directly to the ODE model Eqs. (60)–(64), the above Eq. (69) will provide an indication of viral growth rate, as intended. However, given that Eq. (69) does not account for the diffusion present in Eqs. (39)–(43), its effectiveness as an indicator of viral reproduction for the PDE model is unclear and will be examined during the numerical simulations.

Numerical simulations

In this section, we present two numerical simulation studies of the COVID-19 model in 1D and 2D, respectively, to examine the behavior of the model in Eqs. (39)–(43) in detail.

1D simulation study

In this section, we perform a series of simulations using a one-dimensional version of the model in Eqs. (39)–(43). We aim at examining the impact of various numerical solution techniques. In particular, we analyze the spatial and temporal convergence of the computed solutions over various discretization schemes. We also examine the model dynamics more generally and evaluate the efficacy of the definition Eq.  (69) for the PDE model.

Problem setup

We consider the spatial domain given by [0, L] and a time interval [0, T], with days. In the simulations presented in this section, we normalize in space with respect to the characteristic length L of the spatial domain. Hence, we denote . The domain is populated with a population distribution with the unit “Persons.” One may interpret it as denoting a generic normalized population, as we have done with the length scale. The units and values for the relevant space-normalized parameters for the simulations are accordingly presented in Table 1.

Table 1

Parameter values for the 1D simulations

Parameter	Units	Value
	Days^-1	1/8
	Persons^-1· Days^-1	1/2
	Persons^-1· Days^-1	1/2
	Days^-1	1/24
	Days^-1	1/6
	Days^-1	1/160
μ	Days^-1	0
α	Days^-1	0
ν¯s∗	Persons^-1· Days ^-1	5·10^-5
ν¯e∗	Persons^-1· Days ^-1	1·10^-3
ν¯i∗	Persons^-1· Days ^-1	1·10^-10
ν¯r∗	Persons^-1· Days ^-1	5·10^-5

Note all values have been normalized in space by a characteristic length scale L, with this normalization reflected in the units

For the initial conditions, we set and as follows

Figure 2 shows these initial conditions. We further set , , and . Qualitatively, these initial conditions represent a large population center around with no exposed persons and a small population center around with some exposed individuals. We also enforce homogeneous Neumann boundary conditions at and a zero-population Dirichlet boundary condition at for all model compartments. The latter represents a non-populated area at .

Fig. 2

Initial values for susceptible compartment and exposed compartment for the 1D simulations

Additionally, to assess mesh and time integration convergence, we will analyze the total infected population I(t), defined previously as Eq. (56), and the analogously-defined total deceased population D(t). We will also study the time evolution of the total susceptible population S(t), the total exposed population E(t) and the total recovered population R(t), all defined analogously to I(t).

Numerical methods

We use linear finite elements to discretize the spatial domain and we integrate in time using either a second-order implicit (BDF2) or first-order implicit Backward Euler scheme. Each time step is solved fully implicitly using a Picard linearization. All linear systems are solved using GMRES with a Jacobi preconditioner. We employ mass-lumping on all reaction terms.

Mesh convergence

In this analysis, we compare numerical solutions computed on successively refined uniform grids with mesh size =1/500, 1/1000, 1/2000, and 1/4000. Time integration in this study is performed exclusively with a BDF2 scheme using a constant time step days.

In Table 2, we assess mesh convergence using the peak infection date , the peak total infected population , and the final total deceased population D(T). As the peak infection date for 1/4000 is 118 days, we also evaluate I(118) for each level of spatial resolution. We observe a steady increase in all these metrics as is refined and they all progressively approach the corresponding result for the finest mesh. In particular, the quantities reported in Table 2 vary less than 1% between =1/2000 and =1/4000, which suggests a good level of spatial convergence for =1/2000.

Table 2

Mesh convergence of 1D simulations in terms of the peak infection date , the peak total infected population , the total infected population at peak date of the finest mesh I(118) , and the final total deceased population D(T)

Δx^∗	t^	I(t^)	I(118)	D(T)
1/500	122	.038401	.037923	.01265
1/1000	119	.038556	.038482	.012804
1/2000	119	.038667	.038662	.012875
1/4000	118	.038738	.038738	.012910

The relative difference of all these metrics between the cases =1/2000 and =1/4000 is inferior to 1%

Figure 3a–e show plots of the total populations S(t), E(t), I(t), R(t), and D(t) for all the mesh sizes considered in this study. Additionally, Figs. 4, 5, and 6 respectively present plots of s(x, t), i(x, t) and d(x, t) for the different spatial resolutions. Qualitatively, these plots confirm the existence of mesh convergence, as the difference in the plotted variables progressively reduces as we refine the mesh. Indeed, the change between the results for and cases is negligible.

Fig. 3

Mesh convergence analysis in the 1D simulation study. a Total susceptible population S(t). b Total exposed population E(t). c Total infected population I(t) . d Total recovered population R(t). e Total deceased population D(t). f Percent change in norm with successive refinement. These plots show evidence of mesh convergence, with the solutions for =1/2000 and =1/4000 showing minimal differences

Fig. 4

Evolution of the susceptible population compartment over time for varying mesh sizes in 1D. a days. b days. c days. d days. e days. f days. We see similar results across the different meshes, with some noticeable transient discrepancy occurring at t=90 and t=110 days. This indicates that the coarser mesh resolutions cause dispersion error, in which the phase of the solution is affected. In this instance, the solution on the coarse meshes appears delayed

Fig. 5

Evolution of the infected population compartment over time for varying mesh sizes in 1D. a days. b days. c days. d days. e days. f days. We see noticeable transient discrepancy occurring at t=90, t=110, and days, again suggesting dispersion error arising from the coarse discretizations

Fig. 6

Evolution of the deceased population compartment over time for varying mesh sizes in 1D. a days. b days. c days. d days. e days. f days. We see similar results across the different meshes, with some noticeable transient discrepancy occurring at t=90 and t=110 days, where once again the dispersion error on the coarse meshes is apparent

We further assess mesh convergence with an operator that evaluates the percent change in norm for each model compartment c when a given mesh resolution is refined by a factor of r:In Fig. 3f, we plot the values of this operator for all compartments and =1/500, 1/1000, 1/2000 (note the refinement ratio r=2 for all cases). Again, we observe good evidence of mesh convergence, as is notably smaller than and for all compartments c in the model.

An interesting phenomenon we observe is that the largest source of error does not seem to come from over-diffusion or an underestimation of peaks. In fact, peak quantities are predicted similarly across schemes with only slight variation; instead, dispersion error, in which the primary source of error is not the magnitude but instead the phase of the solution, seems the largest problem here. This is particularly apparent looking at Fig. 3, where the cases of appear similar to the more refined simulations, but with a delay in their occurrence. This is further supported by the predictions of shown in Table 2. Referring to Figs. 4c, d, 5c, d, and 6c, d, one may see this effect in time across various compartments.

Temporal convergence

In this analysis, we examine the impact of time integration and time-step size on the numerical approximation of the model solution. We consider both the Backward Euler and BDF2 time integration schemes with time step sizes =0.25, 0.125, and 0.0625 days. As the results in Sect. 4.1.3 suggested =1/2000 was a sufficiently fine spatial discretization, we utilize this mesh resolution here. Table 3 reports the peak infection day , the peak total infection population , and final total deceased population D(T) for each and time integration scheme. As we reduce , these quantities slightly vary for the Backward Euler scheme, while the changes are negligible for the BDF2 schemes. Additionally, we plot the time evolution of the total population in each model compartment in Fig. 7 for all time steps considered in this analysis and for both time integration algorithms. These plots also show that the results for the Backward Euler method exhibit small but perceptible difference, while the solutions obtained with the BDF2 scheme are virtually the same for all time steps.

Table 3

Temporal convergence of 1D simulations in terms of the peak infection date , the peak total infected population , and the final total deceased population D(T)

Δt	Scheme	t^	I(t^)	D(T)
0.25	Backward Euler	116	.0384732	.0129437
0.125	Backward Euler	118	.0385545	.0129038
0.0625	Backward Euler	118	.0386006	.0128836
0.25	BDF2	119	.0386668	.0128755
0.125	BDF2	119	.0386587	.0128688
0.0625	BDF2	119	.0386536	.0128658

As we reduce , the selected metrics show a slight variation for the Backward Euler method, while the changes are negligible for the BDF2 scheme

Fig. 7

Temporal convergence analysis in the 1D simulation study. a Total susceptible population S(t). b Total exposed population E(t). c Total infected population I(t) . d Total recovered population R(t). e Total deceased population D(t). The model solutions obtained with the Backward Euler method change appreciably when the time step is reduced. In contrast, the BDF2 solutions appear well-resolved in time and change minimally as we refine the time step

We also define relevant error quantities for a compartment c using a related but distinct notation to Eq. (72). Some adjustments must be made owing to the fact that we now consider not only one point of comparison as before (the spatial resolution resolution in the case of Eq. (72)) but two (both time step size and time integration scheme). For a compartment c, this quantity reads:where c gives the compartment, the time step, and SCHEME the time integration scheme. In all instances, we compare with the case computed using BDF2 with . So, for example, to quantify relative error of the solution of c computed with the Backward-Euler scheme (BE) using a given , is defined as:Figure 8 plots the temporal convergence in terms of Eq. (74). We observe that the norm difference decreases as is refined for all model compartments, which indicates temporal convergence. The solutions obtained with the Backward Euler method differ noticeably at the coarser time steps, but this difference reduces as we refine . In contrast, BDF2 appears to be well-resolved in time even at the coarsest time step , with the refinement to showing minimal decrease in norm. Thus, the BDF2 scheme provides satisfactory time resolution, even for large time steps.

Fig. 8

Percent difference in the norm between 1D solutions obtained with the Backward Euler (dashed lines) and BDF2 methods (dotted lines) for each . All cases are compared to the BDF2 solution with , with the formal of definition in Eq. (74). The decreasing trends in both plots show temporal convergence. The BDF2 appears well-resolved in time for even the coarsest time step =.25 days. The Backward Euler method requires a fine time step to render results with comparable accuracy to the BDF2 scheme.

Model dynamics

This analysis focuses on the general model behavior, which we examine in a simulation using the BDF2 scheme with =.0625 days, =1/2000. The results for all model compartments are shown in Fig. 9. The infection begins localized in a small population center around and remains localized for the first part of the simulation. At day , the virus reaches the large population center at , and the number of infections begins to increase dramatically. By day , nearly all of the population near has been exposed to the virus. Eventually, due to lack of susceptible individuals, the virus spread ceases.

Fig. 9

Evolution of all model compartments and as defined by Eq. (69) in time and space in 1D. At days (a), we see an initial exposed population centered around . As time progresses to days (b), the outbreak around has grown, with increasing numbers of infected, recovered, and deceased individuals in that region. By days (c), we begin to see the infection reach the large population center around , and by days (d), the outbreak severity in the areas around and are similar. By days, the outbreak around has died down, with the area around now the most affected region; owever, the around indicates that the epidemic may begin to subside. This is indeed the case, and by days (f), we see decreases in infections and increases in recoveries near

In Fig. 10, we compare as defined by Eq. (69) with the exposed and infected compartments. Although the definition in Eq. (69) does not account for diffusion, we observe that still predicts model behavior reasonably well, with the point where corresponding almost exactly with the decrease in new exposures. This is further corroborated by the results depicted in Fig. 9, where the regions where at ultimately show very little contagion, and indeed a distinct ‘hitch’ forms in the distribution infections between the two population centers. Although there is some slight discrepancy owing to the diffusion, we find the definition of given by (69) to be a reliable predictor of the viral behavior for this 1D simulation scenario.

Fig. 10

Evolution in time of as defined by Eq. (69) as well as the total exposed and total infected populations in 1D. We see that is in good agreement with the observed model dynamics, with the decrease of new exposures corresponding nearly exactly to the point where (indicated with the dotted horizontal and vertical lines for ease of visualization). The presence of diffusion, not accounted for in Eq. (69), is likely the source of the slight discrepancy

The model dynamics shown in the 1D simulation in Fig. 9 is similar to that shown for Lombardy in [22] and in the following section. Indeed, for sustained spread of the disease, a certain level of population density is required. Although the disease contagion will diffuse through low-density regions, the growth in those areas tends to be small. Though there have been some notable exceptions, this behavior pattern is similar to what has been observed worldwide, where low population-density regions have largely avoided the catastrophic contagion found in high-density areas [26] .

2D simulation study

The primary difference between the PDE version and the ODE version of the COVID-19 model lies in the influence of the diffusive term. The impact of diffusion on disease spread is a priori difficult to quantify. Increased diffusion leads to a faster and wider dispersion of the virus. However, it also has regularizing effects and may reduce peaks in general. Therefore, exploring such dynamics in detail is important for a full understanding of the model.

In this section, we examine the role of diffusion using the Italian region of Lombardy as our test geometry, using both qualitative analysis and the formally derived sensitivity equation shown in Eq. (45). The problem configuration is identical to the one given in [22] for the simulation scenario labelled ‘Global Reopening B’. This simulation is intended to model the spread of the COVID-19 epidemic in Lombardy, beginning on February 27, accounting for various governmental restrictions and relaxations as they occur. We report the relevant parameter values in Table 4. Good agreement between the presented simulation setup and measured data was shown in [22], and we refer the reader to [22] for a detailed comparison of simulated and measured data. The problem was solved using linear finite elements on an unstructured triangular mesh. The time integration was performed with a Backward Euler scheme, with a Picard-type linearization used to solve the nonlinear system at each time step. All linear systems were solved with GMRES using a Jacobi preconditioner.

Table 4

Parameter values for the 2D Lombardy simulations

Parameter	Units	Feb.27-Mar.9	Mar.9-22	Mar.22-28	Mar.28-May3	May3-
	Days^-1	1/7	1/7	1/7	1/7	1/7
	Persons^-1· Days^-1	3.3·10^-4	8.5·10^-5	6.275·10^-5	4.125·10^-5	6.6·10^-5
	Persons^-1· Days^-1	3.3·10^-4	8.5·10^-5	6.275·10^-5	4.125·10^-5	6.6·10^-5
	Days^-1	1/24	1/24	1/24	1/24	1/24
	Days^-1	1/6	1/6	1/6	1/6	1/6
	Days^-1	1/160	1/160	1/160	1/160	1/160
ν¯_s	km^2· Persons^-1· Days^-1	4.35·10^-2	1.98·10^-2	0.9·10^-2	0.75·10^-2	2.175·10^-2
ν¯_e	km^2· Persons^-1· Days^-1	4.35·10^-2	1.98·10^-2	0.9·10^-2	0.75·10^-2	2.175·10^-2
ν¯_i	km^2· Persons^-1· Days^-1	1.0·10^-4	1.0·10^-4	1.0·10^-4	1.0·10^-4	1.0·10^-4
ν¯_r	km^2· Persons^-1· Days^-1	4.35·10^-2	1.98·10^-2	0.9·10^-2	0.75·10^-2	2.175·10^-2
A¯	Persons	1.0·10^3	1.0·10^3	1.0·10^3	1.0·10^3	1.0·10^3

The values change with date as these correspond to various restrictions (or relaxtions) taken by the government during the epidemic. We note that these parameters are not normalized in space

In addition to the simulation shown in [22], we now examine two additional cases: one in which the values of , , and are doubled, and another in which they are halved. We also consider a case in which , , and are doubled but is halved. This is similar to the parameter setup in the 1D simulations. The main motivation is to avoid possibly nonphysical diffusion among the susceptible population, causing reduced population density in general.

Figures 11 and 12 show the spatial distribution of infected individuals at t=14 and t=30 days, respectively. We see that larger diffusion leads to a wider geographic range of affected areas. This is particularly noticeable in the southeastern clusters in Fig. 11. There, the double-diffusion case produces a homogeneous, continuous region of infection. In contrast, the half-diffusion case shows more localized dynamics, and a clear separation into distinct regions. This separation is maintained in Fig. 12 at t=30 days, whereas the baseline and double-diffusion cases predict a single, larger area of infection. The simulation case in which , , and are doubled but is halved produces intermediate results between the double-diffusion and half-diffusion cases. In all cases, we note that the outbreak path follows regions of high population density, which is the expected behavior given the constitutive relation defined by Eqs. (18)–(19).

Fig. 11

Spatial distribution of the infected population at days in the 2D simulations over the Italian region of Lombardy. a Baseline scenario. b Half-diffusion case. c Double-diffusion case. d Simulation case in which the baseline , , and are doubled but is halved. With halved diffusion (b), we see that outbreaks are more severe, but also concentrated in smaller regions, particularly apparent in the southwest. In contrast, increased diffusion (c) show a less intense peak over a greater overall area. In d, where the diffusion among susceptibles is decreased but increased in other compartments, outbreak severity seems similar to the baseline in a, but covering a slightly larger area (again, most apparent in the southwest)

Fig. 12

Spatial distribution of the infected population at days in the 2D simulations over the Italian region of Lombardy. a Baseline scenario. b Half-diffusion case. c Double-diffusion case. d Simulation case in which the baseline , , and are doubled but is halved. In b, we see both increased severity and interesting localization dynamics; in a, c, and d there appear to be three primary epicenters of infection, while in the case of b there appear to be four. The outbreak in c is much less severe than the other cases, owing to the increased diffusion. In the case of d, we see a larger overall infected area and similar intensity of infection to the baseline (a)

In Fig. 13, we plot the time evolution of the total active infections I(t) throughout the entire region of Lombardy. Both the baseline and half-diffusion simulations show a distinct long-term growth trend that is not observed in the double-diffusion case. While it is tempting to say that increased diffusion leads to reduced outbreak severity, the reality is more complex. Indeed, the case in which is halved while , , and are doubled shows a higher peak and slightly faster growth when compared to the baseline simulation, although the long-term growth more closely resembles the baseline than either the double-diffusion and half-diffusion cases. This makes intuitive sense, as the low diffusion among the susceptible population leads to higher population densities and more contagion, while increasing diffusivity among the exposed and infected compartments accelerates the speed and area of propagation. As discussed in [22], the spatial pattern predicted by the heterogeneous diffusion shows generally good agreement with reality; however, nonlocal transmission is not possible using the model given by Eqs. 39-43 and the addition of nonlocal operators, such as fractional diffusion operators [27], is an area for future development.

Fig. 13

Time evolution of the total active infections I(t) throughout the entire region of Lombardy, showing that diffusion has a strong influence on the dynamics of disease infection. The double-diffusion case has a distinctly different qualitative pattern, with no substantial increase after , while the baseline and half-diffusion cases increase significantly. The dynamics of I(t) for the case in which is halved while , , and are doubled suggests that varying each of these diffusion parameters may induce dramatically different changes in the evolution of the outbreak. In the particular scenario considered here, the number of total infected cases grows slighlty faster and has a higher peak when compared to the baseline case

In Fig. 14, we plot the computed sensitivity parameters found using Eq. (45) at days. The shown plots quantify sensitivity to (left) and (right). The regions most sensitive to are regions currently affected by the outbreak. However, the sensitivity to shows larger values in more highly-populated areas. At the time shown, the area around Milan (in the west of the shown region, the most populated area of Lombardy) was not experiencing a large outbreak in cases. This is reflected in its relatively low sensitivity to . However, its high sensitivity to indicates its vulnerability, irrespective of its current outbreak status. Indeed, the Milan area was ultimately heavily impacted by the epidemic [28].

Fig. 14

Sensitivity of the computed baseline solution at days for sensitivity to (left) and (right). The sensitivity to is based primarily on currently affected regions, reflecting the state of epidemic progression. The sensitivity to , corresponds primarily to highly populated regions. Even though the number of exposed and infected patients is low in certain heavily populated regions (particularly the area around Milan, in the west), the high susceptible sensitivity shown here indicates the region’s vulnerability to the pandemic (which does eventually occur)

Finally, we find the definition of given by Eq. (69) to be less useful as a predictor of disease spread than in the 1D simulation. This is likely due to the increased role of diffusion. We observe increase in disease exposure and infected individuals in areas where both locally and globally, particularly around Milan (as shown in Fig. 15). This suggests the need to revise the definition of in Eq. (69) for the PDE version of the model to account for the influence of diffusion.

Fig. 15

Comparison between value and infected population. Even though globally, we still observe growth in some regions, suggesting that the definition (69) of does always not hold in the presence of diffusion

Conclusions

In this work, we introduced a new notational framework for understanding reaction–diffusion compartmental models by interpreting them as balance equations similar to those found in continuum mechanics. We first used this system to derive and explain a simple two-compartment Lotka–Volterra model as a simple example. We then examined a more complex compartmental system: the model of COVID-19 spatiotemporal contagion dynamics introduced in [22]. We showed that this model may be regarded as a sort of conservation law, further justifying the continuum-mechanics type interpretation.

We proceeded to formally derive the model’s sensitivity to diffusion, describe its growth and decay, and establish its stability in the norm. We then looked at an ODE version of the model, using it to derive a basic reproduction number as well as analyzing its spectrum. Additionally, we performed a series of numerical simulations, showcasing the role that numerical methods, diffusion, and play in the behavior of the system. We found that implicit models are effective in describing the temporal dynamics of the system, and second-order in-time methods in particular. We also found that the ODE-based is not consistently reliable as applied to the PDE model, as it worked well for the 1D simulations but did not for the corresponding 2D simulations.

For future work on the COVID-19 model, we would like to extend the diffusion to model the effects of geographic features like roads, rivers, and mountains. We would also like to examine the effectiveness of the model over larger geometries and longer time intervals against measured data. To render the model more effective to decision-makers, incorporating an age-structured population is important for accurately evaluating aspects such as hospitalizations and mortality. The model may also be extended to account for the effects of vaccination on adults, by introducing movement between the susceptible and recovered population [1]. More generally, we would like to apply the continuum mechanics framework shown here to a larger class of compartmental models. In particular, in the field of mathematical epidemiology alone, there are many variants of the SIR-models shown here. For instance, the framework established in the present work may be used for susceptible-infected-susceptible models (such as those used for the common cold), or Maternal-Susceptible-Exposed-Infected-Recovered models in which immunity is inherited from the mother [1, 5, 20].