Reaction–Diffusion Equations in Mathematical Models Arising in Epidemiology

Abstract

The review is devoted to an analysis of mathematical models used for describing epidemic processes. Our main focus is on the models that are based on partial differential equations (PDEs), especially those that were developed and used for the COVID-19 pandemic modeling. Most of our attention is given to the studies in which not only results of numerical simulations are presented but analytical results as well. In particular, traveling fronts (waves), exact solutions, and the estimation of key epidemic parameters of the epidemic models with governing PDEs (typically reaction–diffusion equations) are discussed. The review may serve as a valuable resource for researchers and practitioners in the field of mathematical modeling in epidemiology.

1. Introduction

At the present time, there are numerous mathematical models describing epidemic processes that can be found in several books dedicated to mathematical modeling in life sciences (see [1,2,3,4,5,6,7] and papers cited therein). It is widely thought that the Kermack–McKendrick study [8] is a pioneering paper in this direction. The authors created a model based on three ordinary differential equations (ODEs). Nowadays, their model is called the Susceptible–Infectious–Recovered (SIR) model. There are several generalizations of the SIR model such as the Susceptible–Exposed–Infectious–Recovered (SEIR) model suggested in [9,10] and the Susceptible–Infectious–Recovered–Deceased (SIRD) model [7,11,12]. Some other models (see, e.g., [13,14,15]) were developed after the outbreak of the COVID-19 coronavirus, because this novel pandemic has attracted the extensive attention of many mathematicians working in the field of mathematical modeling.

The classical pandemic models, first of all SIR, SEIR and SIRD, have been extensively used in the modeling of the COVID-19 pandemic in order to understand and predict the transmission dynamics of the disease (see, e.g., [16,17,18,19,20,21]). In particular, these models and their variants can be used to estimate key epidemiological parameters of COVID-19, such as the basic reproduction number $R_0$ (see, e.g., [18,22,23]). The number $R_0$ represents an average number of new infections caused by a single infected individual in the susceptible population.

Typically, numerical simulations are used for solving the above-mentioned models [24,25,26]. It should be noted that the numerical solving of ODEs is not a cutting-edge problem at the present time, because there are many computer program packages adopted for these purposes. However, exact solving of the well-known models (SIR, SEIR, SIRD, etc.) is a nontrivial problem because the relevant ODEs are nonlinear. There are several studies focused on investigating and constructing exact solutions of the classical SIR model as well as its generalizations, modifications, and simplification (see, [14,15,27,28,29,30] and papers cited therein). There are also studies focused on solving such ODE systems using approximate techniques. For instance, the SIR epidemic model is solved by the homotopy analysis method, and solutions are derived in the form of series involving exponents in [31]. In [32], the SEIR model is studied, and approximate solutions are obtained using a so-called optimal auxiliary functions method.

It is worth noting that the propagation of numerous epidemic processes, including the COVID-19 pandemic, often exhibits significant spatial heterogeneity. For example, an essential spatial heterogeneity was observed in large European countries (Italy is a typical example) and the USA during the first pandemic wave. This observation can be interpreted in various manners, but the prevailing approach involves partitioning the larger spatial domain (such as a country) into multiple smaller sub-domains (regional divisions) and employing standard ODE-based models to each sub-domain. Nevertheless, an alternative method exists—using the reaction–diffusion equations—to model the diffusion-based spread of the infected population. Reaction–diffusion systems are used to describe the dynamics of spatially extended systems in which the interactions between components involve both reaction and diffusion processes. These models have been widely applied in various fields, including physics, chemistry, biology, and epidemiology, in order to understand the spread and behavior of populations or substances. When considering epidemic modeling, reaction–diffusion systems can be used to study the spatial spread of infectious diseases (see pioneering [33,34,35,36,37] and recent works [38,39,40,41,42,43,44,45,46,47,48]).

There are also epidemic models that involve more complicated equations, such as models based on reaction–diffusion systems with cross-diffusion [7,34,49,50,51]. It should be noted that cross-diffusion phenomena occur in other biological processes and were introduced in the 1970s, in particular, for mathematical modeling in chemotaxis [52] and population dynamics [53]. Nowadays, such mathematical models are intensively studied by different mathematical techniques (see, e.g., [54] and references therein).

Age-structured epidemic models should be mentioned as well. Such models are based on integro-differential equations in order to incorporate both spatial and age dimensions for modeling the spread of infectious diseases (see [55,56,57,58] and papers cited therein). These models aim to understand the complex dynamics of disease transmission in relation to age-specific factors and spatial patterns.

This review is organized as follows. In Section 2, we briefly present classical models (SIR, SEIR, etc.) describing epidemic processes with the stress on the integrability of such models. In Section 3, we present epidemic models based on systems of reaction–diffusion equations. Some of them were developed before the COVID-19 pandemic; others are suggested very recently. In Section 4, models involving PDEs with convective terms and those with cross-diffusion are discussed. In Section 5, age-structured epidemic models are considered. Finally, we present conclusions and briefly discuss unsolved problems in epidemic modeling in the last section.

2. Integrability of the Classical Epidemic Models

As mentioned above, the SIR model is a pioneering model for describing epidemic processes. The model is based on the three-component system of ODEs that reads as [8]

$$\begin{array}{c}\frac{dS}{dt}=-\alpha SI,\hfill \\ \frac{dI}{dt}=\alpha SI-\beta I,\hfill \\ \frac{dR}{dt}=\beta I.\hfill \end{array}$$

(1)

The model divides the total population into three subpopulations (compartments): susceptible (S), infectious (I) (with infectious capacity and not yet recovered) and recovered (R) (recovered and neither infectious nor infected once again). In (1), $\alpha >0$ is the transmission rate (represents the average rate at which susceptible individuals become infected when they come into contact with infectious individuals); $\beta >0$ is the recovery rate (represents the average rate at which infectious individuals recover from the disease and gain immunity). The main assumption of the SIR model consists of the conservation law of the total population. Mathematically, it directly follows from (1) if one takes the sum of all equations

$$\frac{dS}{dt}+\frac{dI}{dt}+\frac{dR}{dt}=0,$$

(2)

i.e., the total population $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)=N_0$ ($N_0$ is a fixed number). It means that the epidemic process does not lead to deaths, i.e., the death rate is zero (each member of the inflicted subpopulation will survive). Obviously, this assumption is rather unrealistic for such processes as the COVID-19 pandemic. It should be also noted that the SIR model neglects the natural birth/death rate.

The general solution of the nonlinear ODE system (1) cannot be expressed explicitly, but it is well known in parametric form (see, e.g., [28]):

$$S\left(t\right)=S_0\tau ,I\left(t\right)=\frac{\beta }{\alpha }ln\tau -S_0\tau -\frac{I_0}{\alpha },R\left(t\right)=-\frac{\beta }{\alpha }ln\tau ,$$

where $S_0$ and $I_0$ are integration constants, while the parameter $\tau$ is defined by the integral

$$\int \frac{d\tau }{\tau \left(I_0-\beta ln\tau +\alpha S_0\tau \right)}=t.$$

(3)

The above integral leads to special functions; therefore, a transcendent functional equation is obtained for finding the parameter $\tau$. Obviously, $S_0$ and $I_0$ allow us to satisfy initial conditions for the functions $S\left(t\right)$ and $I\left(t\right)$ (the third constant $R_0$ springs up from the above integral).

In paper [28], the authors also consider the SIR model with birth and death rates (see (4)–(6) [28]) and show that the relevant nonlinear system of ODEs can be reduced to the Abel equation, which is not solvable. However, the authors studied the equation obtained using semianalytical/numerical methods.

The most known generalization of the SIR model that assumes nonzero death rate for the infected subpopulation is the SIRD model. The SIRD model takes into account the following assumptions: there is no cure or immunity and some infected members of the subpopulation I will die at a certain rate; the recovered subpopulation will not remain immune and can be again infected (in contrast to the SIR model where the recovered group obtains immunity from the disease). Consequently, the SIRD model is described by the following system of differential equations (see, e.g., [12]):

$$\begin{array}{c}\frac{dS}{dt}=-\alpha SI,\hfill \\ \frac{dI}{dt}=\alpha SI-\beta I-\mu I,\hfill \\ \frac{dR}{dt}=\beta I,\hfill \\ \frac{dD}{dt}=\mu I,\hfill \end{array}$$

(4)

where D represents the number of individuals who have died due to the disease; new parameter $\mu$ is the death rate.

In [29], an exact solution of the SIRD epidemic model (4) is constructed for arbitrary initial conditions

$$S\left(0\right)=S_0,I\left(0\right)=I_0,R\left(0\right)=R_0,D\left(0\right)=D_0,S_0+I_0+R_0+D_0=N_0$$

in the parametric form

$$\begin{array}{c}S\left(t\right)=S_{0}exp\left(\frac{\alpha }{\beta }{R}_0\right)\tau ,\hfill \\ I\left(t\right)=N_0-D_0+\frac{\mu }{\beta }{R}_0-S_{0}exp\left(\frac{\alpha }{\beta }{R}_0\right)\tau +\frac{\beta +\mu }{\alpha }ln\tau ,\hfill \\ R\left(t\right)=-\frac{\beta }{\alpha }ln\tau ,\hfill \\ D\left(t\right)=D_0-\frac{\mu }{\alpha }ln\tau -\frac{\mu }{\beta }{R}_0,\hfill \end{array}$$

(5)

where the parameter $\tau \in \left(exp\left(-\frac{\alpha }{\beta }R(\infty )\right),exp\left(-\frac{\alpha }{\beta }{R}_0\right)\right],$

$$t={\int }_{\tau }^{exp\left(-\frac{\alpha }{\beta }{R}_0\right)}\frac{dz}{z\psi \left(z\right)},$$

