Delay differential equations for the spatially resolved simulation of epidemics with specific application to COVID-19.

In the wake of the 2020 COVID-19 epidemic, much work has been performed on the development of mathematical models for the simulation of the epidemic and of disease models generally. Most works follow the susceptible-infected-removed (SIR) compartmental framework, modeling the epidemic with a system of ordinary differential equations. Alternative formulations using a partial differential equation (PDE) to incorporate both spatial and temporal resolution have also been introduced, with their numerical results showing potentially powerful descriptive and predictive capacity. In the present work, we introduce a new variation to such models by using delay differential equations (DDEs). The dynamics of many infectious diseases, including COVID-19, exhibit delays due to incubation periods and related phenomena. Accordingly, DDE models allow for a natural representation of the problem dynamics, in addition to offering advantages in terms of computational time and modeling, as they eliminate the need for additional, difficult-to-estimate, compartments (such as exposed individuals) to incorporate time delays. In the present work, we introduce a DDE epidemic model in both an ordinary and partial differential equation framework. We present a series of mathematical results assessing the stability of the formulation. We then perform several numerical experiments, validating both the mathematical results and establishing model's ability to reproduce measured data on realistic problems.

INTRODUCTION

The worldwide outbreak of COVID‐19 in 2020 has caused unprecedented disruption, leading to massive damage in terms of both economic cost and human lives. Much of the economic damage in particular has been due to government efforts designed to retard the spread of the disease, while undoubtedly effective, the cost of such measures is enormous. In recent months, much research has focused on the mathematical modeling of the epidemic, and of epidemics generally, in the hope that such models may ultimately prove useful to decision‐makers and help to inform more targeted, less‐disruptive interventions.

Many modeling approaches have been proposed, with some combining differential equations and empirical approaches in order to evaluate the effectiveness of various social‐distancing measures. , , , , In order to incorporate spatial variation across different regions, many of these models discretize various regions along geopolitical (or similar) lines, using a network structure to represent movement between the populations in different areas. , ,

In contrast to these approaches, in previous studies, , , , , , , , , the authors instead modeled the spatial diffusion of the disease using partial differential equation (PDE) models. The implemented models followed the compartmental framework but also incorporated nonlinear heterogeneous diffusion terms, giving a reaction‐diffusion system of equations. While the computational cost of such an approach is much higher than an ODE model, different numerical experiments showed that this approach carries several advantages. Notably, the timing of different dynamics across different areas can be resolved in a continuous manner, offering a richer description of the spatiotemporal evolution. Indeed, Boscheri et al and Bertaglia and Pareschi also included hyperbolic commuting terms, to model nonlocal movement of persons in addition to the local dynamics governed by the diffusion.

In the following, we propose an alternative formulation of the model introduced in Viguerie et al and analyzed further and extended in previous studies. , , , In particular, we seek to model many of the dynamics, and in particular the incubation period, with a delay differential equation model. Ordinary delay differential equation models have been extensively used for the study of epidemics, as well other types of biological models, such as predator‐prey equations. , , , , , Delayed models using partial differential equation (PDE) to study epidemics have been discussed and analyzed in previous works. , , , , , Many of these works are restricted to mathematical analysis, though simple numerical tests were also carried out in Pei et al and Zhao et al. Moreover, in Buonomo et al, information‐related delays are investigated, and in Messina, Volterra integral equations are employed in a SIS model. A delay PDE simulation of a realistic problem over a nontrivial geometry has not, to the authors' knowledge, been carried out.

There are several advantages in using a delay formulation rather than the system shown in Viguerie et al. Notably, the exposed compartment, responsible for the incubation period, may be eliminated without losing the relevant dynamics. For a PDE model in which problem size is relevant, obtaining the same dynamics with fewer compartments is desirable. However, we do not expect that the dynamics are exactly the same; we believe in fact that the delay‐equation formulation may better capture the “lag” effect incurred by the introduction of new measures (as seen during the COVID‐19 pandemic), in which there is a delay of several days between the introduction of a new public health ordinance and when its effects are fist observed. These lags may also change depending on the epidemic stage, leading to state‐dependent delays. Though we will not examine such a case here, a thorough understanding of the constant‐delay case is necessary first and is the objective of the present work. Such a formulation is also interesting from the mathematical and computational point of view, and examining such a model is worthwhile, we believe, in and of itself.

