Partial differential equation (PDE) model
In this paper, we adopt a compartmental modeling approach that may be traced back to the seminal work of Kermack and Mckendrick27 in mathematical epidemiology. Essentially, this approach partitions a population subject to an epidemic outbreak into a distinct number of disjoint compartments; in the simplest models, these are the Susceptible (S), Infected (I), and Recovered (R) compartments. Infection spread is governed by a transition rule representing the mapping of individuals from the Susceptible to the Infected compartment (i.e. the process of susceptible individuals becoming infected). Recovery is likewise represented by a transition rule that maps individuals from the Infected to the Recovered compartment. Moreover, the overall dynamics of epidemic evolution are represented by a system of coupled differential equations for the population densities S, I, and R, with the coupling terms representing the corresponding compartmental transitions. As an epidemic evolves—typically through multiple outbreaks—I correspondingly records multiple local maxima to ultimately approach and stay bounded around a stable equilibrium value. Attainment of this endemic equilibrium marks the end of the epidemic outbreak. We note that the I need not necessarily vanish at this equilibrium but attain a value below the threshold required to trigger further outbreaks.
In the compartmental models discussed above, there exists a fundamental distinction between the two types of differential equations that characterize a model—ODE and PDE. ODE-based SIR models consider the compartmental density functions (S, I, R) exclusively as functions of time; consequently, such models only account for temporal epidemic dynamics. In other words, ODE models only track changes in the number of individuals in each compartment over time and, therefore, are unable to address the spatial dynamics of epidemic spread28. On the other hand, PDE models account for spatial variations that emerge from the interactions in a population and thereby capture the heterogeneity of disease spread across different locations. Specifically, the coupled PDEs in such a model describe how the distribution of individuals in each compartment changes not only over time but also across space29, thereby describing spatiotemporal disease transmission.
A typical characteristic of infection spread is the latency period when an individual is infected but not infectious yet30. Accounting for this feature motivates the introduction of a Latent (L) compartment into epidemic models. Yet another subtle but significant aspect related to latency and spatiotemporal spread is the spatial displacement of infectious individuals in whom the infection is latent. This aspect may be accounted for by the introduction of a convolution term in PDE models (see, for instance, Li et al.13).
Finally, it is remarkable that diffusive PDE models-viewed through the lens of statistical mechanics—can naturally provide a probabilistic interpretation of epidemic spread dynamics. The underpinnings of this interpretation lie in the fact that the probability density functions that describe diffusive stochastic processes satisfy PDEs similar to those that arise in PDE epidemic models. Given that uncertainties—due to random human behavior and random transmission characteristics of the pathogen that triggers an epidemic—play a critical role in dictating epidemic spread, the probabilistic interpretation further endows PDE models with immense power to analyze stochastic epidemic spread.
In light of the above discussion, we now consider the compartmental epidemic model described by the following coupled system of PDEs:
$$\begin{aligned} & \frac{\partial S }{\partial t} = \lambda + \eta _{S}\nabla ^{a}S – \theta S – \phi IS \end{aligned}$$
(1)
$$\begin{aligned} & \frac{\partial L}{\partial t} = \eta _{L}\nabla ^{a}L +\phi IS – \epsilon \phi \times \left( \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } I(x,y,t – \tau )S(x,y,t-\tau )f_{a}(x,y)dxdy\right) \end{aligned}$$
(2)
$$\begin{aligned} & \frac{\partial I}{\partial t} = \eta _I \nabla ^{a}I- \delta I + \epsilon \phi \times \left( \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } I(x,y,t – \tau )S(x,y,t- \tau )f_{a}(x,y)dxdy\right) \end{aligned}$$
(3)
$$\begin{aligned} & \frac{\partial R }{\partial t } = \eta _R \nabla ^{a}R- \theta R+ \omega I \end{aligned}$$
(4)
\(\lambda\) is the birth rate in the spatial domain. \(N_{x,y}\) is the total population in cell (x, y), remains constant throughout the specified time period of simulation. \(\eta _{S}\) is the diffusion coefficient representing the intensity of random motion of the Susceptible population. \(\eta _{L}\) is the diffusion coefficient representing the intensity of the random motion of the Latent population. \(\eta _{I}\) is the diffusion coefficient representing the intensity of random motion of the Infected population. \(\eta _{R}\) is the diffusion coefficient representing the intensity of random motion of the Recovered population. \(\nabla ^{a}\) is the Laplacian operator representing the Brownian motion of individuals in a system. \(\theta\) is the mortality rate due to natural deaths. \(\phi\) is the infection rate. \(\delta\) is the rate of individuals exiting the infected compartment due to deaths or recovery. \(\omega\) is the recovery rate of individuals. \(\epsilon\) is the fraction of the infected population that survives the latency period and enters the infected category. \(f_{a}(x,y)\) is the Gaussian kernel defining the extent of the mobility of the latent population. \(\tau\) is the latency period of the model. S, L, I & R is the population densities of Susceptible, Latent, Infected, and Recovered individuals at a specific location and time.
Focusing now on the details of the model equations—which are adapted from Li et al.13—firstly, we reiterate that the population is categorized into disjoint subsets representing the population densities of Susceptible (S), Latent (L), Infectious (I) and Recovered (R) individuals. Note that the population densities are simultaneously functions of the (Cartesian) spatial variables (x, y) as well as time t, with the spatiotemporal dynamics of an epidemic dictated by the system of PDEs as individuals transition through the Susceptible Latent, Infected, and Recovered compartments. Secondly, we note that the diffusion terms in the equations (the Laplacian operator \(\nabla ^a\)) represent the random spatial motion of individuals in the respective compartments. This feature is inherited from the fact that invoking the theory of Brownian motion, under appropriate normalization conditions, the population densities in the model may be interpreted as probability density functions describing diffusive stochastic processes. Thirdly, in our previous work, the efficacy of this PDE model to predict real-world epidemic spread was validated using empirical COVID-19 data from both a larger spatial domain (the State of Ohio, USA) as well as a smaller one (Hamilton County, Ohio, USA)31. Indeed, the present effort is motivated by and grounded upon the aforementioned validation analysis.
PDE model: numerical framework
Our objective of computing reproduction numbers from the PDE model involves numerically solving the PDEs. To this end, guided by our previous work31, we now discretize the system of PDEs using Euler’s forward algorithm as follows.
$$\begin{aligned} & S^{T+1}_{x,y} = N\lambda + \eta _{S}(S^{T}_{x,y + 1} + S^T_{x,y -1} + S^{T}_{x+1,y} + S^{T}_{x-1,y}) + (1-4\eta _{S} – \theta – \frac{\phi }{N}I^{T}_{x,y})S^{T}_{x,y} \end{aligned}$$
(5)
$$\begin{aligned} & L^{T+1}_{x,y} = \eta _{L}(L^{T}_{x,y + 1} + L^T_{x,y -1} + L^{T}_{x+1,y} + L^{T}_{x-1,y}) + (1-4\eta _{L})L^{T}_{x,y} +\frac{\phi }{N}I^{T}_{x,y}S^{T}_{x,y} – \epsilon \frac{\phi }{N}\sigma ^{T}_{x,y} \end{aligned}$$
(6)
$$\begin{aligned} & I^{T+1}_{x,y} = \eta _{I}(I^{T}_{x,y + 1} + I^T_{x,y -1} + I^{T}_{x+1,y} + I^{T}_{x-1,y}) + (1-4\eta _{I} – \delta )I^{T}_{x,y} + \epsilon \frac{\phi }{N}\sigma ^{T}_{x,y} \end{aligned}$$
(7)
$$\begin{aligned} & R^{T+1}_{x,y} = \eta _{R}(R^{T}_{x,y + 1} + R^T_{x,y -1} + R^{T}_{x+1,y} + R^{T}_{x-1,y}) + (1-4\eta _{R} – \theta )R^{T}_{x,y} + \omega I^{T}_{x,y} \end{aligned}$$
(8)
$$\begin{aligned} & \sigma ^{T}_{x,y} = \sum _{i=1}^{\mathscr {N}}\sum _{j=1}^{\mathscr {N}}I^{T – \tau }_{x,y}S^{T-\tau }_{x,y}f_{a} l ^{2} \end{aligned}$$
(9)
$$\begin{aligned} & f_{a}(x,y) = \frac{1}{\sqrt{4\pi \alpha }}e^{-\frac{(i-x)^{2}+ (j-y)^{2}}{4\alpha }} \end{aligned}$$
(10)
Here \(X^{T+1} _{x,y}\) denotes the value of variable X (\({X \in \{S,L,I,R\}}\)), on day \((T + 1)\) in the cell with discretized spatial coordinates (x, y) in the two-dimensional domain under consideration, as a function of the values of the corresponding variables on the previous day T. In the discretization, \(\mathscr {N}\) is the number of cells in each dimension of the domain, and \(\sigma ^{T}\) denotes the latent population in a cell on day (T). Two points regarding how latency is taken into account are in order. Firstly, \(\tau\) represents the latency period; therefore, \(\sigma ^{T}\) is computed for the day T based on \(I^{T-\tau }\) and \(S^{T-\tau }\)—the values of I and S from \((T-\tau )\) days ago. Secondly, we consider the realistic scenario of spatial migration of the infected population during the latency period contributing to the spread. This is represented by the convolution kernel \(f_a\) on the RHS of Eq. (9), the explicit form of which is given by the Gaussian distribution in Eq. (10). We note that the Gaussian distribution captures the diminishing probability of people moving to a farther location compared to a closer one13.
We consider no-flux boundary conditions in the numerical solutions of the PDEs. In other words, we assume no traffic of individuals in either direction at the boundary of the spatial domain. Moreover, this no-flux condition is enforced in the numerical solution by having“reflecting” boundaries for the grid, i.e. attempted transgression of the boundary will result in individuals bouncing back into the cells that define the boundary. This assumption allows us to maintain the law of conservation with respect to the population of the system.
Basic reproduction number: analytical formulation
The basic reproduction number \(\mathscr {R}_0\) is defined as the average expected number of secondary infections that occur in a fully susceptible population due to contact with an infectious individual. A disease spreads if the reproduction number is greater than 1; contrapositively, the spread ceases if the reproduction number is less than 1. The infection becomes an epidemic outbreak if \(\mathscr {R}_0\) stays persistently greater than 110. On the other hand, if \(\mathscr {R}_0\) stays consistently less than 1 for a considerable amount of time, it indicates a waning phase of an epidemic32.
A key idea in this work is to numerically compute the reproduction number \(\mathscr {R}_0\) and hence \(\mathscr {R}_e\) for a spatiotemporal domain from a discretized version of the PDE model. To establish the mathematical basis for this computation, we first note the fundamental result from ODE-based epidemic models that \(\mathscr {R}_0\) is the principal or largest eigenvalue (i.e., the spectral radius) of the next-generation matrix derived from an ODE model4,6,32,33. Whilst this result resists straightforward generalization to the PDE case, approaches such as employing a variational formula have been reported for basic Susceptible-Infected-Susceptible (S, I, S) models with diffusion (see, for instance, Allen et al.34). In particular, Wang et. al.19,35 established \(\mathscr {R}_0\) for reaction-diffusion PDE systems as the spectral radius of the next-infection operator, using the theory of principal eigenvalues. We find that this approach can be fruitfully adapted for our purposes and now provide an outline of the development, referring to the previous literature19,24,35 for details.
At the outset, we note that out of the four compartments in the model that we consider, the Latent and Infected compartments form a subsystem of infected states, whereas the Susceptible and Recovered compartments are the uninfected states. Also, note that the initial infection-free state corresponds to \(L = I = R = 0\).
Consider a compartmental model comprising n compartments, where n is the disjoint union of an arbitrary number of infected and susceptible compartments.
Denoting the population density in the \(i^{th}\) compartment by \(u_i(t,\textbf{x})\), let us collect the populations in all compartments in a column vector \(u(t,\textbf{x}) =(u_1(t,\textbf{x}), \dots , u_n(t,\textbf{x}))^T\), where each \(u_i\) is non-negative, \(\textbf{x}=[x, y]\) is the position vector in 2-dimensional space \(\textbf{R}^2\) represented using Cartesian coordinates, and T is the standard transpose operator. Next, the total number of compartments is categorized into two distinct (disjoint) groups: m infected compartments, denoted by \(i=1,\ldots ,m,\), and the remaining uninfected compartments, identified by \(i=m+1,\ldots ,n\). A generic reaction-diffusion system of PDEs (along with the necessary boundary conditions) representing the aforementioned compartmental model can now be written as:
$$\begin{aligned} \left\{ \begin{array}{ll} \frac{\partial u_i}{\partial t} = \nabla \cdot (d_i(\textbf{x}) \nabla u_i) +\mathscr {F}_i(\textbf{x},u) – \mathscr {V}_i(\textbf{x},u), & \\ \\ \hspace{35mm}1 \le i \le n, \, t> 0, \, \textbf{x} \in \Omega , & \\ \\ \frac{\partial u_i}{\partial \nu } = 0 \quad \forall 1 \le i \le n ; d_i> 0, \, t > 0, \, \textbf{x} \in \partial \Omega .& \end{array} \right. \end{aligned}$$
(11)
In the above equation, \(u_i\) is the density of individuals in the \(i^{th}\) compartment, while \(d_i(\textbf{x})\) is the corresponding diffusion coefficient (representing the random movement of the population \(u_i\)). The terms \(\mathscr {F}_i\) and \(\mathscr {V}_i\) are the reaction terms in the \(i^{th}\) compartment (accounting for transitions between compartments induced by the disease spread). Specifically, \(\mathscr {F}_i(\textbf{x},u)\) characterizes the rate at which the newly infected population is introduced into a given compartment i. Also,
$$\begin{aligned} & \mathscr {V}(\textbf{x},u) = (\mathscr {V}^-(\textbf{x},u)-\mathscr {V}^+(\textbf{x},u))= (\mathscr {V}_1^-(\textbf{x},u),\ldots ,\mathscr {V}_n^-(\textbf{x},u))^T – (\mathscr {V}_1^+(\textbf{x},u),\ldots \mathscr {V}_n^+(\textbf{x},u))^T \\ & \quad = (\mathscr {V}_1(\textbf{x},u),\ldots ,\mathscr {V}_n(\textbf{x},u))^T \end{aligned}$$
where \(\mathscr {V}^+\) denotes the transfer rate of individuals into this compartment through all means, and \(\mathscr {V}^-\) is the transfer rate of individuals out of the compartments. We note that the PDEs are defined in the spatial domain \(\Omega\) (a subset of \(\mathbb {R}^2\)), which is assumed to have a smooth (i.e. continuously differentiable) boundary \(\partial \Omega\). The symbol \(\nu\) represents the outward-facing unit normal vector along the boundary \(\partial \Omega\). The second of the equations 11 is the no-flux boundary condition and reflects the physical assumption that the boundary is impermeable to individual movement in both the inward and outward directions.
It is useful here to write the system of Eq. (11) in an equivalent form by expanding the first term of the RHS of the PDE using the product rule for differentiation as:
$$\begin{aligned} \left\{ \begin{array}{ll} \frac{\partial u_i}{\partial t} = d_i(\textbf{x}) \nabla ^2 u_i – c_i(\textbf{x}) \nabla u_i +\mathscr {F}_i(\textbf{x},u) – \mathscr {V}_i(\textbf{x},u), \le i \le n, \, t> 0, \, \textbf{x} \in \Omega , & \\ \\ \frac{\partial u_i}{\partial \nu } = 0 \quad \forall 1 \le i \le n ; d_i> 0, \, t > 0, \, \textbf{x} \in \partial \Omega . & \end{array} \right. \end{aligned}$$
(12)
where,
$$\begin{aligned} c_i(\textbf{x})= – \frac{d}{d\textbf{x}}d_i(\textbf{x}) \end{aligned}$$
(13)
Before proceeding to define the basic reproduction number \(\mathscr {R}_0\) corresponding to the PDE Eq. (12), we now briefly discuss the key ideas behind the Next Generation Operator defined in terms of the key operators \(\mathscr {F}_i\) and \(\mathscr {V}_i\). Firstly we note that these operators were first defined in the context of ODE models as matrices32. Accordingly, the Next Generation Matrix (NGM), is obtained in terms of the F and V matrices as \(FV^{-1}\). The F matrix is called the “new infection matrix” or the “transmission matrix”. It represents the rate at which new infections are generated by the currently infected individuals and is an \(m\times m\) matrix (considering only the infected group) with its \((p,q)^{th}\) entry given by:
$$\begin{aligned} F_{pq} = \frac{\partial \mathscr {F}_p}{\partial u_q} (u^0). \end{aligned}$$
(14)
Here \(u^0\) represents the disease-free equilibrium point. These matrices are evaluated at the Disease Free Equilibrium (DFE), where no infections are present, serving as a crucial starting point for determining \(\mathscr {R}_0\), which quantifies the potential for disease spread in a fully susceptible population. Moreover, each element of the matrix \(F_{pq}\) represents the expected number of new infections in the compartment p produced by an infected individual in compartment q over a single infectious period in a fully susceptible population. On the other hand, the V matrix, which is known as the “transition matrix” or “transfer matrix”, characterizes the movement of individuals between different compartments due to recovery, progression of the disease, or death. It is also an \(m\times m\) matrix (once again considering only the infected group) with its \((p,q)^{th}\) entry given by:
$$\begin{aligned} V_{pq} = \frac{\partial \mathscr {V}_p}{\partial u_q} (u^0). \end{aligned}$$
(15)
The basic reproduction number \(\mathscr {R}_0\) is then established as the spectral radius of the NGM. The underlying rationale is the fact that the largest absolute eigenvalue—which is the spectral radius—of the NGM matrix represents the peak average number of new infections produced by a single infected individual. This takes into account the rate of transfer among various stages of the disease within a single infection cycle. Essentially, this value reflects the most significant rate of disease transmission in a population, considering the detailed progression and the contributing interactions.
Refocusing our attention on the PDE, and referring to the discussion in19,24 for more details, consider the following vector-valued reaction–diffusion equation. We note that this is identical to the previous equations in 12, except that here only all the infected compartments are considered. Moreover, the corresponding population densities are collected in a column vector. The vector-valued PDE is:
$$\begin{aligned} \left\{ \begin{array}{ll} \frac{\partial u_\mathbb {I}}{\partial t} = d_\mathbb {I}(\textbf{x}) \nabla ^2 u_\mathbb {I} – c_\mathbb {I}(\textbf{x}) \nabla u_\mathbb {I}- V(\textbf{x})u_\mathbb {I}, \quad t> 0, \textbf{x} \in \Omega ,& \\ \\ \frac{\partial u_\mathbb {I}}{\partial \nu } = 0, \quad d_i> 0, \, \forall 1 \le i \le m, \, t > 0, \, \textbf{x} \in \partial \Omega , & \end{array} \right. \end{aligned}$$
(16)
Where \(u_\mathbb {I}\)= \((u_1,…,u_m)^T\) and \(d_\mathbb {I}(\textbf{x})= diag(d_1(\textbf{x}),\ldots , d_m(\textbf{x}))\) and \(c_\mathbb {I}(\textbf{x}) = diag(c_1(\textbf{x}),\ldots , c_m(\textbf{x}))\).
Now, let \(\mathscr {T}(t)\) be the operator that temporally propagates an initial spatial distribution that satisfies the PDE Eqn. 16. If we take the initial infection distribution (arrangement of infections at the beginning of the observation period) to be \(\zeta (\textbf{x})\)= \((u_1(t=0,\textbf{x}),\dots ,u_m(t=0,\textbf{x}))^T\), then, \(\mathscr {T}(t)(\zeta (\textbf{x}))\) represents the distribution of the infections for any \(t > 0\). Therefore, the spread of new infections at any time \(t > 0\) can be given as \(F\mathscr {T}(t)(\zeta (\textbf{x}))\). The distribution of total new infections can then be represented by the following equation:
$$\begin{aligned} \int _{0}^{+\infty } F\mathscr {T}(t)(\zeta (\textbf{x})) \,dt. \end{aligned}$$
(17)
Next, consider the NGM operator K that takes the starting distribution of infections and describes how it spreads to the total number of infected individuals over time, as explained in33. It can be defined as:
$$\begin{aligned} K: \zeta (\textbf{x})\rightarrow F \int _{0}^{+\infty } \mathscr {T}(t)(\zeta (\textbf{x})) \,dt. \end{aligned}$$
(18)
Following19,24,35, the spectral radius of K yields the basic reproduction number of the PDE model:
$$\begin{aligned} \mathscr {R}_0 ^ {\text {PDE}} = \rho (K), \end{aligned}$$
(19)
where \(\rho\) is the spectral radius of the matrix. Moreover, Theorem 3.1 of Wang et al.19, (which follows from Theorem 3.12 of Thieme et al.36), provides the following equation:
$$\begin{aligned} K(\zeta )= F \int _{0}^{+\infty } \mathscr {T}(t)(\zeta ) \,dt = -F\Gamma ^{-1}(\zeta ), \end{aligned}$$
(20)
where \(\Gamma := \nabla . (d_\mathbb {I}\nabla V) – V\). The term \(\nabla . (d_\mathbb {I}\nabla V)\) represents the divergence of the gradient of V multiplied by a diffusion coefficient \(d_\mathbb {I}\).
Basic reproduction number: numerical formulation
We now present the numerical formulation—within a discretized setting—corresponding to the mathematical description of \(\mathscr {R}_0\) provided in “Basic reproduction number: analytical formulationR0spsanalytical”. Let \(\mu\) be an eigenvalue of K in Eq. (17) such that \(K(\zeta (\textbf{x}))= \mu \zeta (\textbf{x})\) with corresponding eigenvector \(\zeta (\textbf{x}) = (\zeta _1(\textbf{x}),…, \zeta _m(\textbf{x}))^T\). Then,
$$\begin{aligned} -F\Gamma ^{-1}(\zeta (\textbf{x}))= \mu \zeta (\textbf{x}). \end{aligned}$$
(21)
Suppose that \(V= diag(\nu _1,\ldots ,\nu _m)\), \(\psi (\textbf{x})\)= \(- \Gamma ^ {-1} (\zeta (\textbf{x}))\), where \(\psi (\textbf{x}) = (\psi _1(\textbf{x}),…, \psi _m(\textbf{x}))^T\), then \(-\Gamma (\psi (\textbf{x})) = \zeta (\textbf{x})\); i.e.,
$$\begin{aligned} \begin{aligned} – \left( d_i(\textbf{x}) \nabla ^2 \psi _i(\textbf{x}) – c_i(\textbf{x}) \nabla \psi _i(\textbf{x}) – \nu _i\psi _i(\textbf{x})\right) = \zeta _i(\textbf{x}), \quad 1 \le i \le m. \end{aligned} \end{aligned}$$
(22)
For numerical computation, we consider a square grid of size \(\mathscr {N} \times \mathscr {N}\). We consider a sufficiently large integer \(\mathscr {N}>0\) that represents the grid size of the computational domain. In the formulation below, we compute the reproduction number for each row indexed by \(\texttt {y}\). For each row \(\texttt {y}\), we consider a one-dimensional discretized spatial domain (normalized) along the x-axis [0, 1]. Here, \(x_j = j/\mathscr {N}\), \(d_{ij} = d_i(x_j)\), \(c_{ij} = c_i(x_j)\), \(\psi _{ij} = \psi _i(x_j)\) and \(\zeta _{ij} = \zeta _i(x_j)\) for \(j= 1,\ldots ,\mathscr {N}\). By applying a central difference scheme, we rewrite Eq. (22) as:
$$\begin{aligned} ( d_{i,j} \frac{\psi _{i,j+1}-2\psi _{i,j}+\psi _{i,j-1}}{1/\mathscr {N}^2} – c_{i,j}\frac{\psi _{i,j+1}-\psi _{i,j-1}}{2/\mathscr {N}} – \nu _i \psi _{i,j}) \approx \zeta _{i,j} \end{aligned}$$
(23)
for all \(1\le\)j\(\le \mathscr {N}\), with \(\psi _{i,1}= \psi _{i,2}, \psi _{i,\mathscr {N}}= \psi _{i ,\mathscr {N}-1}\) according to the no-flux, Neumann boundary conditions. By combining the \({\mathscr {N}}\) terms in a matrix, we obtain:
$$\begin{aligned} A_{i}\Psi _{i} \approx \zeta _{i}. \end{aligned}$$
(24)
where \(A= diag(A_1, \ldots , A_m)\), and:
$$\begin{aligned} \left\{ \begin{array}{ll} A_i & = \begin{bmatrix} a_{i,1} – 2d_{i,1} \mathscr {N}^2 & 0 & \cdots & 0 & 0\\ b_{i,2}^+a_{i,2} & b_{i,2}^- & \ddots & \vdots & \vdots \\ \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & \cdots & b_{i,\mathscr {N}-1}^+ & a_{i,\mathscr {N}-1} & b_{i,\mathscr {N}-1}^- \\ 0 & \cdots & 0 & -2d_{i,\mathscr {N}} \mathscr {N}^2 & a_{i,\mathscr {N}} \\ \end{bmatrix}, \\ \Psi _i & = \begin{bmatrix} \psi _{i,1} \\ \psi _{i,2} \\ \vdots \\ \psi _{i,\mathscr {N}-1} \\ \psi _{i,\mathscr {N}} \end{bmatrix}, \quad \text {and} \quad \zeta _i = \begin{bmatrix} \zeta _{i,1} \\ \zeta _{i,2} \\ \vdots \\ \zeta _{i,\mathscr {N}-1} \\ \zeta _{i,\mathscr {N}} \end{bmatrix} \end{array} \right\} \end{aligned}$$
(25)
with \(b_{i,j}^+ = -\mathscr {N}(d_{i,j}\mathscr {N} + \frac{c_{i,j}}{2}), \quad b_{i,j}^- = -\mathscr {N}(d_{i,j}\mathscr {N} – \frac{c_{i,j}}{2}), \quad \text {and} \quad a_{i,j} = 2d_{i,j}\mathscr {N}^2 + \nu _i; for \quad all \quad 0\le j \le \mathscr {N}\).
From Theorems 3.1 and 3.2 in Wang et al.24, it follows that:
$$\begin{aligned} F\psi (x) = -F\Gamma ^{-1}(\zeta (x))= \mu \zeta (x). \end{aligned}$$
(26)
The above equation may also be written as:
$$\begin{aligned} (F \bigotimes \mathscr {I}_{\mathscr {N}})\Psi = \mu \zeta . \end{aligned}$$
(27)
where \(\mathscr {I}_{\mathscr {N}}\) is the identity matrix of dimensions \((\mathscr {N}) \times (\mathscr {N})\) and \(\bigotimes\) is the Kronecker product. The standard definition of the Kronecker product is as follows: Given an \(r \times s\) matrix M with elements \(m_{ij}\) and another matrix Q of size \(p \times q\),
$$\begin{aligned} M \otimes Q = \begin{bmatrix} m_{11}Q & \dots & m_{1s}Q \\ \vdots & \ddots & \vdots \\ m_{r1}Q & \dots & m_{rs}Q \end{bmatrix} \end{aligned}$$
Substituting \(\Psi \approx A^{-1}\Phi\) in the above equation yields:
$$\begin{aligned} (F \bigotimes \mathscr {I}_{\mathscr {N}})A^{-1}\zeta \approx \mu \zeta . \end{aligned}$$
(28)
Therefore, we can now state the following central result:
$$\begin{aligned} \mathscr {R}_{0(\texttt {y})}^{PDE}= \rho ((F \bigotimes \mathscr {I}_{\mathscr {N}})A^{-1}). \end{aligned}$$
(29)
We note the introduction of row index \(\texttt {y}\) in the subscript above to emphasize that the above equation computes the reproduction number for row \(\texttt {y}\).
As discussed previously, at the onset of an epidemic, \(\mathscr {R}_0\) acts as a marker of the spread intensity, predicated on the assumption that the entire population is susceptible. However, as the epidemic progresses, this assumption becomes tenuous, owing to factors including the acquisition of immunity, behavioural modifications, and spatial heterogeneities. In other words, \(\mathscr {R}_0\) is an inadequate measure of the intensity of spread beyond the initial stages of an active epidemic37. This necessitates the introduction of \(\mathscr {R}_e\), which reflects the actual average number of secondary cases per infectious case at a given time and space for an evolving epidemic. Calculating \(\mathscr {R}_e\) is crucial for understanding the current state of epidemic spread at any point in time, as well as for assessing the impact of public health interventions and making informed decisions about future strategies38. It helps capture the dynamic nature of the epidemic by accounting for changes in susceptibility, contact patterns, and mobility, thus providing a more timely and context-specific metric. Referring to1,5,39, we note that:
$$\begin{aligned} \mathscr {R}_{e(\texttt {y})}^{PDE} = \mathscr {R}_{0(\texttt {y})}^{PDE} \times (S/P). \end{aligned}$$
(30)
where S = The total number of susceptible individuals at a point in time, and P = The total population at that time.
Reproduction number for the discretized SLIR model
In the light of the discussion in the previous section, we now focus on the reproduction number \(\mathscr {R}_e\) for the discretized version (Eqs. 5–8) of our underlying SLIR model. At the outset, we note that the two infected compartments in our model are the Latent and the Infected ones given in (6) and (7), respectively. We recall that when an individual gets infected, they first migrate into the Latent compartment (where they are considered to be infected but not infectious). Once they are fully infectious after a specific incubation period \((\tau )\), they transition into the Infected compartment.
Next, we recall that in the model, all new infections are accounted for in the F matrix. The F and V matrices can be derived by linearizing the discretized equations of the model; the linearization may be achieved by the standard analytical method of computing the Jacobian (first-order partial derivatives of a vector function) of each equation relative to the others at an infection-free equilibrium state denoted by \((S^0,L^0,I^0,R^0)= (S^0,0,0,0)\). The infection-free steady state, as discussed previously, represents the initial scenario when the entire population is susceptible (denoted by \(S^0\)), and, therefore, there are no individuals in the Latent, Infected or Recovered states. The linearization procedure facilitates the estimation of \(\mathscr {R}_0\) as described analytically in the previous sections. Next, we determine the F and V matrices for our underlying model, as outlined in the Eqs. (6) and (7). As the first step in this process, we compute the \(\mathscr {F}\) and \(\mathscr {V}\) matrices, defined in Eq. (12), as:
$$\begin{aligned} & \mathscr {F}= \begin{bmatrix} \frac{\phi }{N}I^{T}_{x,y}S^{T}_{x,y} \\ 0 \\ \end{bmatrix} \end{aligned}$$
(31)
$$\begin{aligned} & \mathscr {V}= \begin{bmatrix} \eta _{L}(L^{T}_{x,y + 1} + L^T_{x,y -1} + L^{T}_{x+1,y} + L^{T}_{x-1,y})+\\ (1-4\eta _{L})L^{T}_{x,y} – \epsilon \frac{\phi }{N}\sum _{x=1}^{\mathscr {N}}\sum _{y=1}^{\mathscr {N}}I^{T – \tau }_{x,y}S^{T-\tau }_{x,y}f_{a}l^{2} \\ \\ \eta _{I}(I^{T}_{x,y + 1} + I^T_{x,y -1} + I^{T}_{x+1,y} + I^{T}_{x-1,y})+\\ (1-4\eta _{I} – \delta )I^{T}_{x,y} + \epsilon \frac{\phi }{N}\sum _{x=1}^{\mathscr {N}}\sum _{y=1}^{\mathscr {N}}I^{T – \tau }_{x,y}S^{T-\tau }_{x,y}f_{a}l^{2}\\ \end{bmatrix} \end{aligned}$$
(32)
From Eqns. 14 and 15 , we obtain the matrices F and V as:
$$\begin{aligned} & F= \begin{bmatrix} 0& \frac{\phi }{N} S^0\\ 0& 0 \end{bmatrix} \end{aligned}$$
(33)
$$\begin{aligned} & V= \begin{bmatrix} 1-4\eta _{L}& – \epsilon \frac{\phi }{N}\sum _{x=1}^{\mathscr {N}}\sum _{y=1}^{\mathscr {N}}S^{T-\tau }_{x,y}f_{a}l^{2}\\ \\ 0& (1-4\eta _{I} – \delta ) + \epsilon \frac{\phi }{N}\sum _{x=1}^{\mathscr {N}}\sum _{y=1}^{\mathscr {N}}S^{T-\tau }_{x,y}f_{a}l ^{2} \end{bmatrix} \end{aligned}$$
(34)
We note that the disease-free equilibrium, \(S^0\), represents a steady state where the number of individuals becoming susceptible (e.g., through birth, represented by \(\lambda )\) is balanced by the number of individuals leaving the susceptible state (e.g., through natural death, which is represented by \(\theta\)). The equilibrium value of \(S^0\) in Eq. (33) is therefore determined by the ratio of these rates10,24,40:
$$\begin{aligned} S^0 = \lambda / \theta . \end{aligned}$$
(35)
The numerical formulation presented in “Basic reproduction number: numerical formulation” is now used in a grid framework where x and y represent the row and column indices. \(L_{x,y}\) represents the Latent population for the grid element positioned at row x and column y. We now apply Eq. (23) to obtain the Latent population \(L_{x,y}\) and arrive at the following equation:
$$\begin{aligned} \begin{array}{c} -\mathscr {N}(d_L*\mathscr {N})L_{x,y-1}+(2(d_L)\mathscr {N}^2 – V_{L})L_{x,y} -\mathscr {N}(d_L*\mathscr {N})L_{x,y+1} \approx L_{x,y} \end{array} \end{aligned}$$
(36)
where \(V_L\), which represents the contributions to V (see Eqn: 34) from the Latent compartment, is given by:
$$\begin{aligned} V_{L} = (1-4\eta _{L}) – \epsilon \frac{\phi }{N}\sum _{x=1}^{\mathscr {N}}\sum _{y=1}^{\mathscr {N}}S^{T-\tau }_{x,y}f_{a}l^{2} \end{aligned}$$
(37)
We also note, in deriving 36, that we assume the diffusion coefficient for the Latent population to be spatially uniform and given by \(d_L\). From Eq. (23), we note that \(c_L\) vanishes under the above assumption. Therefore, the \(A_L\) matrix is obtained as:
$$\begin{aligned} A_L = \begin{bmatrix} a_{L,1} & -2d_{L,1}\mathscr {N}^{2} & & \\ b_{L,2}^{+} & a_{L,2} & b_{L,2}^{-} & & \\ \ddots & \ddots & \ddots & & & \\ & & b_{L,\mathscr {N}-1}^{+} & a_{L,\mathscr {N}-1} & b_{L,\mathscr {N}-1}^{-} \\ & & & b_{L,\mathscr {N}}^{+} & a_{L,\mathscr {N}} \\ \end{bmatrix} \end{aligned}$$
(38)
where,
$$\begin{aligned} \begin{array}{c} b_{L,y}^{+} = -\mathscr {N}(d_{L,y} * \mathscr {N}) \\ a_{L,y} = (2d_{L,y}\mathscr {N}^{2} + V_{L} ) \\ b_{L,y}^{-} = -\mathscr {N}(d_{L,y} * \mathscr {N}) \\ \end{array} \end{aligned}$$
(39)
Similarly, we now determine the matrix \(A_I\). First, we compute the \(I_{x,y}\) using the Eq. (23):
$$\begin{aligned} \begin{array}{c} -\mathscr {N}(d_I*\mathscr {N})I_{x,y-1}+(2(d_I)\mathscr {N}^2 – V_{I})I_{x,y} -\mathscr {N}(d_I*\mathscr {N})I_{x,y+1} \approx I_{x,y}, \end{array} \end{aligned}$$
(40)
where \(V_I\), which represents the contributions in V (see Eq. (34)) from Infected compartment, is given by:
$$\begin{aligned} V_{I} = (1-4\eta _{I} – \delta ) + \epsilon \frac{\phi }{N}\sum _{x=1}^{\mathscr {N}}\sum _{y=1}^{\mathscr {N}}S^{T-\tau }_{x,y}f_{a}l ^{2} \end{aligned}$$
(41)
In the derivation of Eq. (40), we note that the diffusion coefficient for the Infected compartment (\(d_I\)) is assumed to be spatially uniform. This assumption leads to the elimination of the term \(c_I\) in Eq. (23). Therefore, the \(A_I\) matrix is given by:
$$\begin{aligned} A_I = \begin{bmatrix} a_{I,1} & -2d_{I,1}\mathscr {N}^{2} & & \\ b_{I,2}^{+} & a_{I,2} & b_{I,2}^{-} & & \\ \ddots & \ddots & \ddots \\ & & b_{I,\mathscr {N}-1}^{+} & a_{I,\mathscr {N}-1} & b_{I,\mathscr {N}-1}^{-} \\ & & & b_{I,\mathscr {N}}^{+} & a_{I,\mathscr {N}} \\ \end{bmatrix} \end{aligned}$$
(42)
where,
$$\begin{aligned} \begin{array}{c} b_{I,y}^{+} = -\mathscr {N}(d_{I,y} * \mathscr {N}) \\ a_{I,y} = (2d_{I,y}\mathscr {N}^{2} + V_{I}) \\ b_{I,y}^{-} = -\mathscr {N}(d_{I,y} * \mathscr {N}) \\ \end{array} \end{aligned}$$
(43)
Therefore, the matrix A is given by:
$$\begin{aligned} A= \begin{bmatrix} A_L & \\ & A_I \\ \end{bmatrix} \end{aligned}$$
(44)
By substituting the obtained expressions for the matrix A (Eq. 44) and the matrix F (Eq. 33) into Eq. (28) and solving the eigenvalue problem described by Eq. (28), one obtains the principal eigenvalue \(\rho\) in the equation (Eq. 29). From there, we find the \(\mathscr {R}_{e(\texttt {y})}^{PDE}\) of the system by using Eq. (30) that provides the reproduction number of for row \(\texttt {y}\) at any specified time. The overall reproduction number of the entire region is then computed by averaging over all \(\mathscr {N}\) rows using:
$$\begin{aligned} \mathscr {R}_e = \frac{\sum _{\texttt {y}=1}^{\mathscr {N}} \mathscr {R}_{e(\texttt {y})}^{PDE}}{\mathscr {N}} \end{aligned}$$
(45)
where, ’\(\mathscr {N}\)’ is the grid size of the spatial domain, \(\mathscr {R}_{e(\texttt {y})}^{PDE}\) is the reproduction number for the row \((\texttt {y})\), and \(\mathscr {R}_e\) is the reproduction number for the entire domain.
The Wallinga–Teunis method
To carry out a comparative study with respect to established (albeit retrospective) methods for calculating the reproduction number, we first choose the Wallinga & Teunis (W*-T) method25. To determine the reproduction number for a primary case, the W–T method employs the probability distribution for the generation interval (also called the serial interval), which is the time from symptom onset in a primary case to symptom onset in a secondary case. By using this distribution, the method estimates the relative likelihood that a secondary case has been infected by a specific primary case, given the observed difference in their time of symptom onset. The effective reproduction number is then calculated by summing over all observed secondary cases, each weighted by its relative likelihood of being infected by the particular primary case. Our calculations applying the W–T method were carried out through the R package, \(EpiEstim v2.2-4\)26. This framework assumes that there is a complete recording of cases and that there are no asymptomatic cases. The W–T method is completely retrospective, meaning it contains a sliding window that takes in the data for the previous 30 days as input to provide the \(\mathscr {R}_e\) estimates for a given day.
We choose the Wallinga–Teunis (W–T) method as a validation tool for our study as it is a well-established method that has been widely used in the field of epidemiology and by public health agencies, including the Center for Disease Control (CDC), the World Health Organization (WHO), as well in academic research. In addition to being used in multiple studies during the COVID-19 pandemic (see, for instance,41,42), it was also used during the previous SARS epidemic and for other infectious respiratory disease epidemics to gain a comprehensive understanding of the transmission dynamics over time. For our purposes, one of the main strengths of the W–T method is that it does not assume a constant serial interval, allowing for dynamic changes over time resulting from interventions and changes in the pathogen (virus). It is also robust when the incidence is low, and it considers the entire epidemic curve, thus providing a good benchmark41,42. Moreover, it relies on retrospectively collected data. For those reasons, we employ the W–T method to validate the prediction performance of our framework.
Cori method
To further buttress the validation, we also employ the Cori method—yet another widely adopted, retrospective method in epidemiology for reproduction number computations26. The Cori method focuses on obtaining instantaneous values of\(\mathscr {R}_{e}\). This method assumes that the random spread of disease over time follows a Poisson process with mean \(R_t \sum _{s=1}^{t} I_{t-s} w_s\), where \(I_{t-s}\) is the incidence infected in time step \((t-s\)) and \(w_s\) is a probability distribution describing the average infectiousness profile after infection, approximated by the distribution of the serial interval.
The infectiousness profile of individuals-indicating how contagious they are after becoming infected-is modeled using a distribution of the serial interval. The model assumes that the transmissibility of the disease (i.e., the ease of spread) remains constant over a fixed period. To estimate the reproduction number (\(\mathscr {R}_{e}\)), this method follows a Bayesian framework. In this framework, prior information about transmissibility is represented by a Gamma prior distribution, which is based on previous knowledge or appropriate assumptions.
Post observing new data, the model updates its estimate of the reproduction number using this Bayesian estimation approach. The updated estimate (posterior distribution) for the reproduction number also follows a Gamma distribution but with new parameters that reflect the observed data.
Similar to the W–T method, the latest R package of \(EpiEstim v2.2-4\)26, developed by the Cori team at Imperial College, London, UK43 was employed. We ran the program on R4.2.2 (64 bit), using the function “estimate_R”for the Cori method. We calculated the \(\mathscr {R}_{e}\) values over sliding weekly windows in each timeframe.