while $\psi \left(z\right)=\alpha (N_0-D_0)+\frac{\alpha \mu }{\beta }{R}_0-\beta S_{0}exp\left(\frac{\alpha }{\beta }{R}_0\right)z+(\beta +\mu )lnz.$ Notably, the above formula for the parameter $\tau$ is similar to (3).

One notes that the sum of all components in (5) is $N_0$, i.e., the conservation law is again preserved. However, the total population in this case (in contrast to the SIR model) includes those who died from the epidemic disease.

The natural generalization of the SIR model is the SEIR one

$$\begin{array}{c}\frac{dS}{dt}=-\alpha SI,\hfill \\ \frac{dE}{dt}=\alpha SI-\gamma E,\hfill \\ \frac{dI}{dt}=\gamma E-\beta I,\hfill \\ \frac{dR}{dt}=\beta I.\hfill \end{array}$$

(6)

Model (6) includes a new function $E\left(t\right)$ for the exposed subpopulation, representing individuals who have been exposed to the infection but are currently in a latent period and not yet capable of transmitting the disease. In (6), the coefficient $\gamma >0$ denotes the transition rate from exposed individuals to the infected one and determines an incubation period $\frac{1}{\gamma }$ that represents the average duration from the time of infection to the point at which an individual becomes capable of transmitting the disease to others.

Clearly, the SEIR model is developed also under the assumption that the conservation law of the total population takes place. One may say that the SEIR model is a straightforward generalization of the SIR model.

In recent work [30], an exact solution of the SEIR model (6) with the initial conditions

$$S\left(0\right)=S_0,E\left(0\right)=E_0,I\left(0\right)=I_0,R\left(0\right)=R_0,S_0+E_0+I_0+R_0=N_0$$

is constructed. Similarly to the SIR model case, the exact solution is found in the parametric form:

$$\begin{array}{c}S\left(t\right)=S_{0}exp\left(\frac{\alpha }{\beta }{R}_0\right)\tau ,\hfill \\ E\left(t\right)=E_{0}exp\left(-\gamma \phi \left(\tau \right)\right)+S_{0}exp\left(\frac{\alpha }{\beta }{R}_0-\gamma \phi \left(\tau \right)\right){\int }_{\tau }^{exp\left(-\frac{\alpha }{\beta }{R}_0\right)}exp\left(\gamma \phi \left(z\right)\right)dz,\hfill \\ I\left(t\right)=N_0-S_{0}exp\left(\frac{\alpha }{\beta }{R}_0\right)\tau +\frac{\beta }{\alpha }ln\tau -E_{0}exp\left(-\gamma \phi \left(\tau \right)\right)-\hfill \\ S_{0}exp\left(\frac{\alpha }{\beta }{R}_0-\gamma \phi \left(\tau \right)\right){\int }_{\tau }^{exp\left(-\frac{\alpha }{\beta }{R}_0\right)}exp\left(\gamma \phi \left(z\right)\right)dz,\hfill \\ R\left(t\right)=-\frac{\beta }{\alpha }ln\tau ,\hfill \end{array}$$

(7)

where the parameter $\tau \in \left(exp\left(-\frac{\alpha }{\beta }R(\infty )\right),exp\left(-\frac{\alpha }{\beta }{R}_0\right)\right],$

$$t=\phi \left(\tau \right)={\int }_{\tau }^{exp\left(-\frac{\alpha }{\beta }{R}_0\right)}\frac{dz}{z\psi \left(z\right)}.$$

Here, the function $\psi \left(z\right)$ is the solution of the Abel equation of the second kind

$$z\psi {\psi }^{\prime }-(\beta +\gamma )\psi +\beta \gamma lnz-\alpha \gamma exp\left(\frac{\alpha }{\beta }{R}_0\right)z+\alpha \gamma N_0=0,$$

on the interval

$$\left(exp\left(-\frac{\alpha }{\beta }R(\infty )\right),exp\left(-\frac{\alpha }{\beta }{R}_0\right)\right),$$

that satisfies the boundary conditions

$$\begin{array}{c}\psi \left(exp\left(-\frac{\alpha }{\beta }{R}_0\right)\right)=\beta I_0,\hfill \\ \underset{z\to exp\left(-\frac{\alpha }{\beta }R(\infty )\right)+0}{lim}\psi \left(z\right)=0.\hfill \end{array}$$

Moreover, a natural requirement about the positivity of the solution should be fulfilled

$$\psi \left(z\right)>0\mathtt{in}\left(exp\left(-\frac{\alpha }{\beta }R(\infty )\right),exp\left(-\frac{\alpha }{\beta }{R}_0\right)\right].$$

An applicability of the exact solution (7) for practical applications is questionable because that is too cumbersome.

The natural generalization of the SEIR model (6) that assumes a nonzero death rate for the infected subpopulation is the Susceptible–Exposed–Infected–Recovered–Deceased (SEIRD) model

$$\begin{array}{c}\frac{dS}{dt}=-\alpha SI,\hfill \\ \frac{dE}{dt}=\alpha SI-\gamma E,\hfill \\ \frac{dI}{dt}=\gamma E-\beta I-\mu I,\hfill \\ \frac{dR}{dt}=\beta I,\hfill \\ \frac{dD}{dt}=\mu I,\hfill \end{array}$$

where all parameters have the same interpretations as above. To the best of our knowledge, there are no papers devoted to the search for exact solutions of the SEIRD-type models. However, these models were used in several papers (see, e.g., [26,59]) for modeling the COVID-19 pandemic. Typically, results of numerical simulations are presented in such papers that describe the dynamics of the pandemic and suggest control strategies.

There are some other models describing epidemic processes that are based on ODEs systems; however, those cannot be considered as direct generalizations of the SIR model. Interestingly, some of them are integrable, and here we present examples.

In [27], the authors investigate the model

$$\begin{array}{c}\frac{dS}{dt}=-\alpha SI-\mu S+\beta I+\mu K,\hfill \\ \frac{dI}{dt}=\alpha SI-(\mu +\beta )I,\hfill \end{array}$$

(8)

that can be thought of as a simplification of the classical SIR model. However, it is rather difficult to identify assumptions that reduce the SIR model to the ODE system (8). The authors assume that the parameter $\mu$ is the proportionate death rate, while the term $\mu K$ represents a constant birth rate. The integrability of this model is proved by means of the Painlevé analysis. Moreover, using the Lie symmetry analysis, the model was completely integrated. As a result, the following exact solution in terms of elementary functions was constructed

$$\begin{array}{c}S\left(t\right)=\frac{\mu +\beta }{\alpha }+I\left(t\right)\left(\mu +\beta -\alpha K+I\left(t\right)-\mu B{e}^{-\mu t}\right),\hfill \\ I\left(t\right)=\frac{1}{\alpha }\frac{Aexp\left[-(\mu +\beta -\alpha K)t+B{e}^{-\mu t}\right]}{A\mu \int exp\left[-(\mu +\beta -\alpha K)t+B{e}^{-\mu t}\right]dt+C}.\hfill \end{array}$$

Here $A,B$ and C are arbitrary constants. However, it can be noted that one of them can be skipped because three constants cannot vanish simultaneously. Two remaining constants can be used in order to satisfy initial conditions.

In [14,15], we proposed the model

$$\begin{array}{c}\frac{du}{dt}=u(a-b{u}^{\kappa }),\hfill \\ \frac{dv}{dt}=k\left(t\right)u,\hfill \\ u\left(0\right)=u_0\ge 0,v\left(0\right)=v_0\ge 0,\hfill \end{array}$$

(9)

for a quantitative description of the outbreak of the COVID-19 pandemic. In (9), a smooth function $u\left(t\right)$ presents the total number of the COVID-19 cases identified up to day (the time moment) t; $v\left(t\right)$ is the total number of deaths up to the time moment t. Typically, t is an integer number, but we assume that $u\left(t\right)$ and $v\left(t\right)$ are continuous functions similarly to the functions $S,E,I,R$ and D used above. One may also note that the relations between $u\left(t\right)$ and $v\left(t\right)$:

$$u\left(t\right)={\int }_0^{t}I\left(\tau \right)d\tau ,v\left(t\right)={\int }_0^{t}D\left(\tau \right)d\tau ,$$

assuming that the pandemic started at the moment $t=0$.

In (9), $a>0$ is the coefficient for the virus transmission mechanism; $b>0$ is the coefficient for the effectiveness of the government restrictions (quarantine rules); $\kappa >0$ is the exponent, which guarantees that the total number of the COVID-19 cases is bounded in time; the smooth function $k\left(t\right)>0$ is the coefficient for effectiveness of the healthcare system during the epidemic process. From a mathematical point of view, coefficient $k\left(t\right)$ should have the asymptotic behavior $k\left(t\right)\to 0$, if $t\to \infty$; otherwise, all infected people will die. It is assumed that $k\left(t\right)=k_{0}exp(-\alpha t),\alpha >0$. Note that in the case where $\kappa =1$, the first equation of model (9) coincides with the classical logistic equation [60] that occurs naturally in epidemiology, as it was shown under some general assumptions in [1].

The general solution of model (9) is constructed explicitly in the form

$$\begin{array}{c}u\left(t\right)=a^{1/\kappa }{u}_0{e}^{at}{\left(a+b{u}_0^{\kappa }(e^{a\kappa t}-1)\right)}^{-1/\kappa },\hfill \\ v\left(t\right)=a^{1/\kappa }{k}_0{u}_0{\int }_0^t{e}^{(a-\alpha )\tau }{\left(a+b{u}_0^{\kappa }(e^{a\kappa \tau }-1)\right)}^{-1/\kappa }d\tau +v_0.\hfill \end{array}$$

(10)