This paper is outlined as follows. We begin by recalling the model shown in Viguerie et al, along with some of its basic properties and notation. We will then proceed to introduce the delay‐equation formulation of the model for both ordinary and partial differential equation variants, highlighting important points of difference. Following this introduction, we will present several mathematical results of the delay‐differential equation models including the equilibria solutions and a stability analysis. We will then perform a series of numerical examples using the ODE and an idealized 1D problem for the PDE (inspired by Viguerie et al and Grave & Coutinho ) to qualitatively analyze the model behavior and confirm the mathematical results. We finish our numerical tests with a simulation over the Italian region of Lombardy using real data, in order to validate the model's ability to reproduce real‐life data on realistic problems, before concluding with several suggestions for future research in this area.

MODEL

The COVID‐19 PDE model presented in Viguerie et al and further analyzed and extended in previous studies , , , , reads while a similar ODE variant, neglecting diffusion, reads

The mechanism of the model is diagrammed in Figure 1 and operates in the following way: The susceptible population s is exposed to the disease by contact with exposed individuals in compartment e or infected patients in compartment i at rates β and β , respectively. After an incubation period σ, exposed individuals develop symptoms and move to the infected subgroup i. A fraction of symptomatic patients recover at a rate ϕ , moving into the recovered subgroup r. However, the remaining infected patients eventually die at a rate ϕ . The model also features asymptomatic transmission, which has been considered a key driving force in the COVID‐19 pandemic. To this end, we include a fraction of the exposed population e that directly moves to the recovered subgroup r, without ever entering in the symptomatic infected compartment i. Note that denotes the entire living population.

FIGURE 1

Flow chart describing the evolution of the various compartments and parameters in the model equations (1)–(5) [Colour figure can be viewed at wileyonlinelibrary.com]

We briefly make a few additional remarks regarding this model. The first is that this model operates on the principle of mass‐action, with the contact terms β non‐normalized and hence dependent on local population densities, reflected in their units of 1/(Time·Persons). , , , , The spatial dependence of contagion is further augmented by the addition of the Allee term A, which accounts for the tendency of COVID‐19 cases to cluster in areas where n  > > A. This term has been used extensively in other settings, with the form used above inspired directly by applications in cancer modeling. , The Allee term works to reduce transmission in areas where the population density is under a given threshold A, by bringing the population in the exposed compartment to the susceptible compartment. Consequently, in such areas, the population in compartments e and i tends to zero, eventually canceling out the transfer term. We observe that as s,  i, and e are all less than n by definition, we do not expect blowup of this term, even for very small n. We note that the Allee term does not appear in the ODE model (6)–(10). While one may technically include it, the lack of spatial variation in population density limits its usefulness from the modeling point of view, and conceptually, its inclusion only makes sense in the PDE model (1)–(5) for this reason. The diffusion terms in (1)–(5) are weighted by living population n, as the model hypothesizes that diffusion of individuals is not homogeneous, but preferential and directly proportional to the population.

This model was shown in Viguerie et al to exhibit reasonably good agreement with measured data for the region of Lombardy, Italy, with later works , , showing good agreement in other regions. Further numerical and mathematical aspects were investigated in Viguerie et al, where it was also shown that models of this type can be put in the framework of continuum mechanics, and interpreted a balance of forces. Although the model demonstrates acceptable agreement with reality and operates under sound physical assumptions, it relies extensively on unknown data. In particular, the exposed compartment (which also corresponds to the asymptomatic compartment) e is difficult, if not impossible, to quantify with accuracy. We therefore propose the following modified model, which uses a delay differential equation (DDE) formulation:

The corresponding ODE version of (12)–(15) reads

We acknowledge a slight abuse of notation as, strictly speaking, σ is the inverse of the corresponding value in system (1)–(5), (6)–(10). In general, the delay may be state‐dependent; however, for the current work, we will assume that they are constant in order to simplify our analysis and computations. Note also that the definition of n is now

As one may observe, the first major difference between the systems (1)–(5) and (12)–(15) is in the influence of the incubation period. Rather than including the exposed compartment e, the incubation period is incorporated into the system as a delay term on the infected compartment. A result of this choice is that in (12)–(15), asymptomatic individuals are no longer specifically accounted for; all infected persons are considered equally. For the specific case of COVID‐19, this may be a more reasonable assumption at this point in time, as testing protocols have improved and larger portions of asympomatic patients are now detected. ,

The second major difference between (1)–(5) and (12)–(15) is the evolution of the recovered compartment r and deceased compartment d. As formulated in (1)–(5), all members of the infected compartment i are equally likely to die or recover at the same time; it does not make any distinction on these patients based on time of infection. In contrast, the recovery and mortality rates in (12)–(14) are delay‐dependent, evolving according to the infected population at a previous point in time. This is a more realistic representation of epidemic dynamics and may be useful when considering the allocation of public health resources.

Relationship between the PDE and ODE

The presented delayed PDE model (12)–(15) and ODE model (16)–(19) are related, since they are used to describe the same phenomenon but differ due to the presence of diffusion and the Allee term A. This spatial information obviously gives the PDE model a richer descriptive capacity; however, it is also the case that the mathematical analysis and numerical simulations using the PDE require significantly more effort, both in terms of computational time and the complexity of the model implementation. For this reason, the question of when the ODE model may provide a reasonable surrogate for the PDE is an important one.

As our analysis in the following sections will demonstrate, close to the zero equilibrium; that is, when all quantities are small, the spatial (diffusive) terms, which are quadratic, are negligible and the PDE becomes an ODE with delay terms. In the case when , which is considered for example in our 1D simulations, a complete stability analysis of the equation allows to obtain sharp rigorous stability bounds which emphasize the dependence of stability of the steady state with respect to the delay. In the case when A is nonzero, such stability bounds can be interpreted in an approximate way. For reasonably small variations in local population densities, the ODE solution may still well approximate the spatially integrated PDE solution. Thus, if the spatial transients are not considered important, as may be the case in certain applications, the ODE model may be a more practical choice than the PDE model for its computational convenience. We will discuss the relationship between the two models formally in the analysis section.

Moreover, we will qualitatively compare the behavior of the solutions obtained with the two models, in order to illustrate the theoretical results in the numerical experiments sections.

ANALYSIS

In this section, we will analyze the DDE models (12)–(15), (16)–(19) mathematically. In particular, we examine the equilibrium solutions of (16)–(19) and (12)–(15) for the case and their stability properties. We then proceed to analyze the scalar linear equation associated to (16)–(19), (12)–(15) with , deriving stability conditions in terms of the physical parameters. We then examine the general case of (12)–(15) for A ≠ 0 and analyze the impact of the Allee term on the stability behavior.

Equilibria and their stability

It is straightforward to note that the only equilibrium of (16)–(19) is The linearized system is given by

By adding Equations (22), (23), (24) to 21 and making use of the variable n (as defined in (20)), instead of s, we may replace (21) with

Due to the last equation (24) of the DDE system, no contractivity arguments can be used to infer about the asymptotic stability of the solution.

Stability of the scalar linear delay equation

In order to state the stability theorem, we need to consider the scalar equation with the delay τ an arbitrary but fixed positive constant.

By the change of variables , we are led to the equation

The analysis of the characteristic equation, which for b ≠ 0 possesses infinitely many solutions , gives indications about the asymptotic behavior of the solution of (27). The general solution is then obtained as a sum of exponentials: The zeros of (28) are plotted in the left picture of Figure 2 for the case . Notice that they all lie in the left half‐plane, so that the solution will tend to zero for t → ∞ (despite a positive a).

FIGURE 2

Asymptotic stability region of Equation 27 in the (a, b)‐plane (right); roots of the characteristic equation (left) for [Colour figure can be viewed at wileyonlinelibrary.com]

We provide the stability region in the right picture of Figure 2 (see, e.g., Bellen & Zennaro ). Note that if , then the zero solution is asymptotically stable for b ∈ (− π/2, 0) and equivalently–referring to Equation 26—it is asymptotically stable if bτ ∈ (− π/2, 0). The transcendental curve bounding the stability region is expressed in parametric form as which might be expressed as a monotonically increasing function b(a). Asymptotically, the curve approaches the line ; as ϕ → 0, it tends to the point (1, − 1), and at , it crosses the point (0, − π/2).

Assume

α − μ < 0 or ;

.

Then, the zero equilibrium of (21)–(24), (25) is stable.

Consider first the DDE which has the form (26). By non‐negativity of μ, we have that if the characteristic equation has all roots in the negative complex plane so that . Moreover, it can be shown that the decay is exponential (see Diekmann et al. ).

Next, look at from which we obtain boundedness of d as a consequence of the exponential decay of i(·). Similarly, looking at we have exponential decay if μ > 0 and simply boundedness if . Finally, for what concerns the equation we conclude in a similar way, which means that if α − μ < 0 n(t) → 0 as t → ∞; otherwise, if , the stability statement holds trivially. □

Stability of the equilibrium of the PDE

The analysis of the PDE (12)–(15) when leads to the same stability Theorem 3.1. This follows immediately from the fact that, in the linearization of (12), (13), and (14), the terms and are quadratic and hence do not influence the analysis. As a consequence, the linearized system is formally analogous to the one obtained for the simpler DDE model (16)–(19).

More in detail, looking at the PDE‐based model, we have that, at the equilibrium, . Neglecting second‐order terms gives the system