The integral in (10) can be expressed via special functions for arbitrary parameters a, $\alpha$ and $\kappa$. However, one can be expressed in terms of elementary functions in some specific cases. For example, one obtains

$$v\left(t\right)=\frac{2k_0\sqrt{u_0}}{\sqrt{b(a-b{u}_0)}}\left(arctan\left(\frac{\sqrt{b{u}_0}}{\sqrt{a-b{u}_0}}{e}^{\frac{at}{2}}\right)-arctan\left(\frac{\sqrt{b{u}_0}}{\sqrt{a-b{u}_0}}\right)\right)+v_0$$

in the case $2\alpha =a,\kappa =1$.

In [14,15], it was demonstrated that the nonlinear system (9) with correctly specified parameters and given initial conditions can be successfully used for describing the first wave of the COVID-19 pandemic in several countries (China, Austria, and France). In particular, it was established that the exponent $\kappa$ takes different values in different countries. For example, $\kappa =1$ (i.e., the case of the logistic equation in (9)) leads to very good correspondence between the exact solution (10) and measured data taken from [61]. However, the parameter $\kappa$ was essentially smaller for many countries in Europe during the first wave of the COVID-19 pandemic, for example, $\kappa =0.4$ for Austria and France.

Now, we point out that the model (9) was constructed under essential simplifications of the epidemic process in question. In particular, the model implicitly admits that $u\gg v$. On the other hand, it is well known that the COVID-19 outbreak in several countries was so severe that the mortality rate was rather high, i.e., the assumption $u\gg v$ is not true. In such cases, the model (9) can be generalized as follows

$$\begin{array}{c}\frac{du}{dt}=(u-v)(a-b{(u-v)}^{\kappa }),\hfill \\ \frac{dv}{dt}=k\left(t\right)(u-v),\hfill \\ u\left(0\right)=u_0\ge 0,v\left(0\right)=v_0\ge 0,\hfill \end{array}$$

(11)

In fact, the time evolution of the function u cannot depend on the infected persons who already died. Similarly, the number of new deaths cannot depend on the people who already died. Taking into account the equality $u-v=w$, where w is the total number of recovered persons, the nonlinear model (11) is reducible to the form

$$\begin{array}{c}\frac{dw}{dt}=w(a-k\left(t\right)-b{w}^{\kappa }),\hfill \\ \frac{dv}{dt}=k\left(t\right)w,\hfill \\ w\left(0\right)=w_0=u_0-v_0\ge 0,v\left(0\right)=v_0\ge 0,\hfill \end{array}$$

and its exact solution can be expressed in the explicit form

$$\begin{array}{c}w\left(t\right)=exp\left(at-{\int }_0^{t}k\left(\tau \right)\mathrm{d}\tau \right){\left(w_0^{-\kappa }+b\kappa {\int }_0^{t}exp\left(\kappa a\tau -\kappa {\int }_0^{\tau }k\left(z\right)\mathrm{d}z\right)\mathrm{d}\tau \right)}^{-\frac{1}{\kappa }},\hfill \\ v\left(t\right)=v_0+{\int }_0^{t}k\left(\tau \right)w\left(\tau \right)\mathrm{d}\tau .\hfill \end{array}$$

(12)

Obviously, solution (12) leads to the solution

$$\begin{array}{c}u\left(t\right)=v_0+w\left(t\right)+{\int }_0^{t}k\left(\tau \right)w\left(\tau \right)\mathrm{d}\tau ,\hfill \\ v\left(t\right)=v_0+{\int }_0^{t}k\left(\tau \right)w\left(\tau \right)\mathrm{d}\tau .\hfill \end{array}$$

of the model (11).

3. Classical Epidemic Models with Diffusion in Space

A natural generalization of the classical Kermack–McKendric model that takes into account the diffusion process reads as

$$\begin{array}{c}\frac{\partial s}{\partial t}=d_s\Delta s-\alpha si,\hfill \\ \frac{\partial i}{\partial t}=d_i\Delta i+\alpha si-\beta i,\hfill \\ \frac{\partial r}{\partial t}=d_r\Delta r+\beta i,\hfill \end{array}$$

(13)

where $s(t,x)$, $i(t,x)$ and $r(t,x)$ are the susceptible, infective and recovered population densities at time t in position $x\in {\mathbb{R}}^n,n=1,2,3$, respectively. Here, $\Delta$ denotes the Laplace operator, and $d_s,d_i$ and $d_r$ are diffusion constants. Obviously, having the densities $s(t,x)$, $i(t,x)$ and $r(t,x)$ and a domain $\Omega \subset {\mathbb{R}}^n$ in which an epidemic is spread, one can calculate the numbers of each subpopulation using the formulae

$$S\left(t\right)={\int }_{\Omega }s(t,x)dx,I\left(t\right)={\int }_{\Omega }i(t,x)dx,R\left(t\right)={\int }_{\Omega }r(t,x)dx.$$

(14)

Because the three-component nonlinear reaction–diffusion system is a complicated object, typically, the systems involving only the first two equations are under study. In fact, assuming that the conservation law (2) is still valid, the function $R\left(t\right)$ can be easily found using (14). Pioneering works in which extensive research has been conducted to explore the dynamics of traveling waves of the two-component epidemic models were published in the 1970s–1980s [35,36,37]. In [35], the system (with $d_s=d_i$ and $n=1,2$)

$$\begin{array}{c}\frac{\partial s}{\partial t}=d_s\Delta s-\alpha si,\hfill \\ \frac{\partial i}{\partial t}=d_i\Delta i+\alpha si-\beta i,\hfill \end{array}$$

(15)