Observe the following:

If , we formally reobtain the same system (21)–(24) so that Theorem 3.1 applies unchanged. In fact, in such a case, we may view the linearized system (21)–(24) as the system (29)–(32) integrated in space.

The case A ≠ 0 is more involved.

By adding Equations (30), (31), (32) to 29 and making use of the variable n instead of s, we may replace (29) with Since and s, i and r are non‐negative, we have that which allows us to rewrite (29)–(32) as

Looking at the second equation in (34), we recognize a DDE of the form A well‐known asymptotic stability condition for (35) is given by (see, e.g., Bellen & Zennaro ) which yields In the usual case where β ≥β , we obtain that the condition implies asymptotic stability (of contractive type) of the solution i independently of the delay σ.

We remark that this is a sufficient and not necessary condition for a stronger type of asymptotic stability, namely, contractivity, of the solution.

Under this condition, we have that . The analysis of the remaining equations for the variables r, d, and n is analogous to the one provided in the proof of Theorem 3.1.

Note that if we could treat δ(t) as a constant (see (33)), we would get the same sufficient conditions to asymptotic stability provided by Theorem 3.1, i.e.,

α − μ < 0 or or (ii)     ,

thus depending on the delay σ. In a regime where δ(t) ∈ [0, 1] does not exhibit big variations, we expect that such conditions continue to hold true, at least approximately.

Comments

In some cases, we observe non‐physical slightly negative values of the modeled quantities. This is due to the fact that when the rightmost roots of the characteristic equation are complex conjugate. This is easily seen observing that the equation has no real roots if .

Instead, when , a real root dominates. However, since the oscillations occur when the solution approaches the steady state in the asymptotically stable regime, we may overlook the potential misbehavior.

If μ < 0, then the stability bound can be made larger; that is, M increases as μ increases (see Figure 2).

We may understand the bound (ii) of Theorem 3.1 physically as a relationship between the removal rate, as governed by the parameters ϕ and ϕ , and the time delay σ. This bound states that, for the equations to be stable, the rate of recovery and/or mortality from the disease must occur over a time scale sufficiently longer than the time delay σ (recall that ϕ and ϕ have units 1/Time while σ has units Time).

NUMERICAL IMPLEMENTATION AND EXPERIMENTS

In this section, we will perform several numerical tests to evaluate various characteristics of the model and its numerical solution. In particular, we will perform the following experiments:

Several examples using the ODE version of the model (16)–(19). Here, we observe the impact of the delay on aspects of the physical solution, including the effects on contagion and lockdown measures. We also examine the derived stability bounds and influence of different parameters.

A one‐dimensional example using the PDE model (12)–(15) based on the simulation performed in Viguerie et al and Grave and Coutinho for different values of σ and problem parameters. These examples seek to examine the solution characteristics in both quantitative and qualitative aspects on an artificial problem which shares many characteristics with a real‐world problem but remains tractable and sufficiently simple to analyze in detail. We also seek to confirm the correspondence between simulations using the ODE and spatially integrated solutions of the PDE.

A two‐dimensional example using the COVID‐19 outbreak in Lombardy, Italy, employing the PDE model (12)–(15). This simulation is similar to the ones carried out in Viguerie et al, , which were well validated against the measured data at the time of publication. This example is designed to show the viability of the delay‐equation formulation in reproducing real‐world data, as well as its performance when compared to non‐delay models.

ODE model

In order to perform the simulations for the ODE model (16)–(19), we employ the Matlab solver DDE23. As initial conditions, we set the total population n = 1,000, and as a historic function, we choose for t ∈ [− σ, 0]. Assuming and , we end up with s(0) = n − i. The final time of the simulation is days. The parameter values are reported in Table 1. * To observe the impact of the delay, we run the simulation for different values of σ: and 20 days.

TABLE 1

Parameter values for the ODE and 1D simulations

Parameter	Units	Value
β e	Persons−1· days−1	9/40
β i	Persons−1· days−1	3/32
ϕ r	Days−1	1/32
ϕ e	Days−1	1/8
ϕ d	Days−1	3/640
μ	Days−1	0
α	Days−1	0
ν‾s∗	Persons−1· days −1	3.75·10−5
ν‾e∗	Persons−1· days −1	.75·10−3
ν‾i∗	Persons−1· days −1	.75·10−10
ν‾r∗	Persons−1· days −1	3.75·10−5

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

We note that for increasing values of the delay the number of the infections are higher, i.e., in Figure 3 (left), the infection peak for (black line) is much lower than the peak for (magenta line). Furthermore, we observe that the amplitude of the peak is larger for high values of the delay. On the other hand, if the delay is too high with respect to the parameters, we could obtain non‐physical solutions, as in Figure 3 (left), where the infections becomes negative for . In the following, we will investigate the impact of the government restrictions, i.e., the introduction of lockdowns and the stability of the model from a numerical point of view.

FIGURE 3

Total infected for ϕ  = 1/32, ϕ  = 3/640 for different delay values in the non lockdown (left) lockdown (right) case [Colour figure can be viewed at wileyonlinelibrary.com]

Effect of lockdowns

Due to the relevance of the pandemic on the dailylife routine, we seek to observe the effect of government restrictions, i.e., lockdowns, on the evolution of the compartments. To this end, we run two cases, one without lockdowns, in which the contact rates β are kept constant throughout the simulation, and one with lockdowns, in which we set at days. In fact, the aim of the lockdown is to reduce the contact rate.

In Figure 3, we show the evolution of the infected compartment in time, for the different values of σ without (left) and with (right) restrictions. It is clear that the number of infections increases with the delay even in the lockdown situation, as expected. However, the lockdown restriction reduces the number of infected people by about 20%. Moreover, we observe the same effect focusing on the deceased compartment. Indeed, looking at the deaths peak in Figure 4, the deaths peak is lower in the lockdown situation.

FIGURE 4

Total deceased for ϕ  = 1/32, ϕ  = 3/640 and for different delay values in the non lockdown (left) lockdown (right) case [Colour figure can be viewed at wileyonlinelibrary.com]

We also observe that, with larger values for the delay, the effect of the lockdown measures is less readily observed. One may particularly see this in Figure 3 on the right, where the decrease of infections as a result of the lockdown begins very quickly for and progressively more slowly for larger values of σ. This is consistent with our expectations.

Stability of the ODE model

We now investigate the numerical stability of the ODE model, and in particular, we seek to examine the validity of the bound in Theorem 3.1 and verify it numerically. Looking at the deceased compartment, Figure 4 shows that the solution is stable for , since we choose the parameters according to the bounds (i)–(ii) of the Theorem 3.1. Indeed, , which means that for all our choices of the delay.

On the other hand, modifying the parameters as and , we get an unstable solution for , as shown in Figure 5. In fact, for , the oscillations are increasing in time, while for , they are smearing out. This is also consistent with the analysis, as we may expect oscillations for larger values of σ; however, if the numerical bound is respected, these oscillations should stabilize and not affect the solution asymptotically.

FIGURE 5

Total deceased for ϕ  = 3/32, ϕ  = 1/80, for different delay values [Colour figure can be viewed at wileyonlinelibrary.com]

An interesting behavior that we also observe is that if we choose and , we obtain a periodic, non‐physical behavior of the solution for , Figure 6. As shown in the figure, we observe oscillations that neither increase nor decrease, instead demonstrating what appears to be a true periodic regime. Indeed, the with these parameters, we have This suggests that, near the limit of the stability bound, solutions exhibit a periodic behavior. Considering that the oscillations decrease for ϕ ,  ϕ sufficiently below the stability bound and increase for ϕ ,  ϕ sufficiently large, this behavior is perhaps to be expected. Whether this is a mere mathematical curiosity or perhaps indicates relevant biological information is not clear and is potentially a subject of future investigation.

FIGURE 6

Evolution of susceptible compartment(top‐left), infected compartment (top‐right), recovered compartment (bottom‐left), and cumulative death (bottom‐right) for different delay values for and [Colour figure can be viewed at wileyonlinelibrary.com]

One‐dimensional PDE model

In this example, we follow basic setup inspired by the one‐dimensional example introduced in Viguerie et al and also performed in Grave and Coutinho. We will examine the behavior of the solution under various conditions, as well as the validity of the stability bound for the partial differential equation.

Problem setup

For the initial conditions, we set and e(x ∗, 0) = e 0(x ∗) as follows:

These conditions are plotted in Figure 7. Qualitatively, they correspond to a population distribution with one major population center, one moderate population center, and one lesser population center, with an initial outbreak centered in the lesser population center.

FIGURE 7

Initial conditions for the 1D PDE model. Qualitatively, the setup represents a population distribution with one major population center and one lesser population center, with an initial outbreak centered in the lesser population center [Colour figure can be viewed at wileyonlinelibrary.com]

For the parameters, we use the values shown in Table 1. We discretize in space over the unit interval with Δx = 1/2,000 and advance in time using the BDF2 scheme with . It was demonstrated in Viguerie et al that this spatiotemporal discretization resolves all dynamics satisfactorily. We note also that our choice of time step ensures that we may use previously computed solutions for our delay terms, and there is no need for interpolation.