was suggested in order to describe the spread of the well-known black death pandemic. In particular, the conditions for the existence of traveling waves (nowadays, the terminology `traveling fronts’ is used) solutions are analyzed, and the numerical solution in the one-dimensional case $n=1$ is presented.

In [36], the authors consider system (15) with $d_s=0$ and $n=1$ for the spatial spread of rabies. In [37], a comprehensive analysis of this model is provided, including a proof of the existence of traveling waves and the conditions under which the waves travel at a minimal speed.

It should be pointed out that (15) is a particular case of the diffusive Lotka–Voltera system (see, e.g., [4,62])

$$\begin{array}{c}\frac{\partial u}{\partial t}=d_u\Delta u+u(a_1+b_{1}u+c_{1}v),\hfill \\ \frac{\partial v}{\partial t}=d_v\Delta v+v(a_2+b_{2}u+c_{2}v),\hfill \end{array}$$

(16)

where $u=u(t,x)$ and $v=v(t,x)$ are to-be-found functions, which usually represent densities, while $a_k,b_k$ and $c_k(k=1,2)$ are given parameters. Depending on the signs of the above parameters, system (16) describes several types of interactions between species, cells, chemicals, etc. In particular, the prey–predator model

$$\begin{array}{c}\frac{\partial u}{\partial t}=d_u\Delta u+u(a_1-c_{1}v),\hfill \\ \frac{\partial v}{\partial t}=d_v\Delta v+v(-a_2+b_{2}u)\hfill \end{array}$$

(17)

is obtained (here, $a_1,c_1,a_2$ and $b_2$ are positive constants). Obviously, system (17) with $a_1=0$ coincides with (16) up to notations.

It should be stressed that the construction of traveling fronts in explicit forms of the prey–predator system (17) is a highly nontrivial problem. To the best of our knowledge, the first examples of such types of solutions were presented in the recent review [63]. Therefore, it is not surprising that in the above-cited works devoted to the epidemic model (15), traveling fronts were not found.

The generalization of system (13) in the case when the parameters $\alpha$ and $\beta$ are assumed to be positive periodic continuous functions in t is considered in paper [38]. In this paper, the existence of periodic traveling wave solutions of the form $(s,i,r)=\left(S(x+\nu t),I(x+\nu t),R(x+\nu t)\right)$ (here, $\nu$ is the wave speed) is also analyzed. The existence of a traveling wave essentially depends on the basic reproduction number $R_0$. It is defined by the formula $R_0=\frac{\alpha }{\beta }S(-\infty )$ in the case of constant parameters $\alpha$ and $\beta$ and by the formula $R_0=\frac{{\int }_0^T\alpha \left(t\right)dt}{{\int }_0^T\beta \left(t\right)dt}S(-\infty )$ in the case of nonconstant parameters. Note that in the case $R_0\le 1$ (see, Theorem 3.1 [38]), there are no traveling wave solutions. This is in agreement with the general statement that a reaction–diffusion system can describe an epidemic process and a relevant traveling wave exists if $R_0>1$. This means that the disease can propagate and sustain itself in the population. If $R_0\le 1$, then the disease will die out and there will be no traveling wave solutions in the system.

In [40], the main focus is on the problem of how human behavior can affect the COVID-19 spread using a diffusive SEIR epidemic model. The model takes into account contact rate functions, describing different behaviors and interactions among individuals, and it reads as

$$\begin{array}{c}\frac{\partial s}{\partial t}=d_s\frac{{\partial }^{2}s}{\partial x^2}+\mathrm{\Lambda }-{\alpha }_e\left(x\right)\left(a_1-b_1{m}_1\left(e\right)\right)se-{\alpha }_i\left(x\right)\left(a_2-b_2{m}_2\left(i\right)\right)si-{\mu }_{nd}s,\hfill \\ \frac{\partial e}{\partial t}=d_e\frac{{\partial }^{2}e}{\partial x^2}+{\alpha }_e\left(x\right)\left(a_1-b_1{m}_1\left(e\right)\right)se+{\alpha }_i\left(x\right)\left(a_2-b_2{m}_2\left(i\right)\right)si-\gamma e-{\mu }_{nd}e,\hfill \\ \frac{\partial i}{\partial t}=d_i\frac{{\partial }^{2}i}{\partial x^2}+\gamma e-\beta i-\mu i-{\mu }_{nd}i,\hfill \\ \frac{\partial r}{\partial t}=d_r\frac{{\partial }^{2}r}{\partial x^2}+\beta i-{\mu }_{nd}r,\hfill \end{array}$$

(18)

where $\mathrm{\Lambda }$ is the influx rate of susceptible individuals; ${\mu }_{nd}$ is the natural death rate of the human; $\mu$ is the death rate of infected individuals due to COVID-19; $a_1$ and $a_2$ are the direct contact rates of e and i; $b_1$ and $b_2$ are the associated largest reduced rates due to human behavior changes of e and i; the functions ${\alpha }_e\left(x\right)$ and ${\alpha }_i\left(x\right)$ are the direct transmission contribution rates of e and i, which are probabilities that measure the contribution of spatial heterogeneity into direct human-to-human transmission; ${\alpha }_e\left(x\right)$ and ${\alpha }_i\left(x\right)$ are assumed to be non-negative and Hölder continuous functions, e.g., periodic trigonometric functions; and $m_1\left(e\right)$ and $m_2\left(i\right)$ are saturation functions that satisfy the restrictions

$$a_j\ge b_j>0,0\le m_j\le \frac{a_j}{b_j},m_j\in C^1\left([0,\infty )\right),m_j^{\prime }\ge 0,j=1,2.$$

In [40], the basic reproduction number $R_0$ is derived, and a threshold-type result on its global dynamics in terms of $R_0$ is established using the diffusive SEIR model (18). In order to define $R_0$ for the model in question, the authors analyze a linear system near the disease-free steady-state point $E_f$. Obviously, the nonlinear system (18) possesses such a point of the form $E_f=(s^{*},0,0,0)$ with $s^{*}=\frac{\mathrm{\Lambda }}{{\mu }_{nd}}$. The linearized system (18) in the vicinity of the equilibrium point $E_f$ reads as

$$\begin{array}{c}\frac{\partial s}{\partial t}=d_s\frac{{\partial }^{2}s}{\partial x^2}-a_1{\alpha }_e\left(x\right)s^{*}e-a_2{\alpha }_i\left(x\right)s^{*}i-{\mu }_{nd}{s}^{*},\hfill \\ \frac{\partial e}{\partial t}=d_e\frac{{\partial }^{2}e}{\partial x^2}+a_1{\alpha }_e\left(x\right)s^{*}e+a_2{\alpha }_i\left(x\right)s^{*}i-(\gamma +{\mu }_{nd})e,\hfill \\ \frac{\partial i}{\partial t}=d_i\frac{{\partial }^{2}i}{\partial x^2}+\gamma e-(\beta +\mu +{\mu }_{nd})i,\hfill \\ \frac{\partial r}{\partial t}=d_r\frac{{\partial }^{2}r}{\partial x^2}+\beta i-{\mu }_{nd}r.\hfill \end{array}$$

(19)

Since the equations for e and i do not involve s and r, the authors consider the following subsystem:

$$\begin{array}{c}\frac{\partial e}{\partial t}=d_e\frac{{\partial }^{2}e}{\partial x^2}+a_1{\alpha }_e\left(x\right)s^{*}e+a_2{\alpha }_i\left(x\right)s^{*}i-(\gamma +{\mu }_{nd})e,\hfill \\ \frac{\partial i}{\partial t}=d_i\frac{{\partial }^{2}i}{\partial x^2}+\gamma e-(\beta +\mu +{\mu }_{nd})i,\hfill \end{array}$$

(20)

and derive the formula for the basic reproduction number

$$R_0=sup\left\{\right|\lambda |:\lambda \in \sigma \left(L\right)\}.$$

(21)

In (21), $\sigma \left(L\right)$ is the spectral set of operators L defined as

$$L\left(\varphi \right)\left(x\right)=F\left(x\right){\int }_0^{\infty }T\left(t\right)\varphi \left(x\right)dt.$$

(22)

Here, $\varphi =({\varphi }_2,{\varphi }_3)\in C\left(\overline{\Omega },{\mathbb{R}}^2\right)$ represents the initial distribution of the densities e and i, while the function $F\left(x\right)$ and the operator $T\left(t\right)$ in (22) are defined by using the linear system (20).

Furthermore, it is proved that the disease-free state $E_f$ is stable if $R_0\le 1$, meaning the infection is not sustained in the population. However, if $R_0>1$, then a positive stationary solution exists, indicating that the epidemic can persist and spread throughout the population. To investigate further the impact of human behavior, the authors conduct numerical simulations based on their analytical findings. These simulations demonstrate that changes in human behavior can have a positive effect on reducing the infection level by decreasing the number of infected persons. Overall, the study emphasizes the importance of considering human behavior in modeling the spread of COVID-19 and highlights the potential effectiveness of behavior changes in mitigating the infection levels and reducing the overall impact of the pandemic.

In [64], a reaction–diffusion model

$$\begin{array}{c}\frac{\partial s}{\partial t}=d_s\frac{{\partial }^{2}s}{\partial x^2}+{\mu }_{nd}-{\alpha }_i\frac{si}{1+q{i}^2}-{\mu }_{nd}s,\hfill \\ \frac{\partial e}{\partial t}=d_e\frac{{\partial }^{2}e}{\partial x^2}+{\alpha }_i\frac{si}{1+q{i}^2}-(\gamma +{\mu }_e+ϵ+{\mu }_{nd})e,\hfill \\ \frac{\partial i}{\partial t}=d_i\frac{{\partial }^{2}i}{\partial x^2}+\gamma e-(\sigma +{\mu }_i+{\mu }_{nd})i,\hfill \\ \frac{\partial r}{\partial t}=d_r\frac{{\partial }^{2}r}{\partial x^2}+\sigma i+ϵe-{\mu }_{nd}r\hfill \end{array}$$

(23)

for the spread of COVID-19 is investigated. The model is a spatial extension of the SEIR model with nonlinear incidence rates by taking into account the effects of random movements of individuals from different compartments (subpopulations). The diffusive SEIR model (23) can be considered as a very particular case of model (18). Indeed, the authors use the specific function $\frac{1}{1+q{i}^2}$ instead of a general function $m_2\left(i\right)$ (see (18)). However, it can be noted that system (23) contains the new term $ϵe$ that takes into account the COVID-19 immunity of exposed individuals.

In system (23), q is the bilinear incidence rate; $ϵ$ is the immunity rate of exposed individuals; $\sigma$ is the rate of vaccination, quarantine or treatment; ${\mu }_e$ is the mortality rate of exposed individuals due to virus; and ${\mu }_i$ is the death rate of infected individuals due to virus. All other parameters have the same interpretations as above.

The results presented in [64] can be briefly summarized as follows: the equilibrium points and their stability are investigated; stability regions and the influence of key parameters are explored; a structure-preserving finite difference method for simulating the model in question is developed; consistency and stability analysis is provided; and the positivity of solutions is discussed.

In [47], a model for the spread of the COVID-19 epidemic is constructed based on the diffusive SEIR model by incorporating so-called asymptomatic infections. This model consists of three reaction–diffusion equations and two ODEs (see Equation (1) [47])

$$\begin{array}{c}\frac{\partial s}{\partial t}=d_s\left(t\right)▵s-f\left(t\right)\left({\alpha }_{e}e+{\alpha }_{i_s}{i}_s+{\alpha }_{i_a}{i}_a\right)\frac{s}{N},\hfill \\ \frac{\partial e}{\partial t}=d_s\left(t\right)▵e+f\left(t\right)\left({\alpha }_{e}e+{\alpha }_{i_s}{i}_s+{\alpha }_{i_a}{i}_a\right)\frac{s}{N}-\gamma e,\hfill \\ \frac{\partial i_a}{\partial t}=d_s\left(t\right)▵i_a+(1-p)\gamma e-\beta i_a,\hfill \\ \frac{\partial i_s}{\partial t}=p\gamma e-\beta i_s-\mu i_s,\hfill \\ \frac{\partial r}{\partial t}=\beta (i_a+i_s).\hfill \end{array}$$

(24)

In system (24), the functions $i_a$ and $i_s$ are the densities of asymptomatic and symptomatic infected individuals, respectively; the diffusion $d_s\left(t\right)$ has one of the forms presented in Formula (3) [47]; the function $f\left(t\right)$ represents the average number of contacts (see Equation (2) in [47]); the parameter p is the probability of being confirmed, while $(1-p)$ is the probability of being unreported; other parameters have the same interpretations as above.

Using official data, spatial modeling of the density distribution of symptomatic infected individuals in France is performed during the first wave of the pandemic (from 24 January to 16 June 2020). The computational results show good agreement with the official data. It is demonstrated that the total number of cases would be significantly higher without intervention of the government (i.e., without implementing a lockdown).

In [45,46], reaction–diffusion models are proposed to investigate the impact of vaccination and isolation strategies on the progression of the epidemic. In [45], the global asymptotic stability and the persistence of the epidemic are proven using a reaction–diffusion model for the HBV epidemic. Results of numerical simulations are presented as well.

The study [46] begins by exploring the fundamental dynamic properties of the diffusive epidemic system. Subsequently, the asymptotic distributions of the endemic equilibrium under different conditions are analyzed. Overall, this research contributes to the ongoing efforts in epidemic prevention and control by providing insights into the dynamics of COVID-19 and suggesting optimal vaccination and isolation strategies. The relevant model of the COVID-19 epidemic was constructed taking into account the vaccination of patients. The model is based on the system of four reaction–diffusion equations for a heterogeneous medium and has the form

$$\begin{array}{c}\frac{\partial s}{\partial t}=d_s\Delta s+f_1\left(x\right)-{\alpha }_1\left(x\right)\frac{k}{1+a_{1}ki}si-[\sigma \left(x\right)+{\mu }_{nd}\left(x\right)]s+\omega \left(x\right)r,\hfill \\ \frac{\partial v}{\partial t}=d_v\Delta v+f_2\left(x\right)-{\alpha }_2\left(x\right)\frac{k}{1+a_{2}ki}vi+\sigma \left(x\right)s-{\mu }_{nd}\left(x\right)v,\hfill \\ \frac{\partial i}{\partial t}=d_i\Delta i+{\alpha }_1\left(x\right)\frac{k}{1+a_{1}ki}si+{\alpha }_2\left(x\right)\frac{k}{1+a_{2}ki}vi-[\beta \left(x\right)+\mu \left(x\right)+{\mu }_{nd}\left(x\right)]i,\hfill \\ \frac{\partial r}{\partial t}=d_r\Delta r+\beta \left(x\right)i-[\omega \left(x\right)+{\mu }_{nd}\left(x\right)]r,\hfill \end{array}$$

(25)

where the densities $s(t,x)$, $v(t,x)$, $i(t,x)$ and $r(t,x)$ stand for unvaccinated susceptible individuals, vaccinated susceptible individuals, infected individuals, and recovered individuals, respectively; $k=1-m,m\in [0,1)$ (here, m means the fraction of the infected subpopulation). It is assumed that all parameters and functions in (25) are non-negative. The parameters $f_1\left(x\right)$ and $f_2\left(x\right)$ are the inputs to s and v, respectively; $a_1$ and $a_2$ are so-called half-saturation parameters; $\sigma \left(x\right)$ denotes the vaccination rate of s; $\omega \left(x\right)$ indicates the immune loss rate of r; other parameters have the same interpretations as above.

To construct a numerical solution of the reaction–diffusion model (25), the finite difference method in time and space is used. Plots of the numerical solutions are presented in the 1D space approximation. The authors pay special attention to the cases when one or more diffusivities are small. For numerical simulations, the parameters and functions in (25) were specified using available experimental data from a wide range of references. Based on the obtained numerical results, the sensitivity index of the relevant parameters for the basic reproduction number $R_0$ is obtained. It is shown that the functions $f_1\left(x\right)$, ${\beta }_1\left(x\right)$, ${\alpha }_2\left(x\right)$ and $\beta \left(x\right)$ are highly sensitive parameters, while the functions ${\mu }_{nd}\left(x\right),\omega \left(x\right)$ and $\mu \left(x\right)$ are insensitive parameters. The sensitivity of the parameters provides technical guidance for COVID-19 prevention and control. Controlling the high sensitivity parameter, one can better reduce the basic reproduction number $R_0$.

It should be pointed out that all the models presented above involve constant (or time-dependent) diffusivities. Recently, several papers were published [41,42,43,44] by an international group of researchers in which a new SEIRD-type model with nonlinear diffusion terms was introduced and examined by analytical and numerical methods. In the simplest case, the model has the form [42,43]

$$\begin{array}{c}\frac{\partial s}{\partial t}=\nabla ·\left(d_{s}N\nabla e\right)-\frac{\alpha }{N}si,\hfill \\ \frac{\partial e}{\partial t}=\nabla ·\left(d_{e}N\nabla e\right)+\frac{\alpha }{N}si-\gamma e,\hfill \\ \frac{\partial i}{\partial t}=\nabla ·\left(d_{i}N\nabla i\right)+\gamma e-\beta i-\mu i,\hfill \\ \frac{\partial r}{\partial t}=\nabla ·\left(d_{r}N\nabla r\right)+\beta i,\hfill \\ \frac{\partial d}{\partial t}=\mu i.\hfill \end{array}$$

(26)

Hereinafter, $N(t,x)$ denotes the sum of the living population densities, i.e., $N=s+e+i+r$; the diffusion parameters $d_s,d_e,d_i$ and $d_r$ may depend on time and space. Here, all diffusivity coefficients are proportional to the total population and can be locally adjusted to incorporate geographical or human-related inhomogeneities (see Section 3 in [65] for details).

In [43], a further generalization of the above model is suggested:

$$\begin{array}{c}\frac{\partial s}{\partial t}=\nabla ·\left(d_{s}N\nabla s\right)-\left(1-\frac{A}{N}\right){\alpha }_{i}si-\left(1-\frac{A}{N}\right){\alpha }_{e}se,\hfill \\ \frac{\partial e}{\partial t}=\nabla ·\left(d_{e}N\nabla e\right)+\left(1-\frac{A}{N}\right){\alpha }_{i}si+\left(1-\frac{A}{N}\right){\alpha }_{e}se-\gamma e-{\beta }_{e}e,\hfill \\ \frac{\partial i}{\partial t}=\nabla ·\left(d_{i}N\nabla i\right)+\gamma e-{\beta }_{i}i-\mu i,\hfill \\ \frac{\partial r}{\partial t}=\nabla ·\left(d_{r}N\nabla r\right)+{\beta }_{i}i+{\beta }_{e}e,\hfill \\ \frac{\partial d}{\partial t}=\mu i,\hfill \end{array}$$

(27)

where a new parameter A is used to describe the Allee effect (depensation), which serves for describing a tendency of outbreaks to cluster toward small population centers. The Allee effect is a phenomenon in population dynamics that describes a decrease in the per capita growth rate of a population as it approaches lower population densities. In the context of disease outbreaks, the Allee effect is often used to model the tendency for outbreaks to occur more frequently and cluster in smaller population centers or areas with lower population densities.

The diffusive SEIRD system (27) is solved using finite element methods [43]. To verify the model accuracy, the obtained results are cross-referenced with the research outcomes of other authors.

In [41,42], a modification of the diffusive SEIRD system (27) by incorporating the general (nondisease) mortality rates is studied. In the modified model, the linear terms $-{\mu }_{nd}s,-{\mu }_{nd}e,-{\mu }_{nd}i$ and $-{\mu }_{nd}r$ are present on the right-hand side of the first four equations of system (27), respectively.

In [41], the diffusive SEIRD model is used to describe the spatiotemporal spread of the COVID-19 pandemic and aims to capture dynamics based on human behaviors and geographical features. To validate the model, numerical results are compared with measurement data from the Lombardy region in Italy, which was severely affected by the crisis between February and April 2020. The obtained results demonstrate qualitative agreement between the modeled spatiotemporal spread of COVID-19 in Lombardy and epidemiological data. It is concluded that the numerical results can be used to inform healthcare authorities in developing effective measures to mitigate the pandemic and anticipate the geographical distribution of critical medical resources.

In [42], the authors analized an ODE version of the diffusive SEIRD model to derive a basic reproduction number $R_0$. Additionally, the authors explored the role of diffusion and $R_0$ in shaping the behavior of the diffusive SEIRD model. Through numerical simulations, they investigated how the interplay between these factors affects the dynamics of the epidemic.

In order to show the applicability of the model, the role of diffusion is demonstrated in the case of the Lombardy region (Italy). It is shown that the mathematical model (26) reproduces the COVID-19 epidemic spread in Lombardy, beginning on 27 February 2020. In addition to earlier simulations presented in [41], two additional cases are investigated in [42]. In the first case, the values of $d_s$, $d_e$, $d_i$, and $d_r$ are doubled, while in the second case, they are halved. Furthermore, a scenario is examined in which $d_r$, $d_e$, and $d_i$ are doubled, but $d_s$ is halved. This choice closely resembles the parameter configuration used in the 2D simulations. The primary objective is to prevent the potential occurrence of nonphysical diffusion within the susceptible population, which may result in a general decrease of the population density. It is also observed that a wider geographic range of affected areas is produced by larger diffusion, which is in agreement with the physical meaning of diffusion. This effect is particularly evident in the southeastern clusters (see Figure 11 in [42]). In the case of double diffusion, a homogeneous and continuous region of infection is generated. In contrast, more localized dynamics are observed in the case of half diffusion, resulting in a clear separation into distinct regions. The case in which $d_r$, $d_e$, and $d_i$ are doubled and $d_s$ is halved produces intermediate results between the double-diffusion and half-diffusion cases.

4. Other Epidemic Models Taking into Account Spatial Heterogeneity

There are several models for simulating epidemic spreads in time and space that cannot be treated as direct generalizations of the classical models presented in Section 2 by adding diffusion terms. Some of them are presented in this section. It should be stressed that almost all of them were developed after the COVID-19 outbreak.

We start from the simplest model whose type was introduced in [39]. The model reads as

$$\begin{array}{c}\frac{\partial s}{dt}=-\alpha si,\hfill \\ \frac{\partial i}{dt}+\nabla ·\left(iV\right)=\alpha si-(\beta +\mu )i,\hfill \\ \frac{\partial r}{dt}=\beta i,\hfill \end{array}$$

(28)

where the functions $s,i,$ and $r$ and the parameters $\alpha ,\beta$ and $\mu$ have the same interpretation as in the previous section. Spatial movement for the infected population relative to the medium is investigated by introducing the function $V(t,x)$ satisfying the Euler equation

$$i\frac{\partial V}{\partial t}+iV·\nabla V=-\nabla p,$$

(29)

where p is the pressure, which is an internal force of the fluid and assumed to be a known smooth function. It should be noted that the model (28)–(29) consists of the first-order equations in contrast to the models presented in Section 3. In particular, it means that much simpler equations are obtained for searching for traveling waves. Moreover, it can be noted that an autonomous system for finding the functions i and V consisting of the first-order PDEs (29) and

$$\nabla ·\left(iV\right)=-\mu i$$

is obtained providing the conservation law (2) is preserved for the densities $s,i,$ and $r$.

In models (28)–(29), the motion of the infected population, which characterizes the spatial spread of the epidemic, is described as an inviscid fluid. It is important to note that there is no requirement to assume the susceptible population remains stationary in this modeling framework. Instead of explicitly focusing on the movement of susceptible individuals, the susceptible population is treated as a medium, and the motion of the infected fluid is analyzed in relation to this medium. While individuals may exhibit diverse and random movements at the individual level, this model is based on the underlying assumption that the spread of an epidemic can be approximated by an inviscid flow at the macroscopic level. As an illustrative demonstration of model applicability, the spread of the COVID-19 epidemic within Wuhan, China 2020, was used.

To find the numerical solution of model (28)–(29), the finite difference method was applied. The difference scheme obtained was tested using a known analytical solution. Having completed this, the numerical simulations were performed taking into account the number of residents and the area of Wuhan. The authors modeled the spread of the epidemic over a square area equal to the area of the city, assuming that the source of infection is located in the central part of the city. The results obtained by numerical simulations are in good agreement with the official data of the COVID-19 disease in Wuhan.

In paper [66], a spatiotemporal model of the spread of the COVID-19 epidemic with moving boundaries was constructed. The governing equations are direct generalization of those presented in model (28)–(29). As a result, a system of four equations of the SEIR model with convective terms was obtained, where the speed of movement of individuals is taken into account. The model reads as

$$\begin{array}{c}\frac{\partial s}{dt}+\nabla ·\left(sV\right)=\mathrm{\Gamma }-{\alpha }_{e}se-{\alpha }_{i}si,\hfill \\ \frac{\partial e}{dt}+\nabla ·\left(eV\right)={\alpha }_{e}se+{\alpha }_{i}si-\gamma e,\hfill \\ \frac{\partial i}{dt}+\nabla ·\left(iV\right)=\gamma e-(\mu +\beta )i,\hfill \\ \frac{\partial r}{dt}+\nabla ·\left(rV\right)=\beta i,\hfill \end{array}$$

where the function V is the speed that characterizes the epidemic flow; the function $\mathrm{\Gamma }$ represents the rate of the density change for the susceptible individuals due to the expansion of the epidemic domain; all other parameters have the same interpretations as above. The functions V and $\mathrm{\Gamma }$ generally depend on both the time and space, while all other parameters are assumed to be constants.

The authors assume that there are no infected individuals outside a 2D space domain with the boundary Z. Because one needs to introduce some assumptions about the moving boundary Z, the authors studied a radially symmetric case. In this case, the model essentially simplifies, in particular, that the boundary Z is nothing else but a circle with the time-dependent radius L. Thus, using rather a standard approach, the equations for L and the speed V were derived (see (2.6) and (2.7) in [66]). The model is verified on official data on the epidemic in Wuhan for the time frame: from the 23th January to the 18th February, 2020. Based on the results of calculations, the radial distribution of speed and infected individuals are found.

A model involving cross-diffusion was studied in [51]. The model is a generalization of the ODE model (see (9) above)) developed in [15] and reads as

$$\begin{array}{c}\frac{\partial u}{\partial t}=d_1\Delta u+u(a-b{u}^{\kappa }),\hfill \\ \frac{\partial v}{\partial t}=d_2\Delta u+k\left(t\right)u.\hfill \end{array}$$

(30)

In (30), the function $u(t,x,y)$ describes the density (rate) of the infected persons (the number of COVID-19 cases) in a vicinity of the point $(x,y)$, while $v(t,x,y)$ means the density of the deaths from COVID-19. The diffusivity coefficients $d_1$ and $d_2$ describe the random movement of the infected persons, which lead to increasing the pandemic spread. Each coefficient in the reactions terms, $a,b,\kappa$ and $k\left(t\right)$, has the clear meaning described above for the ODE model (9).

It turns out that a wide range of exact solutions (including traveling fronts) of the nonlinear system (30) can be constructed using the Lie symmetry analysis. We have also demonstrated that the exact solutions obtained are useful for describing the spread of the COVID-19 pandemic in 1D approximation.

Let us present some details. Taking into account the assumptions that the distribution of the infected persons is one-dimensional in space (i.e., the diffusion with respect to the axis y is very small) and that $d_1\gg d_2$ (i.e., the space diffusion of the infected persons leads mostly to increasing the total number of the COVID-19 cases and not so much to new deaths), system (30) with $k\left(t\right)=k_0{e}^{-\alpha t}$ (hereafter $k_0>0,\alpha >0$) takes the form

$$\begin{array}{c}\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}+u\left(1-u^{\kappa }\right),\hfill \\ \frac{\partial v}{\partial t}=\frac{k_0}{a}exp\left(-\frac{\alpha t}{a}\right)u.\hfill \end{array}$$

(31)

System (31) admits the Lie symmetry operator (see case 4 in Theorem 2.3 [51])

$$X=c{\partial }_x+{\partial }_t-\frac{\alpha }{a}v{\partial }_v,c\in \mathbb{R},$$

that allowed us to construct its exact solution in the following form

$$\begin{array}{c}u(t,x)={\left(1+Aexp\left(\frac{\kappa }{\sqrt{2(\kappa +2)}}\omega \right)\right)}^{-2/\kappa },\omega =\pm x-\frac{\kappa +4}{\sqrt{2(\kappa +2)}}t,A>0,\hfill \\ v(t,x)=\frac{k_0}{a}\int exp\left(-\frac{\alpha t}{a}\right){\left(1+Aexp\left(\frac{\kappa }{\sqrt{2(\kappa +2)}}\omega \right)\right)}^{-2/\kappa }dt+g\left(x\right),\hfill \end{array}$$

(32)

where $g\left(x\right)$ is an arbitrary smooth function.

Remark 1. 

The assumption $d_1\gg d_2$ is not essential, and one can construct exact solutions of the form (32) for (30) in one-dimensional approximation as well. So, setting $d_1=1,d_2=D>0$, the exact solution

$$\begin{array}{c}u(t,x)={\left(1+F\right)}^{-2/\kappa },\hfill \\ v(t,x)=\frac{k_0}{a}\int exp\left(-\frac{\alpha t}{a}\right){\left(1+F\right)}^{-2/\kappa }dt-\frac{D}{2+\kappa }\int F\left(\kappa -2F\right){\left(1+F\right)}^{-2-2/\kappa }dt+g\left(x\right),\hfill \end{array}$$

can be found. Here, $F=Aexp\left(\frac{\kappa }{\sqrt{2(\kappa +2)}}\omega \right),\omega =\pm x-\frac{\kappa +4}{\sqrt{2(\kappa +2)}}t,A>0$.

The integral in (32) can be expressed in the terms of elementary functions in some particular cases. Taking $\alpha =\frac{5}{6}a$ and $\kappa =1$ , for example, one can obtain the exact solution

$$\begin{array}{c}u(t,x)={\left(1+Aexp\left(\frac{1}{\sqrt{6}}\omega \right)\right)}^{-2},\omega =x-\frac{5}{\sqrt{6}}t,\hfill \\ v(t,x)=g\left(x\right)-\frac{6k_0}{5a}exp\left(-\frac{x}{\sqrt{6}}\right){\left(A+exp\left(-\frac{1}{\sqrt{6}}\omega \right)\right)}^{-1},\hfill \end{array}$$

(33)

of system (31).

Note that the functions u and v in (33) should be non-negative for any $t>0$ and $x\in \mathbb{I}$ (here, $\mathbb{I}\subset \mathbb{R}$) because they represent the densities. Obviously, the function u is always positive. It is easily seen that each function $g\left(x\right)$ satisfying the inequality

$$g\left(x\right)\ge \frac{6k_0}{5a}exp\left(-\frac{x}{\sqrt{6}}\right){\left(A+exp\left(-\frac{x}{\sqrt{6}}\right)\right)}^{-1}$$

guarantees also the non-negativity of v. In particular, one may take the function

$$g\left(x\right)=\frac{6k_0}{5a}exp\left(-\frac{x}{\sqrt{6}}\right){\left(A+exp\left(-\frac{x}{\sqrt{6}}\right)\right)}^{-1},$$

which guarantees the zero density of the deaths in the initial time $t=0$, i.e., $v(0,x)=0$.

Examining the space interval $\mathbb{I}=[x_1,x_2],x_1<x_2$, one can calculate the total number of COVID-19 cases and deaths on this interval as follows

$$\begin{array}{c}U\left(t\right)={\int }_{x_1}^{x_2}u(t,x)dx,\hfill \\ V\left(t\right)={\int }_{x_1}^{x_2}v(t,x)dx.\hfill \end{array}$$

(34)

So, substituting solution (33) into (34), we arrive at the formulae

$$\begin{array}{c}U\left(t\right)=(x_2-x_1)-\sqrt{6}[{\left(1+Aexp\left(\frac{x_1-\frac{5}{\sqrt{6}}t}{\sqrt{6}}\right)\right)}^{-1}-{\left(1+Aexp\left(\frac{x_2-\frac{5}{\sqrt{6}}t}{\sqrt{6}}\right)\right)}^{-1}\hfill \\ +ln\left(1+Aexp\left(\frac{x_2-\frac{5}{\sqrt{6}}t}{\sqrt{6}}\right)\right)-ln\left(1+Aexp\left(\frac{x_1-\frac{5}{\sqrt{6}}t}{\sqrt{6}}\right)\right)],\hfill \\ V\left(t\right)={\int }_{x_1}^{x_2}g\left(x\right)dx-\frac{6\sqrt{6}{k}_0}{5a}{e}^{-\frac{5}{6}t}[ln\left(1+Aexp\left(\frac{\frac{5}{\sqrt{6}}t-x_1}{\sqrt{6}}\right)\right)-ln\left(1+Aexp\left(\frac{\frac{5}{\sqrt{6}}t-x_2}{\sqrt{6}}\right)\right].\hfill \end{array}$$

Obviously, the functions $U\left(t\right)$ and $V\left(t\right)$ are increasing and bounded, because

$$(U,V)\to \left((x_2-x_1),{\int }_{x_1}^{x_2}g\left(x\right)dx\right)ast\to +\infty .$$

Moreover, taking the appropriate function $g\left(x\right)$, we can guarantee that

$$U\left(0\right)=U_0\ge 0,V\left(0\right)=V_0\ge 0.$$

Thus, one may claim that the exact solution (33) possesses all necessary properties for the description of the distribution of the COVID-19 cases (total number of infected population as well) and the deaths from this virus in time–space.

Interestingly, the spread of the COVID-19 cases in space has the form of a traveling wave, and this qualitatively coincides with the real situation in many countries during the first wave of the COVID-19 pandemic. In Ukraine, for example, the pandemic started in the western part of the country and then spread to the central and eastern parts of Ukraine (the major exception was only the capital Kyiv, in which the total number of COVID-19 cases was high from the very beginning).

It should be stressed that cross-diffusion is a distinguished peculiarity in some well-known mathematical models arising in biology, ecology and medicine [52,53] that were extensively studied by different mathematical techniques (see, e.g., [54] and references therein) during the last decade. So, we believe that relevant terms with cross-diffusion should naturally arise in epidemic models. To the best of our knowledge, the first epidemic model with cross-diffusion was briefly presented in [7] (see Section 9.2 therein) that reads as

$$\begin{array}{c}\frac{\partial s}{\partial t}=d_s\Delta s-d_{is}s\Delta i-\alpha si+b,\hfill \\ \frac{\partial i}{\partial t}=d_i\Delta i+d_{is}s\Delta i+\alpha si-\beta i,\hfill \end{array}$$

(35)

where $d_{is}>0$ is cross-diffusivity and the term $d_{is}s\Delta i$ reflects infection caused by isotropic movement in space (in contrast to the term $\alpha si$ reflecting a local infection). Here, the researchers also introduced a birth parameter b; however, we believe that this parameter should be a function of s, and the simplest possibility is $b=b_{1}s,b_1>0$. Moreover, a similar term, say $b=b_{2}i,b_2>0$, should be added in the second equation because the infected population can also produce new members. So, the above model with cross-diffusion (35) can be generalized as follows

$$\begin{array}{c}\frac{\partial s}{\partial t}=d_s\Delta s-d_{is}s\Delta i-\alpha si+b_{1}s,\hfill \\ \frac{\partial i}{\partial t}=d_i\Delta i+d_{is}s\Delta i+\alpha si+(b_2-\beta )i.\hfill \end{array}$$

(36)

It can be noted that the nonlinear model (36) can be essentially simplified in the special case $d_s=d_i$ (a plausible assumption) and $b_1=b_2-\beta$ when one of the equations in (36) can be replaced by the linear reaction–diffusion equation

$$\frac{\partial z}{\partial t}=d_s\Delta z+b_{1}z,z=s+i.$$

5. Age-Structured Epidemic Models

In this section, we present some information about so-called age-structured epidemic models. The simplest representative of such types of models reads as

$$\begin{array}{c}\frac{\partial s}{\partial t}+\frac{\partial s}{\partial a}=-\lambda (a,t)s-{\mu }_{nd}\left(a\right)s,\hfill \\ \frac{\partial i}{\partial t}+\frac{\partial i}{\partial a}=\lambda (a,t)s-\beta \left(a\right)i-{\mu }_{nd}\left(a\right)i,\hfill \\ \frac{\partial r}{\partial t}+\frac{\partial r}{\partial a}=\beta \left(a\right)i-{\mu }_{nd}\left(a\right)r,\hfill \end{array}$$

(37)

where $s(a,t),i(a,t)$ and $r(a,t)$ are the age-specific density of susceptible, infective and recovered individuals of age a at time t, respectively; the function $\lambda (a,t)$ is the age-specific force of infection (the probability for a susceptible of age a to become infective in a unit time interval); ${\mu }_{nd}\left(a\right)$ and $\beta \left(a\right)$ are the age-specific death and recovery rates, respectively.

Depending on the form of the function $\lambda (a,t)$, system (37) can be a system of the two-dimensional first-order PDEs or a system of integro-differential equations. System (37) with

$$\lambda (a,t)=k\left(a\right)i(a,t)\mathtt{and}\lambda (a,t)=k\left(a\right){\int }_0^{\infty }i(a,t)da$$

(here, $k\left(a\right)$ is non-negative bounded continuous function on $[0,\infty )$) was suggested in [67,68]. In particular, an endemic threshold criteria is derived and the stability of steady-state solutions of system (37) is determined therein.

System (37) with a more general form of the function $\lambda (a,t)$ on the the interval $[0,a_{max}]$ (one is more realistic) is studied in [69]. Setting

$$\lambda (a,t)={\int }_0^{a_{max}}k(a,\sigma )i(\sigma ,t)d\sigma$$

(here, $k(a,\sigma )$ is the age-dependent transmission coefficient, i.e., the probability that a susceptible person of age a meets an infectious person of age $\sigma$ and becomes infected, per unit of time, $a_{max}$ is a maximal age) conditions that guarantee the existence and uniqueness of nontrivial steady states of the age-structured model are derived, and the local and global stabilities of the steady states are examined.

The age-structured SIR model (37) with $\lambda \left(t\right)={\int }_0^{\infty }k\left(a\right)i(a,t)da$ is studied in [6] (see Section 6.4 therein). In particular, the basic reproduction number $R_0$ is obtained (see formula 6.72 therein). A limiting case of the age-structured SIR model (37) is called the age-structured SIS model, and that was studied in detail. The model reads as

$$\begin{array}{c}\frac{\partial s}{\partial t}+\frac{\partial s}{\partial a}=-\lambda \frac{\alpha \left(a\right)}{N}s+\beta \left(a\right)i-{\mu }_{nd}\left(a\right)s,\hfill \\ \frac{\partial i}{\partial t}+\frac{\partial i}{\partial a}=\lambda \frac{\alpha \left(a\right)}{N}s-\beta \left(a\right)i-{\mu }_{nd}\left(a\right)i,\hfill \\ s(0,t)={\int }_0^{\infty }b\left(a\right)\left(s+i\right)da,i(0,t)=0,\hfill \\ \lambda ={\int }_0^{\infty }k\left(a\right)i(a,t)da,N={\int }_0^{\infty }\left(s+i\right)da.\hfill \end{array}$$

(38)

Here, $N\left(t\right)$ is the total population, and $b\left(a\right)$ is the birth rate. The exact solution of the age-structured SIS model (38) is constructed in the form [6]

$$s(a,t)=S\left(a\right)e^{\kappa t},i(a,t)=I\left(a\right)e^{\kappa t},$$

(39)

where $\kappa$ is a to-be-determined constant. Using ansatz (39), the age-structured SIS model (38) is reduced to the following form

$$\begin{array}{c}\frac{dS}{da}=-{\lambda }^{*}\frac{\alpha \left(a\right)}{N^{*}}S+\beta \left(a\right)I-{\mu }_{nd}\left(a\right)S-\kappa S,\hfill \\ \frac{dI}{da}={\lambda }^{*}\frac{\alpha \left(a\right)}{N^{*}}S-\beta \left(a\right)I-{\mu }_{nd}\left(a\right)I-\kappa I,\hfill \\ S\left(0\right)={\int }_0^{\infty }b\left(a\right)\left(S+I\right)da,I\left(0\right)=0,\hfill \\ {\lambda }^{*}={\int }_0^{\infty }k\left(a\right)I\left(a\right)da,N^{*}={\int }_0^{\infty }\left(S+I\right)da.\hfill \end{array}$$

(40)

Formulae (40) form a boundary-value problem with the governing equations in the form of ordinary integral–differential equations. Using the governing equations from (40), one arrives at the relation

$$S\left(a\right)+I\left(a\right)=Cexp\left(-\kappa a-{\int }_0^a{\mu }_{nd}\left(\sigma \right)d\sigma \right).$$

(41)

Assuming $S\left(0\right)+I\left(0\right)=1$ (without losing a generality), one may set $C=1$. Thus, using (41) and the boundary conditions from (40), the parameter $\kappa$ is defined by the transcendent equation

$${\int }_0^{\infty }b\left(a\right)exp\left(-\kappa a-{\int }_0^a{\mu }_{nd}\left(\sigma \right)d\sigma \right)da=1.$$

(42)

Integrating the second equation from (40) and taking into account the above formulae, the exact solution of the boundary-value problem (40) was found [6]

$$\begin{array}{c}S\left(a\right)=exp\left(-\kappa a-{\int }_0^a{\mu }_{nd}\left(\sigma \right)d\sigma \right)-I\left(a\right),\hfill \\ I\left(a\right)={\lambda }^{*}exp\left(-\kappa a-{\int }_0^a{\mu }_{nd}\left(\sigma \right)d\sigma \right){\int }_0^a\frac{\alpha \left(\sigma \right)}{N^{*}}exp\left(-{\int }_{\sigma }^a\beta \left(\varsigma \right)d\varsigma -{\lambda }^{*}{\int }_{\sigma }^a\frac{\alpha \left(\varsigma \right)}{N^{*}}d\varsigma \right)d\sigma ,\hfill \end{array}$$

(43)

where ${\lambda }^{*}$ satisfies the equation

$${\int }_0^{\infty }k\left(a\right)exp\left(-\kappa a-{\int }_0^a{\mu }_{nd}\left(z\right)dz\right){\int }_0^a\frac{\alpha \left(\sigma \right)}{N^{*}}exp\left(-{\int }_{\sigma }^a\beta \left(\varsigma \right)d\varsigma -{\lambda }^{*}{\int }_{\sigma }^a\frac{\alpha \left(\varsigma \right)}{N^{*}}d\varsigma \right)d\sigma da=1.$$

(44)

Substituting (43) into ansatz (39) and taking into account (42), (44) and formulae ${\lambda }^{*}=e^{-\kappa t}\lambda ,N^{*}=e^{-\kappa t}N$, one obtains the exact solution of the age-structured SIS model (38).

Obviously, age played an essential role during the COVID-19 pandemic because the deaths rate depends very much on the age of infected persons. Thus, the development and application of the age-structured epidemic models in the modeling of the COVID-19 pandemic are important, and several studies are devoted to these aspects [70,71,72]. In [70], the age-structured SEIR model

$$\begin{array}{c}\frac{\partial s}{\partial t}+\frac{\partial s}{\partial a}=-\alpha s\left({\int }_0^{+\infty }i(a,t)da+{\int }_0^{+\infty }e(a,t)da\right),\hfill \\ \frac{\partial e}{\partial t}+\frac{\partial e}{\partial a}=\alpha s\left({\int }_0^{+\infty }i(a,t)da+{\int }_0^{+\infty }e(a,t)da\right)-\gamma e,\hfill \\ \frac{\partial i}{\partial t}+\frac{\partial i}{\partial a}=\gamma e-\beta i-\mu \left(a\right)i,\hfill \\ \frac{dr}{dt}=\beta {\int }_0^{+\infty }i(a,t)da,\hfill \end{array}$$

(45)

is considered in order to predict the epidemic peak outbreak in South Africa, Turkey and Brazil. In (45), all parameters have the same interpretations as above. A generalization of the above model (see (6) [70]) is also suggested by adding a new subpopulation (compartment) for individuals in quarantine.

In [71], another approach was suggested by considering an age-structured model that takes into account two main components of the COVID-19 pandemic: the number of infected individuals requiring hospitalization, which leads to the estimation of required beds, and the potential infection of healthcare personnel. Consequently, the model predicts the timing of the peak and the number of infectious cases at that peak both before and after the implementation of nonpharmaceutical interventions. Additionally, a comparison is made with the scenario of a full lockdown. In [72], an age-structured model is proposed to study the outbreak of the COVID-19 coronavirus in Wuhan, China. The value of the basic reproduction number is computed in order to provide an initial understanding of how contagious or virulent the pandemic is and how one might spread within a population.

In order to take into account both age and spatial heterogeneity, more complicated models were developed [55,56,57,58]. In particular, the age-structured SIS epidemic model with diffusion [55]

$$\begin{array}{c}\frac{\partial s}{\partial t}+\frac{\partial s}{\partial a}=d_s\Delta s-\lambda s+\gamma i-\mu s,\hfill \\ \frac{\partial i}{\partial t}+\frac{\partial i}{\partial a}=d_i\Delta i+\lambda s-\gamma i-\mu i,\hfill \end{array}$$

(46)

was further developed. In (46), $s(a,t,x)$ and $i(a,t,x)$ are the densities of susceptible and infective individuals of age $a\ge 0$ at time t in position $x\in \Omega \subset {\mathbb{R}}^n$, respectively; $\mu (a,x)$ is the mortality of an individual of age a in position x; $\gamma (a,x)$ is the recovery rate of an infective individual of age a in position x; $\lambda (a,t,x)$ is the force of infection to a susceptible individual of age a at time t in position x given by formula

$$\lambda (a,t,x)={\int }_0^{a_{max}}{\int }_{\Omega }k(a,\sigma ,x,y)i(\sigma ,t,y)dyd\sigma ,$$

where $k(a,\sigma ,x,y)$ denotes the rate of disease transmission from an infective individual of age $\sigma$ in position y to a susceptible individual of age a in position x.

In [55,56], the existence of nontrivial steady-state solutions of the epidemic model (46) is investigated; the basic reproduction number $R_0$ for the system in question is estimated; numerical simulations are performed to verify the analytical results obtained.

In [57], the well-posedness of an age-structured SIS model is proved, the existence and uniqueness of the nontrivial steady state corresponding to an endemic state are investigated, and the local and global stability of this nontrivial steady state are studied. Furthermore, the asymptotic properties of the principal eigenvalue and the nontrivial steady state with respect to the nonlocal diffusion rate are discussed.

In [58], the modified age-structured SIS model

$$\begin{array}{c}\frac{\partial s}{\partial t}=d_s\Delta s-f\left(s\right){\int }_0^{a_{max}}k\left(a\right)g\left(i(a,t,x)\right)da,\hfill \\ \frac{\partial i}{\partial t}+\frac{\partial i}{\partial a}=d_i\Delta i-\gamma i-\mu i,\hfill \end{array}$$

is studied. It can be noted that the function s, describing susceptible individuals, does not depend on the variable a (in contrast to the function i) but only on the maximal population age $a_{max}$. Here, the smooth non-negative functions $f\left(s\right)$ and $g\left(i\right)$ satisfy the following conditions:

$$f\left(0\right)=g\left(0\right)=0,\frac{df}{ds}>0,\frac{dg}{di}\ge 0,\frac{d^{2}g}{d{i}^2}\le 0.$$

Moreover, the function $\frac{g\left(i\right)}{i}$ is continuously differentiable and nonincreasing. The existence of traveling fronts of the form $(s(a,t),i(a,t,x))=\left(S(x+\nu t),I(a,x+\nu t)\right)$ (here, $\nu$ is the wave speed) of the above model is proved under the following asymptotic boundary conditions

$$S(-\infty )=S_0,S(+\infty )=S_1,I(a,-\infty )=I(a,+\infty )=0,0\le S_1\le S_0<\infty .$$

The effects of nonlinear functions f and g and age structure on the basic reproduction number and critical wave speed are investigated in this paper as well.

However, the basic equations of the model, multidimensional nonlinear integro-differential equations with nonconstant coefficients, are very complicated. One needs to simplify the model in order to derive results that are important from an applicability point of view.

6. Conclusions

This review provides a comprehensive analysis of mathematical models used for describing epidemic processes. A main focus is placed on the models that were developed and used for modeling the COVID-19 pandemic. A huge number of studies were published since the outbreak of COVID-19 at the end of 2019. So, it is practically impossible to highlight even only the main results of many hundreds of papers. Thus, we mostly paid attention to the studies devoted to the mathematical models based on partial differential equations (typically reaction–diffusion equations) and those involving not only numerical simulations but analytical techniques as well.

The review begins by acknowledging the historical significance of the Kermack–McKendrick work, in which the SIR model was developed. It is widely accepted that this model laid foundations for mathematical modeling in epidemiology. In Section 2, we also present various extensions and generalizations of the SIR model such as the SIRD model, the SEIR model, etc. All models of this type are based on ODEs. Typically, numerical simulations have been the primary methods for solving these models due to nonlinearities in the basic ODEs. However, several studies have focused on finding exact and approximate solutions using relevant analytical techniques (Lie symmetry method, classical methods for integration of nonlinear ODEs, the homotopy analysis method, optimal auxiliary functions method, etc.). This review highlights the importance of the results obtained because exact solutions (even particular those) are very useful for qualitative descriptions of the pandemic spread and for estimating the accuracy of numerical methods. Key epidemiological parameters, like the basic reproduction number $R_0$, have been estimated using analytical results for these models.

During the COVID-19 pandemic, significant spatial heterogeneity was observed in the pandemic spread, emphasizing the importance of considering spatial aspects in mathematical modeling. It is well known that reaction–diffusion equations are applied for the mathematical description of a wide range of biomedical processes in order to take into account the spatial dynamics. Recently, these equations were used to study the spatial spread of the COVID-19 pandemic. The main part of the review, Section 3 and Section 4, is devoted to such types of models because they are the ones that allow us to better understand the spatial dynamics of epidemics.

We start from the natural generalizations of the classical models, in particular, the SIR and SEIR models with diffusion in space. Conditions for the existence of traveling waves, relations with the diffusive Lotka–Volterra system, and estimation of the basic reproduction number are discussed. More complicated models, especially those based on multicomponent systems and/or involving nonconstant diffusivities are presented as well, and their application for numerical estimations of the spatial spread of the COVID-19 pandemic is discussed. Furthermore, epidemic models involving reaction–diffusion equations with cross-diffusion and convective terms are discussed. Cross-diffusion is a known phenomenon observed in various biomedical and ecological processes. Recently, reaction–diffusion systems with cross-diffusion were applied to the mathematical modeling of the COVID-19 spread. In particular, exact solutions in the form of traveling waves were constructed, and their interpretation was provided.

Finally, the review touches upon age-structured epidemic models. Such types of models incorporate age dimension as a second time variable. Usually, relevant governing equations are integro-differential; however, their plausible approximations can be derived in the form of PDEs. In this case, systems of the first-order PDEs are obtained, and exact solutions (at least in implicit forms) can be derived. However, if one needs to take into account also spatial heterogeneity, then relevant models are much more complicated. Recently, some analytical results, especially conditions for existence of traveling fronts, were derived. However, to the best of our knowledge, there are no applications of such models for numerical simulation of the spatial spread of the COVID-19 pandemic.

It should be noted that there are some papers, see, e.g., [73,74], in which the models based on reaction–diffusion equations with time delays are suggested for describing epidemic processes. Actually, such types of models are natural extensions of the models based on ODEs with time delays (see, e.g., Chapter 10 in [1]). However, to the best of our knowledge, there are no direct applications for modeling the COVID-19 spread in those papers.

In summary, this review could serve as a valuable resource for researchers and practitioners in the field of modeling in epidemiology, offering analysis and application of a wide range of mathematical models. The main emphasis is placed on the models taking into account the spatial heterogeneity that was widely observed during the COVID-19 pandemic spread. We presented a comprehensive overview of the studies, especially those published after the outbreak of COVID-19, devoted to epidemic models and based on partial differential equations.