We run the simulation for days for and 20. We seek to examine the effect of σ on contagion, but also on the efficacy of public health interventions. To this end, for each value of σ, we run the full time interval both with and without lockdown measures. For the case with no lockdown, we use the parameters shown in Table 1 over the entire time interval. For the case with lockdowns, at , we multiply all diffusive terms ν by 1/2, simulating restricted mobility, and contact terms β by 1/4, corresponding to measures such as bar closures, mask wearing, etc. We note that the problem setup is designed to resemble the ODE simulations in the preceding subsection.

In Figures 8, 9, 10, 11, we show the total (i.e., integrated in space) susceptible, infected, recovered, and deceased compartments in time, respectively, for each value of σ and lockdown configuration. We see that longer σ is associated with higher infective peaks; however, cumulative deaths, long‐term, are similar for all σ. This is consistent with our expectations, as reducing infective peaks does not necessarily correspond to fewer cases overall.

FIGURE 8

Total susceptible for ϕ  = 1/32, ϕ  = 3/640. We have stable and monotonic behavior in this compartment across all different cases. A lag is observed, as expected, for the different values of σ [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 9

Total infected for ϕ  = 1/32, ϕ  = 3/640. A longer σ is associated with more infections, but also a more dramatic decrease in infection after lockdowns are initiated. The impact of lockdowns in the infected compartment is visible more immediately, with the lag effect being more pronounced in the deceased compartment. We also observe clearly the non‐physical behavior for in both cases, with the infected compartment becoming negative [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 10

Total recovered for ϕ  = 1/32, ϕ  = 3/640. Each case is stable, but we observe nonphysical behavior as the deceased compartment for , with both cases demonstrating noticeable non‐monotonicity. The smaller amount of recovered individuals in the lockdown cases is explained by the reduced contagion overall [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 11

Total deceased for ϕ  = 1/32, ϕ  = 3/640. We see a longer σ leads to more fatalities, as well as delaying the efficacy of public health measures. However, while the effects of intervention are delayed, we note that, once visible, their impact occurs more suddenly. We observe non‐physical behavior for , showing that the model may exhibit non‐physical behaviors for larger σ [Colour figure can be viewed at wileyonlinelibrary.com]

More interesting is the observed effect of σ on lockdown efficacy. Referring to Figure 11, the results show that longer σ leads to a noticeable lag in the delay compartment; indeed, one begins to see a decrease in mortality as a result of the lockdown at a later date for larger σ. In the plot of the infected compartment in Figure 9, the impact of lockdowns appears immediate. This is indeed what we expect, as the definition of this compartment includes pre‐symptomatic patients; thus, while the effects are immediate by this definition, other indicators, such as the aforementioned deceased compartment, will lag in proportion to σ.

The delay σ appears to not only affect the time at which the effect of lockdowns appear, but also the sharpness of such effects. While for larger σ, the deceased compartment begins to decrease later as compared to smaller values, once this decrease begins, it occurs much more rapidly. This is particularly visible in the infected compartment shown in Figure 9; around , the total number of infections is indeed larger for smaller values of σ.

An important dynamic that we notice for large σ is the emergence of non‐physical behavior, similar to those observed in the ODE case. This is particularly apparent for the case of , where we observe the infected compartment becoming negative, and decreases in the recovered and deceased compartments, which should be monotonic. While this behavior is non‐physical, it is not mathematically inconsistent with the model behavior and does not represent instability as such. Indeed, to guarantee positivity, more stringent conditions on the relationship between ϕ ,  ϕ , and σ are likely required.

Stability of the PDE model

We now examine the validity of the bound (ii) in Theorem 3.1, derived for the ODE model, for the PDE model. As mentioned in the analysis section, we expect the results to hold identically in this case. In the preceding, we observe that for the parameters in Table 1, for σ = 5, 10, 15, 20, thus satisfying the condition for all considered σ. While we see non‐physical behaviors for both cases of , the solution does indeed remain stable. Such behaviors appear to be independent of both time‐step size and time‐integration scheme and represent the model behavior. Thus, to guarantee physical behavior, stability is necessary but not sufficient.

For and , again satisfying the stability condition for all σ, we again see further numerical validation of the system stability. As shown in Figure 12, all values of σ are stable. and both show physical behavior, avoiding oscillations and large decreases in the deceased compartment. While is stable, the behavior is clearly nonphysical, and we see oscillations and decreases in the deceased compartment.

FIGURE 12

Total deceased for ϕ  = 3/56, ϕ  = 3/320. We see that with this choice of ϕ, the equations become noticeably nonphysical for , with large decreases in the deceased compartment. However, we nonetheless observe stability for all σ [Colour figure can be viewed at wileyonlinelibrary.com]

For and , for but not for . The behavior of is both physical and stable, and is stable but nonphysical, as shown in Figure 13. violates the stability bound slightly, and we observe unstable behavior, with large increasing oscillations. This establishes not only the validity of the stability bound, but also its strictness. These results confirm the analysis, showing that the stability bound (ii) in Theorem 3.1, holds for the PDE model. However, as the results also show clearly, stability does not imply a physical solution, as one still may observe negative values in the infected compartment and non‐monotonic behavior in the recovered and deceased compartments. For sufficiently small values of σ in comparison to ϕ ,  ϕ , we observe numerically that one may expect physical behavior in the solution. A positivity condition, a much stronger condition than stability, is required to guarantee physical solution behavior; this likely depends on the initial data and is an interesting direction for future work.

FIGURE 13

Total deceased for ϕ  = 3/32, ϕ  = 1/80. We see that with this choice of ϕ, the equations are stable and physical for , nonphysical but stable for , and unstable for , as predicted by the stability condition [Colour figure can be viewed at wileyonlinelibrary.com]

Relationship with ODE solutions

In this case, we observe that the ODE provides a good approximation to the space‐integrated PDE in 1D. This is expected from our analysis, as in this case, we have taken and the variation in population density, while present, is not extreme. As mentioned previously, ODE models are obviously much less demanding than PDE models from the point of view of both implementation and simulation time as well as from analytical/stability analysis point of view. For this reason, if possible, their use may be preferred over PDEs in certain contexts. These simulations confirm that, near equilibria, for small values of A, and relatively small variation in population density, the ODE may provide a surrogate for the PDE, if spatial differences are not considered important and one wishes to instead consider the population independently of space. On the other hand, the PDE model provides more reliable forecast when these factors are important, due to the presence of the spatial information.

Lombardy simulation

Our final numerical test models the outbreak of COVID‐19 in the region of Lombardy, Italy. This test case was first used to validate a SEIRD model in Viguerie et al, in which the simulation results showed good agreement with the measured data. It was then examined in further detail in Viguerie et al, where other aspects of the problem, including the model sensitivity to diffusion, were investigated. We will not discuss such aspects of this model here, as the focus of the present work is on the delay differential equation model. Hence, for a more detailed discussion of the aforementioned results, we kindly refer the reader to Viguerie et al. ,

For the spatial discretization, we utilize an unstructured triangular mesh with 41,625 elements. For the temporal discretization, we use the Backward‐Euler method with a time step days. We solve the nonlinear problem at each time step using a Picard‐style linearization, and the corresponding linear systems are solved using the GMRES algorithm with a Jacobi‐style preconditioner. At all boundaries, we assign no‐flux boundary conditions, corresponding to total isolation. The parameter values are reported in Table 2. We note that these differ from those shown in Viguerie et al.; , this is primarily due to differences resulting from the fact that, in the current delay model, asymptomatic individuals are now considered within the “infected” compartment.

TABLE 2

Parameter values for the 2D Lombardy simulations

Parameter	Units	Feb 27 to Mar.9	Mar 9–22	Mar 22–28	Mar 28 to May 3	May 3–
σ	Days	8	8	8	8	8
β e	Persons−1·days−1	3.75·10−4	3.11·10−5	2.0625·10−5	1.5·10−5	2.75·10−5
β i	Persons−1·days−1	3.75·10−4	3.11·10−5	2.0625·10−5	1.5·10−5	2.75·10−5
ϕ r	Days−1	3/64	3/64	3/64	3/64	3/64
ϕ d	Days−1	3/320	3/320	3/320	3/320	3/320
ν‾_s	km2· 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	km2· Persons−1·days−1	2.175·10−2	1. · 10−2	0.45·10−2	0.325·10−2	1.0625·10−2
ν‾_r	km2· Persons−1·days−1	4.35·10−2	1.98·10−2	0.9·10−2	0.75·10−2	2.175·10−2
Ā	Persons	1.0·103	1.0·103	1.0·103	1.0·103	1.0·103

Note: 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 Figure 14, we show how the results compare to two formulations of the problem shown in Viguerie et al.: the “baseline” case, which features the same diffusive parameters used here, and the case with doubled diffusion. Note that we focus here on the deceased compartment, rather than the infectious compartment, due to the differing definitions “infected” between the two models. We see similar qualitative behavior to the baseline case, with an R 2 correlation coefficient between the of 99.8%. Over the first  40 days, the behavior is nearly identical, with the results beginning to differ somewhat further in time. The agreement with the baseline case, rather than the high‐diffusion case from Viguerie et al, suggests that the delay model does not interfere with the diffusive behavior generally.

FIGURE 14

Deceased individuals for the Lombardy test case. We show the delay PDE results together with two other non‐delay model formulations, validated against measured data. We see that the results exhibit similar qualitative behavior, establishing the potential of the delay PDE model to produce realistic simulation results in time and space [Colour figure can be viewed at wileyonlinelibrary.com]

In Figure 15, we observe the evolution of the epidemic in space. Starting from the top‐left and moving clockwise, we show the density of infected individuals on days 1, 5, 10, and 20. What begins as a small cluster of cases in the south of the region moves northward, into the region's large cities (Milan, Bergamo, and Brescia). While the initially affected regions improve rapidly, in the large cities, the epidemic continues to grow. This is consistent with the observed data, and with the simulations using non‐delay models shown in Viguerie et al.

FIGURE 15

Evolution of the epidemic in Lombardy using the Delay PDE. Clockwise, from top left: Days 1, 5, 10, and 20. The outbreak begins with a clusters of cases in the south of the region, before moving north into the region's large cities [Colour figure can be viewed at wileyonlinelibrary.com]

In order to assess the stability bound, we performed a second numerical simulation. We note that, strictly speaking, as A is nonzero in this case, the assumptions of the stability bound do not strictly hold. However, in other simulations (not shown), the qualitative behavior of the model with was similar to that shown here, and so, we may expect the bound to hold heuristically. To better observe the stability behavior, we set from February 27 to March 4, and for the remainder of the simulation. We then let † . Hence, In accordance with the theory, we do not expect stability for this case. Indeed, this is what we observe. In Figure 16, we plot the total number of active infections for the stable (previously discussed) case, as well as the unstable case. While in the stable case, we observe the expected behavior, in which infections remain positive and grow or decrease in response to the pandemic‐arresting measures in place, we instead observe large oscillations between positive and negative values for the unstable case. We show this behavior in time and space in Figure 17, where one sees the oscillations concentrated in the heavily‐affected outbreak zones within the region. The frequency of these oscillations appears to be nonuniform, and dependent on infection concentration. Whether this provides any important information is, unclear, though it is an interesting observation and perhaps worthy of some investigation.

FIGURE 16

Active infections for the Lombardy simulations in both the stable and unstable regimes. In the stable regime, we see the expected behavior; infections remain positive and change in response to pandemic‐arresting measures. In the unstable regime, we instead see spurious nonphysical oscillations between positive and negative [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 17

The simulation of Lombardy in the unstable regime, clockwise from left: days 20, 40, 60, 80. We see interesting oscillatory behavior between positive and negative throughout the region and the primary infectious zones [Colour figure can be viewed at wileyonlinelibrary.com]

CONCLUSIONS

We have presented a new formulation for epidemic models utilizing delay differential equations in both an ODE and PDE formulation. We have established stability results for the ODE formulation, which were then confirmed with numerical experiments. For the PDE model, we observed interesting dynamics regarding the relationship between the delay time and lockdowns, as well as contagion in general. We further showed, with numerical evidence, that the stability bounds established for the ODE also hold for the PDE and that, in some situations, the ODE may a reasonable surrogate for the PDE. We then concluded with a simulation on a realistic problem, in which we showed that the delay PDE can reproduce reality at a reasonable level, by obtaining results similar to non‐delay models shown in other work. We also analyzed the spatiotemporal behavior of the unstable regime, finding it produced large oscillatory behavior between positive and negative values within the heavily impacted regions.

There are many worthwhile directions for future work on this model and the area generally. One could consider for example other compartments, as vaccinated, and/or structure the population on different age classes, i.e., young, adult, and elderly people, or in a fully continuous manner as done in Murray. , While we provided theoretical and numerical evidence of stability, we also showed for both the PDE and ODE that stability does not necessarily guarantee physical solution behavior. A stronger result, such as a positivity condition, is needed to provide such a result. Lastly, we have restricted ourselves to the constant delay case as a first step, but the most general and realistic models of this type should incorporate state‐dependent delays. Although such models have many theoretical and numerical difficulties, many epidemics and other related phenomena exhibit this type of behavior. , , , , Using the proposed model to solve inverse problems, for example, as the forward model in a procedure to identify a state‐dependent delay parameter, is also a planned subject of future work.

CONFLICT OF INTEREST

The authors have no conflict of interest to disclose.