# Reaction-diffusion waves in biology: new trends, recent developments

V. Volpert, S. Petrovskii

<figure><img src="image_1.png" /></figure>

## To cite this version:

V. Volpert, S. Petrovskii. Reaction-diffusion waves in biology: new trends, recent developments. Physics of Life Reviews, 2025, 52, pp.1-20. 10.1016/j.plrev.2024.11.007. hal-04952506

HAL Id: hal-04952506

https://hal.science/hal-04952506v1

Submitted on 23 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

## Review

# Reaction-diffusion waves in biology: new trends, recent developments

V. Volpert a,c,*, S. Petrovskii b,c

<figure><img src="image_4.png" /></figure>

## ARTICLE INFO

Communicated by Jose Fernando Fontanari

**Keywords:**

Reaction-diffusion waves

Ecology

Evolution

Biomedical applications

## ABSTRACT

Reaction-diffusion systems are widely used in the description of propagation phenomena in biological systems, where chemical and biological processes combine to produce spatial and temporal patterns. This paper explores the recent trends and developments in the study of reaction-diffusion waves, highlighting their relevance to diverse biological contexts such as population dynamics, ecology or biomedical applications. Progress in mathematical techniques and computational methods advances our ability to model these systems, providing deeper insights into wave propagation, stability, and bifurcations. We also discuss novel models and their implications for understanding processes such as biological invasions or disease proliferation. Overall, the integration of modern theoretical frameworks with experimental data continues to push the boundaries of this interdisciplinary field, revealing new applications and mechanisms underlying biological wave dynamics.

## 1. Introduction

Reaction-diffusion waves are a ubiquitous phenomenon observed in many nonlinear systems of different origin and, correspondingly, in many areas of science such as chemistry [40,93,105], biology [64,97,100], ecology [57,89], epidemiology [55,96], neuroscience [94,106], medicine [17,64,86], biological evolution [12,23] and even social dynamics [25,76,107] and economy [11]. Reaction-diffusion wave provides a basic example of self-organisation where the emerging dynamics exhibits a certain symmetry (e.g. translational) that the original system may not have possessed [52,58]. In the travelling wave regime, instead of depending on space and time separately, the system dynamics depends on their linear combination. That reduces system's dimensionality and, generally speaking, its dynamical complexity. In mathematical terms, instead of partial differential equations (PDEs), travelling waves are described by ordinary differential equations (ODEs). Since ODEs are usually easier to study analytically than PDEs, the emergence of travelling waves allows for a more comprehensive investigation of the model; in particular, the speed of the travelling wave can often be found analytically.

As a starting point, consider the dynamics of a reaction-diffusion system in an infinite one dimensional (1D) space described by a single scalar variable (say, $u$) that comes as a solution of the reaction-diffusion equation quantified by

$$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial^2 u(x,t)}{\partial x^2} + F(u), \quad -\infty < x < \infty, \quad t \ge 0, \qquad (1.1) $$

<footnote>a Institut Camille Jordan, UMR 5208 CNRS, University Lyon 1, 69622 Villeurbanne, France</footnote>

<footnote>b School of Computing and Mathematical Sciences, University of Leicester, Leicester, LE1 7RH, UK</footnote>

<footnote>c Peoples Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya Str., 117198 Moscow, Russia</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: volpert@math.univ-lyon1.fr (V. Volpert).</footnote>
where x and t are space and time, respectively, D is the diffusion coefficient, and the nonlinear reaction term F(u) describes the local kinetics of the system. In applications, u often has the meaning of concentration or density (e.g. the concentration of a chemical substance, population density, etc.), so that usually u ≥ 0. Under certain assumptions about function F and the initial and boundary conditions, Eq. (1.1) has a travelling wave solution, so that u(x, t) = g(ξ) where ξ = x − ct is called the travelling wave coordinate, the constant c (that can be of a different sign) has the meaning of the speed of the wave. In particular, in case function F(u) has at least two zeros, say F(0) = F(ū) = 0 (which correspond to steady states of the corresponding nonspatial system) and the initial conditions are given by a piece-wise constant function, u(x, 0) = 0 for x < 0 and u(x, 0) =ū for x ≥ 0, the emerging travelling wave (to which the initial condition promptly converges) has a form of front or kink connecting the two steady states u = 0 and u =ū, its speed and the direction of propagation depending on details of function F(u). This fundamental result can be generalised onto the case of a system of reaction-diffusion equations [99], however in that case the emerging travelling wave may have a shape much more complicated than a kink [100], e.g. exhibiting multiple oscillations at the front [74].

In a general case, the reaction term F can also depend on space and time. The dynamics of the system then becomes much richer. In a strict sense, a stationary travelling wave does not exist anymore; however, in many cases its generalisation such as a non-stationary travelling wave can be considered [14,95,99,100,105]. The solution of the corresponding reaction-diffusion equation takes a somewhat different form, e.g. as u(x, t) = g(x − ct, t), and it can then describe a variety of new phenomena. For instance, the dynamics may be regarded as perturbations (which can be solitary, periodic or irregular) of a stationary travelling wave, or else in the course of time the travelling front can be slowed down, blocked or reversed.

Another way to generalise the traditional reaction-diffusion framework comes from a more advanced analysis of either the transport (movement) processes or the local kinetics of the system. In the former case, the simple Brownian diffusion (as described by the Laplacian term in Eq. (1.1)) may need to be replaced by a more complicated movement mode, e.g. to account for nonlocal diffusion or cross-diffusion or nonlinear advection. In the latter case, the pointwise ‘local’ reaction term may become nonlocal. In either of the above cases, the corresponding reaction-diffusion-type system may exhibit a variety of travelling waves with different properties. One of the motivations behind this paper is to provide an overview of some recent research in those directions.

Reaction-diffusion waves, in spite of being a focus of intense research for several decades, e.g. see [100] and further references there, continue attracting considerable attention. Along with widening the list of applications where reaction-diffusion waves play an important role and broadening the context in which reaction-diffusion waves can occur, over the last fifteen years increasing attention has been paid to further generalisation of the traditional mathematical framework and to, correspondingly, revealing the travelling wave-type dynamics with new properties.

In this review paper, we endeavour to provide an overview of some recent work on reaction-diffusion waves where a traditional analysis and applications of such waves was combined with some new approaches and/or new applications. The paper is structured according to the applications. In Section 2, we focus on reaction-diffusion waves arising in ecological, environmental and social systems. In Section 3, we reveal how the concept of reaction-diffusion waves may provide a possible mechanism for biological evolution and the emergence of new species. Section 4 mostly focuses on biomedical applications.

## 2. Population dynamics, ecology and environment

### 2.1. Travelling wave in a system of competing populations: movement response can change the direction of propagation

Interspecific competition is a common ecological interaction that to a large extent shapes the structure and dynamics of the population community [91]. One important property of a community of competing species is the competitive exclusion principle [37, 43]. Its exact formulation remains a subject of discussion and can be context-specific; in a not very rigorous form, it can be stated as “two species competing for the same resources cannot stably coexist if other ecological factors are constant” [45]. In practical terms, it often means that even a small competitive difference may results in the extermination of the weaker competitor, especially when the species are considered on a longer (evolutionary) timescale.

Because of their apparent high ecological relevance, spatially explicit mathematical models of competing populations have been a focus of attention and intense research, e.g. see [59,65,66,75,80,87]. In the baseline case of two competing populations - say, u(x, t) and v(x, t) where x and t are space and time, respectively - their spatiotemporal dynamics can be described by the following reaction-diffusion system:

$$ \frac{\partial u(x, y, t)}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + r_1 u \left(1 - \frac{u}{K}\right) - a_1 u v, \quad (2.1) $$

$$ \frac{\partial v(x, y, t)}{\partial t} = D \frac{\partial^2 v}{\partial x^2} + r_2 v \left(1 - \frac{v}{K}\right) - a_2 u v, \quad (2.2) $$

(cf. [65]) where it is assumed for the sake of simplicity that the diffusion coefficient D and the carrying capacity K are the same for both populations. System (2.1)-(2.2) is known to have a travelling wave solution (e.g. see [100]), which was successfully applied to model and explain a real-world ecological phenomenon [66].

System (2.1)-(2.2), although having been widely used in many theoretical and modelling studies, apparently has a number of assumptions and limitations. One of them, which we are going to highlight here, is that the movement of the two populations is entirely independent of each other. The interaction between u and v only occurs through the reaction term. In reality, however, this is often not quite true. As just one example, it has recently been observed that top-predators change their movement when close to competing predators; see [80] for details and further references. A question therefore arises as to whether a similar kind of
<figure><img src="image_5.png" /><figcaption>Fig. 1. (left) Wave profiles arising from step-wise initial conditions; (right) The structure of $(a_1, \gamma_2)$ parameter plane where different domains correspond to different direction of the travelling front propagation. Adapted from [80].</figcaption></figure>

behavioural response may change the outcome of the competition, in particular in the parameter range where it would otherwise lead to the extinction of the ‘weaker’ population (i.e., the weaker competitor) due to the competitive exclusion. The key idea is that movement responses may be used by an otherwise weaker population to push back the travelling wave, causing a reversal in the eventual fate of the populations.

Mathematically, the movement responses described above are described by an advection term, turning system (2.1)-(2.2) to a system of reaction–advection–diffusion equations. Taking the above arguments into account and choosing appropriate dimensionless variables, see [80] for details, Eqs. (2.1)-(2.2) (where we additionally assume that $r_1 = r_2 = r$) are transformed into the following system:

$$ \frac{\partial u(x, y, t)}{\partial t} = \frac{\partial^2 u}{\partial x^2} - \gamma_1 \frac{\partial}{\partial x} \left[ u \frac{\partial v}{\partial x} \right] + ru(1-u) - a_1 uv, \quad (2.3) $$

$$ \frac{\partial v(x, y, t)}{\partial t} = \frac{\partial^2 u}{\partial x^2} - \gamma_2 \frac{\partial}{\partial x} \left[ v \frac{\partial u}{\partial x} \right] + rv(1-v) - a_2 uv. \quad (2.4) $$

Consider a biologically inspired choice of the initial conditions where one half of the domain is occupied by species u and the other half by species v. These initial conditions promptly converge to an interface that separates the domains where either u or v dominates; see Fig. 1 (left). Exact form of the emerging interface depends on the properties of the nonspatial system and hence on coefficients $a_1$ and $a_2$; the profiles shown in Fig. 1 (left) correspond to the case where the nonspatial system does not possess a positive ‘coexistence’ steady state but has two boundary states (0, 0) and (0, 1) where only one of the species survives, cf. [59,64]. The interface between the two sub-domains appears to be a travelling wave; in the system without advection terms (i.e., Eqs. (2.3)-(2.4) in the special case $\gamma_1 = \gamma_2 = 0$), the direction of its propagation depends on the relative competition strength. Numerical simulations (supported also by some analytical results, e.g. see [38,47,48]) suggest that in case $a_1/a_2 < 1$ (i.e., species u is a stronger competitor), this travelling wave moves to the right, so that species u wins. However, in the presence of movement responses (modelled by the advection terms), the dynamics can change dramatically. In particular, an increase in $\gamma_2$ (which can biologically be interpreted as a more aggressive behaviour) can reverse the direction of front propagation from right to left, so that species v becomes the winner (see Fig. 1, right).

Interestingly, while the change in the propagation direction observed for higher values of $\gamma_2$ might be intuitively expected, the dependence of the critical advection rate $\gamma_2^*$ (i.e., the value of $\gamma_2$ when the propagation direction changes) on the population dynamics parameter such as the growth rate r and the relative competition strength appears to be less intuitive. An example is shown in Fig. 2.

## 2.2. Travelling waves blocking by a narrow passage

Biological invasion is a biological phenomenon where reaction-diffusion waves arise in the most natural way [57,89,90] and this is indeed the context where the reaction-diffusion waves were originally discovered [36,51]. An alien (exotic) population introduced locally may start proliferating and spreading into space; such spread is well described by a travelling wave. Most theoretical studies into this phenomenon consider such spread in a uniform geographical space (but see [49,72,88]) hence focusing on other aspects such as different patterns of individual movement, etc. Correspondingly, the question how species spread can be affected by an environmental heterogeneity, in particular in habitats of complex shape, remains somewhat under-investigated. In this section, we address this problem by considering biological invasion in an environment of complex “H-shaped” geometry where two large areas are connected by a narrow corridor; see Fig. 3.
<figure><img src="image_6.png" /><figcaption>Fig. 2. Dependence of the critical advection rate $\gamma_2^*$ on the relative population growth rate $r/a_2$ for different values of the relative competition strength $a_1/a_2$. Each of the curves therefore divides parameter plane $(a_1/a_2, \gamma_2)$ into two parts. For parameters from below the curve (i.e., for $\gamma_2 < \gamma_2^*$) the front propagates to the right – species u wins; for parameters from above the curve (for $\gamma_2 > \gamma_2^*$) the front propagates to the left – species v wins. Adapted from [80].</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 3. A sketch of the H-shaped domain consisting of Habitat 1 at the left and Habitat 2 at the right connected by a narrow passage (corridor). The alien species u is introduced in Habitat 1, the initially inhabited area is shown by a small yellow rectangle. The initial distribution of species v (which can be regarded as a biological control agent) has a similar shape.</figcaption></figure>

To describe the population dynamics, we use the following spatially explicit prey-predator system:

$$ \frac{\partial u(x, y, t)}{\partial t} = \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) + \gamma u(u - \beta)(1 - u) - \frac{uv}{1 + \alpha u}, \quad (2.5) $$

$$ \frac{\partial v(x, y, t)}{\partial t} = \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) + \frac{uv}{1 + \alpha u} - \delta v, \quad (2.6) $$

where $\alpha$, $\beta$, $\gamma$ and $\delta$ are dimensionless parameters; see [1] for details. System (2.5)-(2.6) is rather generic and can be applied to a variety of ecological situations; here, we interpret species u as an alien (invading) species and species v as a biological control agent which is deliberately introduced in an attempt to slow down or stop the spread of u.

In the first stage of invasion, the invading population u spreads into space (with species v spreading in the wake) in the same manner as in an infinite space [57]; see Fig. 4 top, left. Once it reaches the entrance to the corridor, the spread slows down. Depending on the size of the corridor, it may take a considerable extra time to penetrate through it (cf. Fig. 4). Eventually, the population spills over to Habitat 2 and invade over the rest of the domain.<fnref n="1" /> However, the latter does not necessarily happen. Depending on the width

<footnote n="1">¹ Note that the pattern of spread on the other side of the passage (i.e., in Habitat 2) is different from Habitat 1: the narrow passage restores the symmetry in the population distribution, so that, instead of the irregular, patchy spread in Habitat 1, the population spread in Habitat 2 takes the form of a travelling wave with apparent circular symmetry.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 4. Snapshots of the spatial distribution of prey at different moments obtained by solving numerically Eqs. (2.5)–(2.6) for parameters $\gamma = 3.9$, $\alpha = 0.1$, $\beta = 0.2$ and $\delta = 0.5$ and compact initial condition (cf. Fig. 3). Parameters of the corridor are $(\lambda, w) = (10, 2.8)$. Yellow for high population density, dark blue for the zero density, green for intermediate density. Adapted from [1].</figcaption></figure>

($w$) and length ($\lambda$) of the corridor, the system can exhibit a different dynamics resulting in an invasion failure. An increase of the corridor width will not change the dynamics much resulting only in a decrease of the time required for the population to penetrate through. However, a decrease of the corridor width will eventually result in invasion blocking: as soon as the width is smaller than a certain critical value, the invading population is blocked in the passage and never spills over to Habitat 2.

It appears possible to obtain an analytical estimation of the critical width. Using the comparison principle for nonlinear partial differential equations [81,98], it is readily seen that the following single-species model:

$$ \frac{\partial u(x, y, t)}{\partial t} = \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) + \gamma u(u - \beta)(1 - u), \qquad (2.7) $$

i.e., Eq. (2.5) without the last term in the right-hand side, is the upper bound for the prey-predator system (2.5)–(2.6). Furthermore, once the front arrives to the end of the passage and starts spilling into Habitat 2, it is of approximately circular shape (as is clearly seen in numerical simulations, e.g. see the bottom-left panel in Fig. 4). For a circular propagating front of the radius $R$ in a 2D space, its speed $c(R)$ is known to be related to the speed $c$ of the corresponding 1D front as

$$ c(R) = c - \frac{D}{R}, \qquad (2.8) $$

(e.g. see [100]) where $D$ is the diffusion coefficient. For the front confined at the exit of the passage, $R = w/2$. Since the front is blocked, $c(R) = 0$. Equation (2.8) then turns into the following (where we recall that in our choice of dimensionless variables $D = 1$):

$$ c - \frac{2}{w} = 0. \qquad (2.9) $$

For a scalar reaction-diffusion equation with the population growth parameterized by a cubic polynomial, the analytical expression for the speed of the front propagation is well-known:

$$ c = \sqrt{\frac{\gamma}{2}(1 - 2\beta)}, \qquad (2.10) $$
e.g. see [57]. From (2.9) and (2.10), we therefore obtain:

$$w_{cr} = \frac{2\sqrt{2}}{(1 - 2\beta)\sqrt{\gamma}} \qquad (2.11)$$

Estimate (2.11) agrees well with the simulation results [1].

A different analysis can be used in case of a long narrow passage. It also appears possible to obtain a simple analytical estimate of the critical corridor width, although the argument and calculations are different. We recall that the single-species model is an upper bound for the prey-predator model, so a condition of population propagation blocking in the model (2.7) will be a sufficient condition of the propagation blocking in model (2.5)-(2.6). Consider the case of a long corridor. Consider the situation where the invasive species has penetrated through the corridor and successfully spilled over to Habitat 2. In this case, in the large time limit the population density profile across the corridor does not depend on its position along the corridor as far as the position is sufficiently far away from the corridor ends. The problem then becomes one-dimensional as the population density depend only on the coordinate across the corridor but not along the corridor. If considered in the corresponding 1D settings, this turns into the classical critical size problem [50] (also see Section 3.1 in [57]), that is, how large the domain must be to ensure the survival of the population inside. For the growth rate described by a cubic polynomial, its analytical solution is unknown. However, one can consider an upper bound of the growth rate as

$$\gamma u(u - \beta)(1 - u) \leq gu \quad \text{for any } u \geq 0. \qquad (2.12)$$

Now, if we choose g in such a way that the linear function gu is the tangent line for the graph of the nonlinear growth rate in the left-hand side of (2.12) (in this case the linear function becomes the exact upper bound), then g appears to be related to γ as $γ = 4g/(1 - \beta))^2$ by the following equation [56]. Using the upper bound instead of the actual nonlinear growth rate, the expression for the critical size in the corresponding linear problem is well known:

$$w_{cr} = \frac{2\pi}{(1 - \beta)\sqrt{\gamma}} \qquad (2.13)$$

(in the dimensionless variables). Estimate (2.13) gives the value $w_{cr}$ which is in a good agreement with numerical results [1].

Note that the mechanism of spread blocking is somewhat different in the two cases. For the zero-density conditions at the boundary of the passage, the population is blocked inside the passage (in a manner similar to how a subcritical size of the inhabited domain results in the population extinction). For the zero-flux conditions, however, the population cannot be blocked inside (except for the trivial case where $\beta > 0.5$); instead, its propagation is blocked at the exit of the passage.

## 2.3. Travelling waves in the plankton-oxygen system: growth of an Oxygen Minimum Zone and onset of global anoxia

Distribution of dissolved oxygen (DO) in the ocean has been attracting considerable attention over the last two decades [67,69,78]. There are several reasons for that. Firstly, DO is an important indicator of a marine ecosystem health, as a decay in the DO availability may result in anoxic condition, which in turn would lead to a mass mortality of marine fauna. Secondly, a considerable part of the DO produced by phytoplankton will go to the air through the ocean-air surface: it is estimated that about two-thirds of the total oxygen budget in the atmosphere originated in the ocean. And thirdly, it was observed in several recent studies that the average level of DO in the ocean showed a clear tendency to decrease, this decrease being attributed to the effect of global warming [26,102].

Importantly, the spatial distribution of DO in the ocean is highly heterogeneous, both in vertical and horizontal directions. This heterogeneity often has the form of a large stable ‘patch’ (or volume) of water with anoxic or even hypoxic conditions inside. Such patch is called an Oxygen Minimum Zone (OMZ). OMZ can be as large as hundreds of meters in depths and tens or even hundreds kilometres in the horizontal direction [69,103]. Some recent field studies indicated that at least some of the OMZs tend to further grow in size (see [2,3] for details and further references). This growth was attributed to the effect of global warming but no specific mechanism was suggested.

The processes resulting in the OMZ formation and maintenance have been a focus of intense study over the last two decades. One line of research endeavours to link the phenomenon of OMZ to the specifics of ocean mixing in a stratified ocean combined with differential oxygen solubility at different water temperatures. Although apparently plausible, this explanation entirely disregards the fact that considerable amount of oxygen is in fact produced in the ocean in phytoplankton photosynthesis. Correspondingly, an alternative (or complimentary) line of research attempts to link the OMZ formation to the nonlinear dynamics of the plankton-oxygen system [70,79].

A mathematical model of coupled plankton-oxygen dynamics was developed and studied in [77,82] and further investigated in [39,68,83–85]. In the baseline nonspatial case (e.g. in a well-mixed photic ocean layer), it consists of the following two equations:

$$\frac{dc}{dt} = \frac{Au}{c + 1} - \frac{\delta u c}{c + c_2} - c, \qquad (2.14)$$

$$\frac{du}{dt} = \left( \frac{Bc}{c + c_1} - u \right) u - \sigma u, \qquad (2.15)$$

(in dimensionless variables, see [77,82] for details) where c is the concentration of oxygen, u is the density of phytoplankton, t is time, and A is the rate of net oxygen production which is known to depend on water temperature.
<figure><img src="image_9.png" /><figcaption>Fig. 5. (Left) The nullclines of the oxygen-phytoplankton system (2.14)-(2.15). Black curves show the first (oxygen) nullcline obtained for $c_2 = 1$, $\delta = 1$ and $A = 3, 5, 6$ and 7 (from left to right, respectively); red curve shows the second (phytoplankton) nullcline obtained for $B = 1.8$, $\sigma = 0.1$ and $c_1 = 0.7$. (Right) The steady state values of $c$ (obtained as the intersection points of the two nullclines) as a function of the controlling parameter $A$. Adapted from [82].</figcaption></figure>

It is readily seen that system (2.14)-(2.15) always possesses the trivial ‘extinction’ steady state (0,0) which is always stable. Furthermore, subject to some restrictions on the parameter values, in particular for the rate of net oxygen production being not too small, $A > A_{cr}$, system (2.14)-(2.15) possesses two positive steady states (a stable node and a saddle); see Fig. 5 and [82] for more details. If $A$ falls below a critical value $A_{cr}$, the system passes a tipping point as the two positive steady states disappear in a saddle-node bifurcation. Thus, for $A < A_{cr}$ the extinction steady state is the global attractor, so that, regardless of the initial conditions the system dynamics inevitably results in the oxygen depletion and plankton extinction. Some field observations indicate that $A$ may decrease with an increase in water temperature. Correspondingly, an increase in the average water temperature in the ocean upper layer over a certain threshold may disrupt phytoplankton photosynthesis and hence eventually lead to a global anoxia. This theoretical prediction originally made in [82] using the toy model (2.14)-(2.15) was later generalised and extended onto a class of more realistic models in [77].

An important question arises here as to how close the ocean may be to the tipping point where the oxygen production by phytoplankton stops entirely. This issue can be addressed by linking the properties of model (2.14)-(2.15) to the observed OMZ growth, in particular by considering the OMZ boundary as a travelling wave [3].

A spatially explicit extension of model (2.14)-(2.15) can be obtained by adding diffusion terms to account for the effect of turbulent mixing:

$$ \frac{\partial c(x,t)}{\partial t} = D \frac{\partial^2 c}{\partial x^2} + \frac{Au}{c+1} - \frac{\delta u c}{c+c_2} - c, \qquad (2.16) $$

$$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + \left( \frac{Bc}{c+c_1} - u \right) u - \sigma u, \qquad (2.17) $$

where the oxygen concentration and the phytoplankton density now also depend on the position in space $x$. Note that, generally speaking, the horizontal dynamics in an ocean layer is two-dimensional. Here, for the sake of simplicity, we consider a 1D system; this can be interpreted as a transect along the ocean surface.

In the reaction-diffusion system (2.16)-(2.17), initial conditions mimicking an OMS (e.g., with low values of DO and plankton density in one half of the domain and high values in the rest of the domain) converge to a travelling front, see Fig. 6 (left). Interestingly, the direction of the front propagation appears to depend on $A$ (Fig. 6, right): there is a certain critical value $A_*$ such that the front propagates to the right for $A > A_*$ but it propagates to the left for $A < A_*$. For our choice of the initial condition, propagation to the right means that the OMZ shrinks, while propagation to the left means that the OMZ grows. Note that $A_* > A_{cr}$, so that the OMZ growth through this travelling front scenario occurs in the parameter range where the corresponding nonspatial system still possesses a stable positive steady state.

Thus, in terms of the above modelling approach, the OMZ growth observed in real ocean may indicate that the system has already passed through the first critical value $A_*$ and hence is likely to be dangerously close to the second critical value $A = A_{cr}$. Recall that $A = A_{cr}$ is a tipping point where the oxygen production stops, surpassing it would inevitably lead to global anoxia.

Model (2.16)-(2.17) is apparently very schematic, leaving out of scope many factors and processes. In particular, in a real ocean the phytoplankton abundance is usually controlled by zooplankton as its main consumer (predator), and this was shown to have a significant effect on the spatiotemporal dynamics of dissolved oxygen, e.g. see [77,82–84]. Zooplankton usually has lower growth rate and longer lifespan, so that altogether its dynamics occurs on a slower timescale. In its turn, zooplankton abundance is controlled by higher trophic layers, e.g. fish. Correspondingly, a somewhat more detailed model was recently considered in [29] consisting of the following equations:

$$ \frac{\partial c}{\partial t} = \frac{\partial^2 c}{\partial x^2} + \frac{Au}{c+1} - \frac{\delta u c}{c+c_2} - \frac{v c v}{c+c_3} - c, \qquad (2.18) $$
<figure><img src="image_10.png" /><figcaption>Fig. 6. (left) Numerical solution of system (2.16)-(2.17) obtained for A = 1.1 showing the travelling front of oxygen (blue) and phytoplankton (red) obtained at different times (from left to right): t = 500, 1000, 1500 and 2000. (right) Speed of the front as a function of parameter A. The direction of front propagation changes at A = A* (shown by the big black dot). In parameter Domain II, the nonspatial system possesses a stable positive steady state; in parameter Domain I, the positive steady state disappears through a saddle-node bifurcation (cf. Fig. 5, right). Adapted from [3].</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7. Spatiotemporal plot of the growth of the initial OMZ simulated in the spatial domain 0 ≤ x ≤ 500 over time 0 < t ≤ 2500, left and right panels obtained for two different parameter sets. Apparently, in both cases the OMZ boundaries propagate as transient travelling waves, eventually leading to oxygen depletion all over the domain (as is shown by deep blue colour in the top part of the panels). Adapted from [29].</figcaption></figure>

$$ \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \left( \frac{Bc}{c+c_1} - u \right) u - \frac{uv}{u+h} - \sigma u, \qquad (2.19) $$

$$ \frac{\partial v}{\partial t} = D \frac{\partial^2 v}{\partial x^2} + \epsilon \left( \frac{\eta c^2}{c^2 + c_4^2} \frac{uv}{u+h} - \mu_1 v - \mu_2 v^2 \right), \qquad (2.20) $$

(in dimensionless units). Here new parameter $\epsilon < 1$ accounts for the difference in the intrinsic timescales of phyto- and zooplankton (hence potentially making Eqs. (2.18)-(2.20) a slow-fast system, see [29] for details), the quadratic term $\mu_2 v^2$ implicitly accounts for the effect of higher trophic levels (cf. [32]). D is the ratio of the diffusion coefficients and it can be larger as well as smaller than unity to account for zooplankton own ability to move.

Questions arise here as to (i) whether predictions from the more complex (and, arguably, more realistic) model (2.18)-(2.20) of the OMZ dynamics are consistent with those obtained in a simple model (2.16)-(2.17), i.e., that OMZ growth may potentially lead to global anoxia, and (ii) whether reaction-diffusion waves remain relevant as a possible scenario of the OMZ dynamics.

A thorough numerical investigation of system (2.18)-(2.20) reveals that the answer to both above questions is yes. Examples of the spatiotemporal evolution of the initial OMZ (located in the middle of the spatial domain) are shown in Fig. 7. It is readily seen that, in both cases, the OMZ growth occurs because its left and right bounds move away as travelling waves. However, differently from the predictions of model (2.16)-(2.17) (see Fig. 6), in either of the two cases presented in Fig. 7, these travelling waves are clearly not stationary. The corresponding system dynamics appears to be essentially transient as the waves suddenly disappear at a certain moment resulting in oxygen depletion over the entire domain.
## 2.4. Transient population waves in the dynamics of street protests

In this section, we will discuss how reaction-diffusion waves may emerge in a somewhat different, less traditional context, namely, in the population dynamics of street protests. Modelling of street protests, manifestations, riots etc., as well as civil unrest more generally, has been attracting increasing attention recently [10,28,31,101] and reaction-diffusion equations were shown to be an efficient mathematical technique to describe this phenomenon [9,76].<fnref n="2" /> Below we will show that, while spatial dynamics of street protests exhibits many features common with reaction-diffusion systems of other origin, the specifics of social interactions can also bring new or more exotic properties such as transient auto-collapsing population waves.

Following [24], we describe the protest by a single variable (say, N). In the nonspatial settings, this is the total number of people actively participating in the street actions. Obviously, N depends on time. In a spatially explicit model, the corresponding variable (say, u) has the meaning of the ‘population density’ of protest participants and depends on time (t) and space (x).

Let us begin with a nonspatial framework. The number N of participants changes with time as a result of the interplay between recruitment (people joining the protest) and withdrawal (people quitting the protest). Recruitment is believed to be a collective, density dependent phenomenon (e.g. see [41]) so that the recruitment rate is a nonlinear function of the number of people currently involved in the event, say F(N). The decision to withdraw from the protest (quantified by the withdrawal rate m) would normally be made individually, however it can also depend on N through some implicit feedbacks (see below).

Correspondingly, we consider the following generic model of the population dynamics of street protests [76]:

$$ (a) \quad \frac{dN(t)}{dt} = F(N) - mN, \quad (b) \quad \frac{dm(t)}{dt} = G(m, N), \qquad (2.21) $$

where the recruitment rate F can be parameterised as follows:

$$ F(N) = \epsilon N + \frac{aN^2}{h^2 + N^2}, \qquad (2.22) $$

where the first and second terms correspond to, respectively, the contribution of individual decision making and collective effects [76]; $\epsilon, a$ and $h$ are nonnegative parameters.

Further analysis essentially depends on the assumptions about the dependence of the withdrawal rate on other factors, i.e., about function G. In the rest of this section, we will focus on two special cases. We first consider the case where collective effects can be neglected, so that G does not depend on N. With regard to its dependence on m, we take into account that, in the modern society, even in the absence of major protests or unrest, a number of small-scale protests are almost always happening. We therefore assume that there is a certain background rate (say, $m_0$) at which the participants of those protests quit the actions. We then further assume that the society is in a stable equilibrium, so that any perturbation from the background withdrawal rate should decay with time. The simplest equation that describes the corresponding evolution of m is given as follows:

$$ \frac{dm(t)}{dt} = b(m_0 - m), \qquad (2.23) $$

where b is a coefficient. Equation (2.23) is of course readily solved resulting in

$$ m(t) = m_1 + (m_0 - m_1)(1 - e^{-bt}), \qquad (2.24) $$

where $m_1 = m(0)$ is the initial value of the withdrawal rate. Note that in this case system (2.21) is reduced to just one non-autonomous equation (2.21a), with m being a known function of time.

In the second special case, we consider the opposite situation where function G does not depend on m but depends only on N. This case may arise as an effect of a ‘collective guilt’. Recall that street actions often bring considerable economic damage which can be direct (e.g. broken shop windows, burnt cars, etc.) and/or indirect (e.g. disrupting the public transport). We assume that the total amount of damage is proportionate to the number of protests participants. We then further assume that participants eventually develop the feeling of guilt for causing damage to the community, and that this feeling is likely to cut short the duration of the participation in the street actions. Correspondingly, that leads to the expression for the withdrawal rate as $G = \beta N$, $\beta$ being a coefficient.

In order to turn the nonspatial model to a spatial one, we have to take into account human movement. It was previously shown [27] that most of the everyday human travel occurs within a short range from their homes. We therefore neglect the long-distance travel and assume that the short-distance human movement (e.g. relocation between different protest sites) can be described by Brownian diffusion. In the first special case considered above, the corresponding spatially explicit model is then given by the following equation:

$$ \frac{\partial u(x, t)}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + \left( \epsilon + \frac{au}{h^2 + u^2} \right) u - m(t)u, \qquad (2.25) $$

(restricting our analysis to 1D space) where D is the diffusion coefficient and the dependence of the withdrawal rate m on time is given by (2.24).

Equation (2.25) is supplemented with compact initial condition, i.e., we assume that protests start at a certain location or inside a small area. Some typical solutions of Eq. (2.25) are shown in Fig. 8. We readily notice that the compact initial condition fast develops into a travelling front (more precisely, in the 1D space, two fronts propagating in opposite directions away from the place of the

<footnote n="2">² For a similar approach but dealing with the dynamics in a discrete space or network, see also [15,25].</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 8. Solution of the model (2.25) shown at different moments of time (as is explained in the figure legend) obtained for $m_1 = 3$ and different value of the final withdrawal rate: (left) for $m_0 = 5.4$, (right) for $m_0 = 5.8$. Other parameters are $\epsilon = 3$, $a = 150$, $h = 35$, $b = 0.01$ and $D = 1$. Adapted from [76].</figcaption></figure>

original protesters location). However, the front is non-stationary as the population density behind the front decreases with time. Moreover, a more careful look at the solution reveals that the speed of the front is not constant but decreases with time as well. As a result, eventually the direction of propagation changes to the opposite: instead of growing, the space occupied by the protesters starts shrinking, which eventually leads in the protests disappearing. We refer to this dynamical behaviour as travelling waves auto-collapse.

In the second special case, the spatially explicit model is given by the following equations:

$$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + (\epsilon - m)u + \frac{au^2}{h^2 + u^2}, \quad (2.26) $$

$$ \frac{\partial m(x,t)}{\partial t} = D_m \frac{\partial^2 m}{\partial x^2} + \beta u. \quad (2.27) $$

For the sake of generality, we assume that the spatial dynamics of the withdrawal rate may occur with a different diffusion coefficient $D_m$.

In case system (2.26)-(2.27) is complemented with compact initial conditions (hence assuming that the protests started locally), reaction-diffusion waves emerge. Some typical results are shown in Fig. 9. Thus, different assumptions about the evolution of the withdrawal rate (i.e., function $G(m, N)$) lead to reaction-diffusion waves with rather different properties. Interestingly, although system (2.26)-(2.27) is autonomous, the emerging travelling waves are non-stationary, as it is readily seen from the shape of the wave profile for variable $m$ (see Fig. 9, bottom).

For the travelling wave shown in Fig. 9, one can obtain an analytical estimate of the wave speed using an approach based on the comparison principle for nonlinear PDEs [81,98]. We first notice that Eq. (2.26) can be written as the following differential inequality:

$$ \frac{\partial u(x,t)}{\partial t} \le \epsilon u + \frac{au^2}{h^2 + u^2} + D \frac{\partial^2 u}{\partial x^2}. \quad (2.28) $$

We now build an auxiliary function $\tilde{F}(u)$ as follows:

$$ \tilde{F}(u) = \alpha u \quad \text{for} \quad 0 \le u \le u_A, \quad (2.29) $$

$$ \tilde{F}(u) = \epsilon u + \frac{au^2}{h^2 + u^2} \quad \text{for} \quad u > u_A, \quad (2.30) $$

where $\alpha$ is the slope of the tangent line to the graph of the reaction term and $u_A$ is the abscissa of the conjunction point. After some simple algebra, we obtain that

$$ \alpha = \frac{a}{2h} + \epsilon. \quad (2.31) $$

Obviously, $\tilde{F}(u)$ is an upper bound for the reaction term in the right-hand side of (2.28), so that the following inequality holds:

$$ \frac{\partial u(x,t)}{\partial t} \le \tilde{F}(u) + D \frac{\partial^2 u}{\partial x^2}. \quad (2.32) $$

By the virtue of the comparison principle, the solution of the following equation

$$ \frac{\partial \tilde{u}(x,t)}{\partial t} = \tilde{F}(\tilde{u}) + D \frac{\partial^2 \tilde{u}}{\partial x^2} \quad (2.33) $$

gives an upper bound for the solution $u(x,t)$ of Eq. (2.26) provided the initial and boundary conditions are the same.

Function $\tilde{F}(u)$ is convex, therefore the linearization at the leading edge can be applied resulting, taking (2.31) into account, in the following upper bound for the speed of the wave:
<figure><img src="image_13.png" /><figcaption>Fig. 9. (top) Population density vs space obtained as a solution of the model (2.26)–(2.27) shown at different moments of time. (bottom) The corresponding distribution of the withdrawal rate over space. Parameters are $\epsilon = 1$, $a = 115$, $h = 20$, $\beta = 0.00025$, $D = 1$ and $D_m = 0.1$. Adapted from [76].</figcaption></figure>

$$c = 2\sqrt{\alpha D} = 2\sqrt{\left(\frac{a}{2h} + \epsilon\right)D} \qquad (2.34)$$

Expression (2.34) appears to be in a good agreement with the values of speed obtained in simulations [76]. Although it somewhat overestimates the actual value (which is of course hardly surprising as (2.34) gives an upper bound), it correctly describes the dependence of the speed on $\epsilon$ showing an approximately linear increase with an increase in $\epsilon$.

In conclusion to this section, we mention that the propagating spike or ‘pulse’ of protesters (qualitatively similar to the one shown in Fig. 9, top) was also observed in a different model of protests using somewhat different mathematical techniques [25]. However, the protests dynamics described by the auto-collapsing travelling front shown in Fig. 8 (see also [76]), to the best of our knowledge, is novel and has not been previously reported in the literature.

## 3. Modelling biological species

### 3.1. Emergence and evolution of species

One of the questions discussed in [100] is related to sympatric speciation characterized by the emergence of biological species without genetic or geographic barriers. This is the keystone of Darwin’s theory largely debated in the literature. Let us recall that it is modelled in [100] with a nonlocal reaction-diffusion equation where the diffusion term characterizes small random variations, and nonlocal logistic term describes reproduction with the intraspecific competition for resources. This competition determines natural selection. The space variable in this model corresponds to animals’ phenotype, according to Darwin’s definition of biological species as a group of individuals with a close phenotype. The simplest model represents the reaction-diffusion equation

$$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + au^k(1 - J(u)) - \sigma u, \qquad (3.1)$$

with the integral term

$$J(u) = \int_{-\infty}^{\infty} \phi(x - y)u(y, t)dy$$

which describes the intraspecific competition by means of nonlocal consumption of resources. The kernel $\phi(x)$ can be considered as a step-wise constant function

$$\phi(x) = \frac{1}{2N} \begin{cases} 1, & |x| < N \\ 0, & |x| \ge N \end{cases}.$$
<figure><img src="image_14.png" /><figcaption>Fig. 10. Darwin's diagram (in the centre) and numerical simulations around it showing the main evolutionary patterns (see the explanation in the text). Reprinted from [22].</figcaption></figure>

In the formal limit $N \to 0$, we obtain conventional (local) reaction-diffusion equation with logistic reproduction term. In this case, the homogeneous in space stationary solution $u = 1$ is asymptotically stable independently of the diffusion coefficient. However, this solution can become unstable for equation (3.1). Indeed, linearising this equation about the stationary solution, and applying the Fourier transform, we obtain the following expression for the spectrum:

$$ \lambda = -D\xi^2 - a \frac{\sin(\xi N)}{\xi N}, \quad \xi \in \mathbb{R} $$

$(k = 1, \sigma = 0)$. In the limit $N \to 0$, we have $\sin(\xi N)/(\xi N) \to 1$, and $\lambda < 0$ for all values of $\xi$. However, for finite and sufficiently large value of $N$ and sufficiently small $D$, $\lambda$ becomes positive in some interval of $\xi$. Therefore, the homogeneous in space stationary solution $u = 1$ loses its stability resulting in bifurcation of periodic in space solutions. In this case, equation (3.1) describes propagation of periodic waves. As discussed in [100], they correspond to the emergence of biological species.

Darwin's diagram shows different evolutionary patterns, such as speciation, competition, extinction, stable persistence (Fig. 10). All of them can be reproduced in the model based on nonlocal reaction-diffusion systems of equations [22,23]. Each line in the diagram corresponds to one species, and it is denoted by a letter. The horizontal axis corresponds to some phenotypical characterization of the species (e.g., animal's size, the beck of the bird, etc.), the vertical axis corresponds to time. So, the diagram illustrates how different species evolve in the phenotype space.

In the model, each species is characterized by a narrow population density distribution (spike) with respect to its phenotype. The spatial coordinate of this spike corresponds to the species location in the phenotype space, while its time dependence describes the species evolution (see [22,23] for more detail).

Consider some examples from Fig. 10. The most frequent pattern in the diagram corresponds to the emergence of several species from one single species. It is shown as several traits starting from a single point. After some time, all of them disappear except for two extremal left and right ones. The biological explanation of this effect is that they have weaker competition and better chances to survive. In the model, in order to reproduce this effect, we consider space-dependent coefficients describing how the reproduction or survival rates depend on the phenotype.

Another example shows extinction of species. As such, species B, C, D extinct because they enter in competition with species A which evolves in the phenotype space and occupies the niches of other species. It is interesting to note that in the simulations (but not in the diagram), these species try to escape the evolutionary pressure changing their phenotype but finally fail and extinct. Species G, H, I, K, L coexist and slowly move from each other, while species E and F do not change their position and do not bifurcate. We will discuss this last case in the next section.

If Darwin considered biological species as groups of individual with a similar phenotype, Meyr defined them as groups of individual that can reproduce between themselves (implying that different species do not give joint offsprings). Modelling allows us to make a link between these two definitions. In this case, instead of equation (3.1) we consider the equation

$$ \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + a(S(u))^2(1 - J(u)) - \sigma u, \quad (3.2) $$

where the integral
$$S(u) = \int_{-\infty}^{\infty} \psi(x-y)u(y,t)dy$$

determines how the reproduction rate depends on the distance in the phenotype space. If $\psi$ is a delta-function, then we recover equation (3.1) describing the case in which only animals with the same phenotype can give offspring. However, if the kernel support is sufficiently large, such that animals from different species can interbreed, then analysis of this equation shows that the species do not emerge. Thus, if interbreeding is allowed in the model, that is, there are no species in Meyr's sense, then the species in Darwin's sense do not emerge neither [12].

### 3.2. Existence and bifurcation of pulses

From the mathematical point of view, as discussed above, biological species can be considered as population density distributions with respect to their phenotype or genotype. Existence and stability of solutions describing such distributions determine the relevance of the corresponding model and of the underlying assumptions. In order to determine the model appropriate to describe biological species, let us begin with conventional reaction-diffusion equation

$$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + au^k(1-u) - \sigma u, \qquad (3.3)$$

where $u(x,t)$ is the population density distribution, the diffusion term describes its variation, the next two terms characterize population natality and mortality rates. For $k = 1$ (asexual reproduction) and $a > \sigma$, this is the KPP-Fisher equation first suggested to describe propagation of dominant gene. Travelling waves described by this equation correspond to population spreading in some habitat but not to stationary density distributions in the genotype space. For $k = 2$ (sexual reproduction) such localized stationary distributions (positive solutions with zero limits at infinity - pulses) can exist, but they are unstable. Thus, this equation is not appropriate to describe biological species considered as localized stable stationary distributions of population density.

Instead of equation (3.3) we consider the following equation

$$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + au^2(1-I(u)) - \sigma u, \qquad (3.4)$$

where the integral

$$I(u) = \int_{-\infty}^{\infty} u(x,t)dx$$

is proportional to the total population density. In this case, the reproduction term takes into account global competition for resources independent of the phenotype $x$. It is different in comparison with equation (3.1) where the integral with kernel $\phi(x-y)$ describes nonlocal competition for resources which depends on the distance in the phenotype space.

The corresponding stationary problem

$$Dw'' + aw^2(1-I(w)) - \sigma w = 0, \quad w(\pm\infty) = 0$$

has positive solutions [12]. Their stability is not proved but it is verified in numerical simulations. Thus, equation (3.4) can have stable pulse solutions which can be interpreted as persistent biological species without bifurcations (lines E, F in Fig. 10).

Bifurcation of pulses from a homogeneous in space stationary solution can be studied for equation (3.4) considered in a bounded interval $[0, L]$ with no-flux boundary conditions. Linearizing the problem about this solution $w_*$, we obtain the eigenvalue problem

$$Dv'' + \sigma v - aw_*^2 I(w) = \lambda w, \quad v'(0) = v'(L) = 0.$$

We consider the eigenfunctions

$$v_n(x) = \cos(n\pi x/L), \quad n = 0, 1, 2, \dots$$

and determine the corresponding eigenvalues

$$\lambda_0 = \sigma - aw_*^2 L, \quad \lambda_n = -D(n\pi/L)^2 + \sigma, \quad n = 0, 1, 2, \dots$$

It can be verified that $\lambda_0 < 0$, while for $D$ sufficiently small, $\lambda_1 > 0$. Therefore, the loss of stability occurs with the eigenfunction $v_1(x)$ leading to the emergence of space-dependent solutions corresponding to pulses in the infinite domain.

This mechanism of the loss of stability of the homogeneous in space stationary solutions of nonlocal reaction-diffusion equations is generic. It is different in comparison with the Turing instability, and it can be realized in different biological and biomedical problems.

### 3.3. Species as distributions in the genotype space

According to the discussion in the previous section, biological species, considered as stable stationary pulses, exist due to the intraspecific competition. This conclusion may be realistic for some species (or quasi-species, see Section 4.2) but it can be too
limiting for some other species. For example, human height distribution (phenotypical parameter) represents a typical bell-shaped function close to the normal distribution. It can be considered as a stable stationary solution if we neglect its slow evolution in a long time scale. However, it would not be justified to assume that human shape distribution has a bell-shaped form due the competition for resources. Hence, this model has its limitations, and the underlying assumptions may not be always justified.

One of the main assumptions of the model is related to the diffusion term describing small random variations from parents to offsprings. It seems natural to assume that offsprings resembles their parents (this is described by the reproduction with the same x) and that small variations indeed occur. If we denote by $x_i$ some given phenotype and by $x_{i±1}$ some close variations, while $u_i$, $u_{i±1}$ are the corresponding population densities, then variations $u_i \Leftrightarrow u_{i±1}$ can be described by the discrete equations

$$ \frac{\partial u_i}{\partial t} = \mu(u_{i+1} - 2u_i + u_{i-1}) + au_i^k(1-u_i) - \sigma u_i, \quad (3.5) $$

where parameter $\mu$ characterizes the frequency of such variations. This equation represents a discrete version of equation (3.3) (or of equation (3.1) in the case of nonlocal consumption of resources).

It is implied in the derivation of equation (3.5) that direct and inverse transitions $u_i \Leftrightarrow u_{i±1}$ have the same rates. This assumption seems quite natural, but it may not be satisfied if we consider population density distributions as a function of their genotype defined as a sequence of loci which can be either in the original state or in the mutated state. Suppose that a sequence $A_k B_{n-k}$ contains $k$ original loci A and $n-k$ mutated loci B. Then a single mutation $A \to B$ converts it in a sequence $A_{k-1} B_{n-k+1}$ with probability $k\nu$, where $\nu$ is the probability of a single mutation and $k$ is the number of loci A where it can occur. However, inverse transition from a sequence $A_{k-1} B_{n-k+1}$ to a sequence $A_k B_{n-k}$ due to mutation $B \to A$ can occur with probability $(n-k+1)\nu$ proportional to the number of loci B. Thus, these transitions have different probabilities, and the corresponding equation becomes, instead of (3.3), as follows [71]:

$$ \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + \frac{\partial(xu)}{\partial x} + au^k(1-u) - \sigma u. \quad (3.6) $$

The distance $x$ in the genotype space is proportional to the number of mutations.

The first-order term, which appears in the right-hand side of this equation, essentially changes the properties of solutions. Instead of travelling waves described by equation (3.3), equation (3.6) describes stationary pulses, that is, positive solutions converging to zero at infinity. An example of such solution is given by normal distribution $u_n(x) = e^{-x^2/D}$ which represents a solution of this equation for $a = \sigma = 0$. If $a > 0$ and $k = 1$, then it is proved that this equation has a family of stationary pulses with polynomial decay $u_p(x) \sim |x|^{-1-a}$ at infinity and a single minimal pulse with fast decay rate $u_0(x) \sim x^a e^{-x^2/D}$ [71]. All these solutions are stable with respect to certain classes of initial conditions. In the case of the minimal solution $u_0(x)$, the class of initial conditions includes all functions with bounded supports that are appropriate in applications.

Thus, population density distributions with respect to their genotype converge to a positive stationary bell-shape function close to the normal distribution but with a slower decay rate at infinity (up to a polynomial factor). In order to make a link between the genotype (sequences of loci or letters) and phenotype (e.g., human height) we assume that the change of phenotype is proportional to the number of mutations, that is, to the distance in the genotype space (space variable in equation (3.6)). Then we conclude that observed population distributions (phenotype) converge to a normal-like distribution. Existence and stability of such distributions are determined by asymmetry of mutations: mutations in the direction of the mean-value phenotype are more frequent since there more sites where such mutations can occur.

Let us finally note that there is some similarity between the waves described by equation (3.3) in the monostable case and the pulses described by equation (3.6). In both cases, there is a family of such solutions; the minimal solution has a faster decay rate, and it determines the convergence of solution of the Cauchy problem for the initial conditions with bounded supports.

# 4. Biomedical applications

Reaction-diffusion waves are widely used in biomedical applications such as heart modelling and excitable media, tumor growth models or calcium waves. In this section we will discuss some other recent applications of reaction-diffusion waves related to propagation of viral infection and inflammation in cell culture or tissue. Nonlocal reaction-diffusion models in biology are reviewed in [13].

## 4.1. Viral infections

Viral infection spreads in cell culture or tissue due to virus replication inside infected cells and its transportation between cells, either as random motion in the extracellular matrix or direct cell-to-cell transmission. The simplest model describing this process includes the equations for the concentration of uninfected cells U, infected cells W and virus V:

$$ \frac{\partial U}{\partial t} = -aUV, \quad (4.1) $$

$$ \frac{\partial W}{\partial t} = aUV - \sigma_1 W, \quad (4.2) $$

$$ \frac{\partial V}{\partial t} = D \frac{\partial^2 V}{\partial x^2} + bW_{\tau} - \sigma_2 V. \quad (4.3) $$
The right-hand side of equation (4.1) describes infection rate of uninfected cells, a similar term with the opposite sign appears in the next equation followed by the death term of infected cells. The right-hand side in equation (4.3) describes virus random motion, its replication and death. The replication rate is proportional to the concentration of infected cells with time delay, $W_τ(x, t) = W(x, t - τ)$ taking into account replication delay $τ$.

System (4.1)-(4.3), considered on the whole axis, can have a travelling wave solution $U(x, t) = u(x - ct)$, $W(x, t) = w(x - ct)$, $V(x, t) = v(x - ct)$ with the limits at infinity

$$x = -\infty : U = U_f, W = 0, V = 0, \quad x = \infty : U = U_0, W = 0, V = 0. \qquad (4.4)$$

The initial concentration of uninfected cells $U_0$ is a given constant, their final concentration $U_f$ is unknown and should be found from the solution of the problem. It can be shown [6] that the ratio $\eta = U_f/U_0$ satisfies the equation

$$\ln \eta = R_v(1 - \eta), \qquad (4.5)$$

where $R_v = abU_0/(\sigma_1\sigma_2)$ is the virus replication number. This equation has a solution in the interval $0 < \eta < 1$ (as it should be according to its biological meaning) if and only if $R_v > 1$. Hence, propagation of infection waves is observed if the virus production rate exceeds its elimination. It is interesting to note that equation (4.5) is similar to the equation for the basic reproduction number in the epidemiological models.

If $R_v > 1$, then infection spreading can be characterized by two main quantities, the speed of propagation and the viral load. The speed of infection wave correlates with virus virulence, and it can be analytically determined by the linearization method [4,6]:

$$c = \min_{\mu > \mu_0} \frac{D\mu\sqrt{\mu + \sigma_1}}{\sqrt{(\mu + \sigma_1)(\mu + \sigma_2) - abU_0e^{-\mu\tau}}}.$$

Here $\mu_0$ is a positive zero of the denominator which exists due to the condition $R_v > 1$. Let us note that this formula gives the minimal speed of propagation, which is more appropriate for the applications. We can conclude from the formula for the wave speed that increase of virus mutation and production rate (D, a, b) augments virus virulence, while time delay in its replication decreases it.

Viral load, defined as $V_L = \int_{-\infty}^{\infty} v(\xi)d\xi$, can also be determined analytically [6]:

$$V_L = \frac{c}{a} \ln \eta \approx \frac{cR_v}{a}.$$

The last approximation follows from (4.5) for large values of $R_v$, for which $\eta \ll 1$.

In the case of respiratory viral infections, viral load in the upper respiratory tract determines infection transmission rate. Together with the wave speed and the virus replication number, these are the main characterizations of viral infections in cell culture. They can be used to compare different viruses or strains, such as Delta and Omicron variants of SARS-CoV-2 infection [8], or, in a more general context virus competition and mutations [5,92]. In the case of viral infection in the human organism, we need to take into account the immune response [7], and other factors [16,17].

## 4.2. Virus mutations and quasi-species

Ensemble of viruses with close genotypes are called virus quasi-species. Similar to biological species, they can be considered as density distributions with respect to their genotypes. Then we can analyse evolution of the existing viruses and the emergence of new ones. Assuming for simplicity that the concentration of uninfected cells in the previous model (4.1)-(4.3) is constant (their depletion can be compensated by the influx of new cells), we introduce the immune response [18]:

$$\frac{\partial W}{\partial t} = aU_0V - (\sigma_1 + \gamma_1 J(C))W, \qquad (4.6)$$

$$\frac{\partial V}{\partial t} = D \frac{\partial^2 V}{\partial x^2} + bW - (\sigma_2 + \gamma_2 J(C))V, \qquad (4.7)$$

$$\frac{\partial C}{\partial t} = kV - \sigma_3 C. \qquad (4.8)$$

Here $C(x, t)$ is the concentration of immune cells which depends on the genotype variable $x$ and time $t$. In the case of adaptive immune response, antigen-specific T-lymphocytes eliminate infected cells, and antigen specific B-lymphocytes produce virus neutralizing antibodies. For simplicity of presentation, we do not distinguish different types of lymphocytes in the model and consider them as a single concentration. Immune cells eliminate infected cells and virus in some range of genotypes implying that the corresponding terms in equations (4.6) and (4.7) include the integral

$$J(C)(x, t) = \int_{-\infty}^{\infty} \phi(x - y)C(y, t)dy,$$

where the kernel $\phi(x - y)$ determines the cross-reactivity of the immune response, $\int_{-\infty}^{\infty} \phi(x)dx = 1$. The right-hand side of the last equation of this system describes production of the immune cells stimulated by antigens and their death or depletion.
<figure><img src="image_15.png" /><figcaption>Fig. 11. The concentration $V(x,t)$ in numerical simulations of system (4.6)-(4.8). Each curve corresponds to some fixed moment of time. Time increase is shown by colour change from blue to red. The initial condition is a step-wise constant function in the middle of the interval. We observe periodic wave propagation from the unstable zero point to the spatial structure resulting from the instability of the positive stationary point. (Figure courtesy of D. Neverova.)</figcaption></figure>

This problem has an infection-free stationary solution $W = V = C = 0$ and, if $R_v > 1$, a positive stationary solution $(W_0, V_0, C_0)$. It is stable with respect to the corresponding ODE system. However, in the presence of the space variable, it can lose its stability resulting in the emergence of spatial structures [18]. The stability boundary ($\lambda = 0$) is given by the equality

$$D\xi^2 = -\frac{k\gamma_2 V_0}{\sigma_3} \left( 1 + \frac{b\gamma_1 W_0^2}{aU_0\gamma_2 V_0^2} \right) \tilde{\phi}(\xi), \quad (4.9)$$

where $\tilde{\phi}(\xi)$ is the Fourier transform of the function $\phi(x)$, $\xi \in \mathbb{R}$. In the case of a step-wise constant cross-reactivity function,

$$\phi(x) = \frac{1}{2N} \begin{cases} 1, & |x| \le N \\ 0, & |x| > N \end{cases},$$

which means that lymphocytes eliminate cells and viruses with sufficiently close antigen-specificity, we have $\tilde{\phi}(\xi) = \sin(\xi N)/(\xi N)$. For $N$ sufficiently small, equation (4.9) with respect to $\xi$ does not have solution, and the instability does not occur. If $N$ is large enough, then this equation has some solutions, and the homogeneous in space stationary solution $(W_0, V_0, C_0)$ loses its stability leading to the emergence of spatial structures with the corresponding spatial frequency $\xi$. An example of numerical simulations is shown in Fig. 11 where a spatial structure emerges behind the propagating periodic wave.

The peaks of virus density distribution in the genotype space correspond to different quasi-species which emerge due to the immune escape. This effect is qualitatively similar to the emergence of biological species due to the intra-specific competition previously studied in the nonlocal reaction-diffusion equation [100] (see Section 3). Some other models were developed to study the existence of virus quasi-species [54,62], their dynamics and emergence to the anti-viral treatment [19–21].

### 4.3. Inflammation

Inflammation is a part of the immune response to pathogens or injury. It can be considered as a driving force of inflammatory diseases, such as atherosclerosis, rheumatoid arthritis or asthma. But it is also strongly involved in many other diseases including cancer, Alzheimer disease or viral infections. Let us begin with a relatively simple model of inflammation in order to present the underlying biological processes and their description in the model. We consider the following system of equations for the concentration of uninflamed cells $U$, inflamed cells $I$, immune cells $M$, and inflammatory cytokines $C$:

$$\frac{\partial U}{\partial t} = -\alpha U C, \quad (4.10)$$

$$\frac{\partial I}{\partial t} = \alpha U C - \beta I, \quad (4.11)$$

$$\frac{\partial M}{\partial t} = \delta I + kC - \gamma_1 M, \quad (4.12)$$

$$\frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} + b_1 I + b_2 M - \gamma_2 C. \quad (4.13)$$
Inflammatory cytokines act on uninflamed cells and lead to their activation converting them into inflamed cells. The corresponding terms are present in the right-hand side of the first two equations. Proliferation (or influx) and death of uninflamed cells are neglected in the model, while death of inflamed cells is described by the second term in the right-hand side of equation (4.11).

Immune cells (macrophages, lymphocytes) are attracted to the inflammation site by the inflammatory cytokines and by inflamed cells. For example, in the case of atherosclerosis inflamed (activated) endothelial cells at the surface of blood vessel walls determine penetration of monocytes from blood flow inside vessel wall. The influx of immune cells is described by the first two terms in the right-hand side of equation (4.12) and their death by the last term. The last equation for the concentration of inflammatory cytokines describes their diffusion in the tissue, their production by inflamed cells and immune cells and their degradation or depletion.

In order to describe inflammation propagation as a reaction-diffusion wave, we set $U(x,t) = u(x-ct)$, $I(x,t) = v(x-ct)$, $M(x,t) = w(x-ct)$, $C(x,t) = z(x-ct)$, where c is the wave speed. These functions satisfy the following system of equations:

$$cu' - \alpha u z = 0, \qquad (4.14)$$

$$cv' + \alpha u z - \beta v = 0, \qquad (4.15)$$

$$cw' + \delta v + kz - \gamma_1 w = 0, \qquad (4.16)$$

$$Dz'' + cz' + b_1 v + b_2 w - \gamma_2 z = 0, \qquad (4.17)$$

and have the following limits at infinity:

$$x = -\infty : u = u_f, v = 0, w = 0, z = 0, x = \infty : u = u_0, v = 0, w = 0, z = 0. \qquad (4.18)$$

The initial concentration of uninflamed cells $u_0$ is given, the final concentration $u_f$ is unknown. It can be determined from the equation

$$R(1 - \omega) = \ln \omega, \qquad (4.19)$$

where $\omega = u_f/u_0$, and parameter

$$R = \frac{\alpha u_0 (b_1 \gamma_1 + b_2 \delta)}{\beta (\gamma_1 \gamma_2 - b_2 k)}$$

is called inflammation number [33]. This equation is similar to equation (4.5) for viral infection where virus replication number $R_v$ is replaced by the inflammation number $R$.

There are three different regimes of inflammation progression depending on the value of inflammation number. If $R > 1$, then equation (4.19) has a solution $0 < \omega < 1$. In this case, inflammation spreads in the tissue as a reaction-diffusion wave with the speed

$$c = \min_{\mu > \mu_0} \frac{\sqrt{D\mu \sqrt{(\mu + \beta)(\mu + \gamma_1)}}}{\sqrt{(\mu + \beta)(\mu + \gamma_1)(\mu + \gamma_2)} - \alpha b_1 u_0 (\mu + \gamma_1) - b_2 k (\mu + \beta) - \alpha \delta b_2 u_0},$$

where $\mu_0 > 0$ is the value for which the denominator vanishes, and which exists for $R > 1$.

If $0 < R < 1$, then there is no inflammation wave, and solution converges to the inflammation-free equilibrium. Condition $R > 0$ in this and in the previous case means that the total death rate of inflammatory factors (immune cells, cytokines) exceeds their production rate. If this condition is not satisfied and $R < 0$, then instead of normal inflammation presented above, the model describes self-amplifying growth of inflammatory factors and excessive inflammation which can occur during cytokine storm and other abnormal immune and inflammatory reactions.

Thus, the model considered here describes different regimes of inflammation progression depending on the inflammation number. We neglect in this model, among other simplifying assumptions, the turnover of uninflamed cells. However, in a longer time scale, their proliferation (or influx) and death should be considered, either for the description of normal physiological conditions or wound healing [46].

**Disease-specific models of inflammation.** The models discussed above can be considered as generic models of inflammation since they are based on some general properties of the immune response, and they do not take into account some possible disease-specific features of these processes. For example, in the case of atherosclerosis, such specific features include cholesterol and low-density lipoproteins (LDL) crucial for the disease initiation and progression, formation of atheroma plaque and its interaction with blood flow. There is a vast literature devoted to atherosclerosis modelling (see [34] and the references therein). Inflammation plays a crucial role in many other diseases including viral infections [63], Alzheimer disease [30] or cancer [53,108]. In many cases, inflammation spreads as a reaction-diffusion wave. Modelling and analysis of such solutions can help in the understanding of the disease progression and in the optimization of medical treatment.

## 5. Concluding remarks

Reaction-diffusion waves have been one of the most fruitful concepts in nonlinear systems. Being originally discovered in 1937 in the biological context of a new gene spread and shortly followed by their independent discovery in chemical physics, the original works by Fisher [36], Kolmogorov et al. [51] and Zeldovich & Frank-Kamenetskii [104] were eventually followed by hundreds of
new studies concerned both with mathematical aspects of reaction-diffusion waves (e.g. see [35,99,105] and further references there) and numerous applications.

In spite of the large amount of research done so far about reaction-diffusion waves properties and applications (in particular in biology, e.g. see [100]), it remains to be a very active field of study. We mention that most of the recent studies on reaction-diffusion waves are essentially interdisciplinary, spanning across all natural sciences, medicine and even social science and economy. A growing list of new applications brings new challenges for the correspondingly mathematical theory. Arguably, some of the results reviewed in this paper are worthy of a further rigorous mathematical investigation. One such result is the existence of transient travelling waves. Classical theory of reaction-diffusion systems usually treats the travelling wave as asymptotic dynamics, so that, once it has emerged, it will persist indefinitely.<fnref n="3" /> However, some recent results [29,73] suggest that a travelling wave can as well be a manifestation of long transient dynamics [44,60,61]. In such a case, a travelling wave only exists for a finite (but very long) time and then disappears without any apparent reason, e.g. see Fig. 7.

We further recall here that reaction-diffusion systems are a paradigmatic example of a nonlinear system with a capability of self-organisation [42,58], i.e., when the properties of the dynamics emerging in the course of time depend little on the initial conditions (that are 'forgotten' after a short relaxation time) but instead are determined by the processes, in particular by the system's local kinetics and/or peculiarities of the spatial transport. The existence of transient reaction-diffusion waves (along with similar observations on transient pattern formation, e.g. see [61]) therefore suggests that a 'transient self-organisation' may in some cases be more relevant than the classical theory. Further development of the corresponding mathematical theory of such 'transient self-organisation' could become a direction of exciting new research.

Declaration of competing interest

The authors declare no conflict of interest.

Acknowledgements

Both V.V. and S.P. were supported by the RUDN University Strategic Academic Leadership Program.

References

[1] Alharbi W, Petrovskii SV. Patterns of invasive species spread in a landscape with a complex geometry. Ecol Complex 2018;33:93-105.

[2] Alhassan Y, Siekmann I, Petrovskii S. Mathematical model of oxygen minimum zones in the vertical distribution of oxygen in the ocean. Sci Rep 2024;14:22248.

[3] Alhassan Y, Petrovskii S. Standing on a cliff: OMZ growth may indicate the approach of a global anoxia. Mathematics 2024. Submitted for publication.

[4] Ait Mahiout L, Bessonov N, Kazmierczak B, Sadaka G, Volpert V. Infection spreading in cell culture as a reaction-diffusion wave. ESAIM: Math Model Numer Anal 2022;56(3):791-814.

[5] Ait Mahiout L, Kazmierczak B, Volpert V. Viral infection spreading and mutation in cell culture. Mathematics 2022;10(2):256.

[6] Ait Mahiout L, Mozokhina A, Tokarev A, Volpert V. Virus replication and competition in a cell culture: application to the SARS-CoV-2 variants. Appl Math Lett 2022;133:108217.

[7] Ait Mahiout L, Mozokhina A, Tokarev A, Volpert V. The influence of immune response on spreading of viral infection. Lobachevskii J Math 2022;43(10):2699-713.

[8] Ait Mahiout L, Bessonov N, Kazmierczak B, Volpert V. Mathematical modeling of respiratory viral infection and applications to SARS-CoV-2 progression. Math Methods Appl Sci 2023;46(2):1740-51.

[9] Alnufaie A. Mathematical modelling of population dynamics of social protests. PhD thesis. Leicester, UK: University of Leicester, 2023.

[10] Alsulami A, Glukhov A, Shishlenin M, Petrovskii S. Dynamical modelling of street protests using the Yellow Vest Movement and Khabarovsk as case studies. Sci Rep 2022;12:20447.

[11] Banerjee M, Petrovskii SV, Volpert V. Nonlocal reaction-diffusion models of heterogeneous wealth distribution. Mathematics 2021;9:351.

[12] Banerjee M, Vougalter V, Volpert V. Doubly nonlocal reaction-diffusion equations and the emergence of species. Appl Math Model 2017;42:511-99.

[13] Banerjee M, Kuznetsov M, Udovenko O, Volpert V. Nonlocal reaction-diffusion equations in biomedical applications. Acta Biotheor 2022;70(2):12.

[14] Berestycki H, Hamel F. Generalized travelling waves for reaction-diffusion equations. In: Perspectives in nonlinear partial differential equations. In honor of H. Brezis. Contemp math, vol. 446. Amer. Math. Soc.; 2007. p. 101-23.

[15] Berestycki H, Nadal J-P, Rodriguez N. A model of riots dynamics: shocks, diffusion and thresholds. Netw Heterog Media 2015;10:443-75.

[16] Bessonov N, Bocharov G, Mozokhina A, Volpert V. Viral infection spreading in cell culture with intracellular regulation. Mathematics 2023;11(6):1526.

[17] Bessonov N, Volpert V. Airway obstruction in respiratory viral infections due to impaired mucociliary clearance. Int J Numer Methods Biomed Eng 2023:e3707.

[18] Bessonov N, Neverova D, Popov V, Volpert V. Emergence and competition of virus variants in respiratory viral infections. Front Immunol 2023;13:945228.

[19] Bessonov N, Bocharov G, Volpert V. Space and genotype-dependent virus distribution during infection progression. Mathematics 2021;10(1):96.

[20] Bessonov N, Bocharov G, Meyerhans A, Popov V, Volpert V. Existence and dynamics of strains in a nonlocal reaction-diffusion model of viral evolution. SIAM J Appl Math 2021;81(1):107-28.

[21] Bessonov N, Bocharov GA, Leon C, Popov V, Volpert V. Genotype-dependent virus distribution and competition of virus strains. Math Mech Complex Syst 2020;8(2):101-26.

[22] Bessonov N, Reinberg N, Volpert V. Mathematics of Darwin's diagram. Math Model Nat Phenom 2014;9(3):5-25.

[23] Bessonov N, Reinberg N, Banerjee M, Volpert V. The origin of species by means of mathematical modelling. Acta Bioteoretica; 2018. p. 1-12.

[24] Biggs M. Size matters: quantifying protest by counting participants. Sociol Methods Res 2018;47(3):351-83.

[25] Bonnasse-Gahot L, Berestycki H, Depuiset M. Epidemiological modelling of the 2005 French riots: a spreading wave and the role of contagion. Sci Rep 2018;8:107.

[26] Breitburg D, et al. Declining oxygen in the global ocean and coastal waters. Science 2018;359:eaam7240.

[27] Brockmann D, Hufnagel L, Geisel T. The scaling laws of human travel. Nature 2006;439:462-5.

<footnote n="3">³ In an unbounded space; in a finite domain, it persists until it reaches the domain boundary.</footnote>
[28] Cartes C, Asahi K, Fernández R. Riots and subways, a relationship moderated by the neighborhood's income level. Sci Rep 2022;12:10557.

[29] Chowdhury PR, Banerjee M, Petrovskii S. A two-timescale model of plankton-oxygen dynamics predicts formation of oxygen minimum zones and global anoxia. J Math Biol 2024;89(1):8.

[30] Ciuperca I, Pujo-Menjouet L, Matar-Tine L, Torres N, Volpert V. A qualitative analysis of an Aβ-monomer model with inflammation processes for Alzheimer's disease. R Soc Open Sci 2024;11(5):231536.

[31] Davies TP, Fry HM, Wilson AG, Bishop SR. A mathematical model of the London riots and their policing. Sci Rep 2013;3:1303.

[32] Edwards AM, Brindley J. Zooplankton mortality and the dynamical behaviour of plankton population models. Bull Math Biol 1999;61:303-39.

[33] El Hajj W, El Khatib N, Volpert V. Inflammation propagation modeled as a reaction-diffusion wave. Math Biosci 2023;365:109074.

[34] El Khatib N, Kafi O, Sequeira A, Simakov S, Vassilevski Y, Volpert V. Mathematical modelling of atherosclerosis. Math Model Nat Phenom 2019;14(6):603.

[35] Fife P. Mathematical aspects of reacting and diffusing systems. New York: Springer; 1979.

[36] Fisher RA. The wave of advance of advantageous genes. Ann Eugen 1937;7:355-69.

[37] Gause GF. The struggle for existence. Baltimore: Williams and Wilkins; 1934.

[38] Girardin L, Nadin G. Travelling waves for diffusive and strongly competitive systems: relative motility and invasion speed. Eur J Appl Math 2015;26:521-34.

[39] Gökce A, Yazar S, Sekerci Y. Delay induced nonlinear dynamics of oxygen-plankton interactions. Chaos Solitons Fractals 2020;141:110327.

[40] Gorodetskii V, et al. Coupling between adjacent crystal planes in heterogeneous catalysis by propagating reaction-diffusion waves. Nature 1994;370:276-9.

[41] Granovetter M. Threshold models of collective behavior. Am J Sociol 1978;83:1420-43.

[42] Haken H. Synergetics. An introduction. Springer series in synergetics, vol. 1. Berlin: Springer-Verlag; 1977.

[43] Hardin G. The competitive exclusion principle. Science 1960;131(3409):1292-7.

[44] Hastings A, Abbott KC, Cuddington K, Francis T, Gellner G, Lai YC, et al. Transient phenomena in ecology. Science 2018;361:eaat6412.

[45] Holt RD. Species coexistence. Reference module in life sciences. Elsevier; 2017.

[46] Hussain S, Sen M, Volpert V. Infection spreading in tissue as a reaction-diffusion wave. Math Med Biol, J IMA 2024;41(3):169-91.

[47] Kan-On Y. Parameter dependence of propagation speed of travelling waves for competition-diffusion equations. SIAM J Math Anal 1995;26(2):340-63.

[48] Kan-On Y. A note on the propagation speed of travelling waves for a Lotka-Volterra competition model with diffusion. J Math Anal Appl 1998;217(2):693-700.

[49] Kinezaki N, Kawasaki K, Takasu F, higesada N. Modeling biological invasions into periodically fragmented environments. Theor Popul Biol 2003;64:291-302.

[50] Kierstead H, Slobodkin LB. The size of water masses containing plankton blooms. J Mar Res 1953;12:141-7.

[51] Kolmogorov AN, Petrovskii IG, Piskunov NS. A study of the diffusion equation with increase in the quantity of matter, and its application to a biological problem. Bull Moscow Univ Math Ser A 1937;1:1-25.

[52] Kondo S, Miura T. Reaction-diffusion model as a framework for understanding biological pattern formation. Science 2010;329:1616-20.

[53] Kuznetsov M, Clairambault J, Volpert V. Perspectives in cancer treatment. Reply to comments on "Improving cancer treatments via dynamical biophysical models". Phys Life Rev 2022.

[54] Leon C, Kutsenko I, Volpert V. Existence of solutions for a nonlocal reaction-diffusion equation in biomedical applications. Isr J Math 2022;248(1):67-93.

[55] Lin J, Andreasen V, Casagrandi R, Levin SA. Traveling waves in a model of influenza A drift. J Theor Biol 2003;222:437-45.

[56] Lewis MA, Kareiva P. Allee dynamics and the spread of invading organisms. Theor Popul Biol 1993;43:141-78.

[57] Lewis MA, Petrovskii SV, Potts J. The mathematics behind biological invasions. Interdisciplinary applied mathematics, vol. 44. Springer; 2016. p. 362p.

[58] Mikhailov AS, Loskutov AY, Hermann Haken H. Foundations of synergetics I: distributed active systems. Springer; 2012.

[59] Mimura M, Kawasaki K. Spatial segregation in competitive interaction-diffusion equations. J Math Biol 1980;9:49-64.

[60] Morozov A, Abbott KC, Cuddington K, Francis T, Gellner G, Hastings A, et al. Long transients in ecology: theory and applications. Phys Life Rev 2020;32:1-40.

[61] Morozov A, Feudel U, Hastings A, Abbott KC, Cuddington K, Petrovskii SV. Long-living transients in ecological models: recent progress and open questions. Phys Life Rev 2024. https://doi.org/10.1016/j.plrev.2024.11.004.

[62] Moussaoui A, Volpert V. The influence of immune cells on the existence of virus quasi-species. Math Biosci Eng 2023;20(9):15936-56.

[63] Mozokhina A, Ait Mahiout L, Volpert V. Modeling of viral infection with inflammation. Mathematics 2023;11(19):4095.

[64] Murray JD. Mathematical biology. Heidelberg: Springer; 1989.

[65] Namba T, Mimura M. Spatial distribution of competing populations. J Theor Biol 1980;87:795-814.

[66] Okubo A, Maini PK, Williamson MH, Murray JD. On the spatial spread of the grey squirrel in Britain. Proc R Soc B 1989;238:113-25.

[67] O'Boyle S. Oxygen depletion in coastal waters and the open ocean. CRC/Taylor and Francis; 2020.

[68] Ozarslan R, Sekerci Y. Fractional order oxygen-plankton system under climate change. Chaos 2020;30:033131.

[69] Paulmier A, Ruiz-Pino D. Oxygen minimum zones (OMZs) in the modern ocean. Prog Oceanogr 2009;80:113-28.

[70] Pena MA, Katsev S, Oguz T, Gilbert D. Modeling dissolved oxygen dynamics and hypoxia. Biogeosciences 2010;7:933-57.

[71] Pena B, Trofimchuk S, Volpert V. Biological populations as stationary distributions in the space of genotypes. J Differ Equ 2023;359:303-29.

[72] Petrovskii SV. Localization of a nonlinear switching wave in an active medium with an isolated inhomogeneity. Tech Phys 1997;42:866-71.

[73] Petrovskii SV, Medvinsky AB, Ivanitskii GR. On the effect of a dynamical "confinement" of a localized initial perturbation in a predator-prey type system. Trans (Dokl) Russ Acad Sci 1997;357:550-3.

[74] Petrovskii SV, Malchow H. Critical phenomena in plankton communities: KISS model revisited. Nonlinear Anal, Real World Appl 2000;1:37-51.

[75] Petrovskii SV, Kawasaki K, Takasu F, Shigesada N. Diffusive waves, dynamical stabilization and spatio-temporal chaos in a community of three competitive species. Jpn J Ind Appl Math 2001;18:459-81.

[76] Petrovskii S, Alharbi W, Alhomairi A, Morozov A. Modelling population dynamics of social protests in time and space: the reaction-diffusion approach. Mathematics 2020;8:78.

[77] Petrovskii SV, Sekerci Y, Venturino E. Regime shifts and ecological catastrophes in a model of plankton-oxygen dynamics under the climate change. J Theor Biol 2017;424:91-109.

[78] Pitcher GC, et al. System controls of coastal and open ocean oxygen depletion. Prog Oceanogr 2021;197:102613.

[79] Pizarro-Koch M, et al. On the interpretation of changes in the subtropical oxygen minimum zone volume off Chile during two La Niña events (2001 and 2007). Front Mar Sci 2023;10:1155932.

[80] Potts JR, Petrovskii SV. Fortune favours the brave: movement responses shape demographic dynamics in strongly competing populations. J Theor Biol 2017;420:190-9.

[81] Protter MH, Weinberger HF. Maximum principles in differential equations. New York: Springer; 1984.

[82] Sekerci Y, Petrovskii SV. Mathematical modelling of plankton-oxygen dynamics under the climate change. Bull Math Biol 2015;77:2325-53.

[83] Sekerci Y, Petrovskii SV. Global warming can lead to depletion of oxygen by disrupting phytoplankton photosynthesis: a mathematical modelling approach. Geosciences 2018;8:201.

[84] Sekerci Y, Petrovskii SV. Pattern formation in a model oxygen-plankton system. Computation 2018;6:59.

[85] Sekerci Y, Ozarslan R. Oxygen-plankton model under the effect of global warming with nonsingular fractional order. Chaos Solitons Fractals 2020;132:109532.

[86] Sherratt JA, Murray JD. Mathematical analysis of a basic model for epidermal wound-healing. J Math Biol 1991;29(5):389-404.

[87] Shigesada N, Kawasaki K, Teramoto E. Spatial segregation of interacting species. J Theor Biol 1979;79:83-99.

[88] Shigesada N, Kawasaki L, Teramoto E. Traveling periodic waves in heterogeneous environments. Theor Popul Biol 1986;30:143-60.

[89] Shigesada N, Kawasaki K. Biological invasions: theory and practice. Oxford: Oxford University Press; 1997.
[90] Skellam JG. Random dispersal in theoretical populations. Biometrika 1951;38(1-2):196-218.

[91] Tilman D. Resource competition and community structure. Princeton, NJ: Princeton University Press; 1982.

[92] Tokarev A, Mozokhina A, Volpert V. Competition of SARS-CoV-2 variants in cell culture and tissue: wins the fastest viral autowave. Vaccines 2022;10(7):995.

[93] Toth A, Horvath D, Siska A. Velocity of propagation in reaction-diffusion fronts of the chlorite-tetrathionate reaction. J Chem Soc Faraday Trans 1997;93(1):73-6.

[94] Trofimchuk S, Volpert V. Traveling waves in delayed reaction-diffusion equations in biology. Math Biosci Eng 2020;17:6487-514.

[95] Vakulenko S, Volpert V. Generalized travelling waves for perturbed monotone reaction-diffusion systems. Nonlinear Anal, Theory Meth Appl 2001;46(6):757-76.

[96] van Panhuis WG, Choisy M, Xiong X, et al. Region-wide synchrony and traveling waves of dengue across eight countries in Southeast Asia. Proc Natl Acad Sci USA 2015;112:13069-74.

[97] Vicker MG. Eukaryotic cell locomotion depends on the propagation of self-organized reaction-diffusion waves and oscillations of actin filament assembly. Exp Cell Res 2002;275(1):54-66.

[98] Volpert AI, Khudyaev SI. Analysis in classes of discontinuous functions and equations of mathematical physics. Dordrecht: Nijoff; 1985.

[99] Volpert A, Volpert Vit, Volpert VI. Traveling wave solutions of parabolic systems. Providence: AMS; 1994.

[100] Volpert V, Petrovskii S. Reaction-diffusion waves in biology. Phys Life Rev 2009;6:267-310.

[101] Volpert V. Mathematical model of repressive response to collective action and protest waves. J Theor Biol 2024;595:111970.

[102] Watson AJ, Lenton TM, Mills B. Ocean deoxygenation, the global phosphorus cycle and the possibility of human-caused large-scale ocean anoxia. Philos Trans A 2017;375:20160318.

[103] Webb P. Introduction to oceanography. Roger Williams University; 2021.

[104] Zeldovich YB, Frank-Kamenetskii DA. Theory of uniform flame propagation. Zh Fiz Khim 1938;12:100.

[105] Zeldovich YB, Barenblatt GI, Librovich VB, Makhviladze GM. The mathematical theory of combustion and explosions. New York: Consultants Bureau; 1985.

[106] Zhang LH. How do synaptic coupling and spatial temporal delay influence traveling waves in nonlinear nonlocal neuronal networks? SIAM J Appl Dyn Syst 2007;6:597-644.

[107] Zhao H, Zhu L. Dynamic analysis of a reaction-diffusion rumor propagation model. Int J Bifurc Chaos 2016;26:1650101.

[108] Zhao Huakan, Wu Lei, Yan Guifang, Chen Yu, Zhou Mingyue, Wu Yongzhong, et al. Inflammation and tumor progression: signaling pathways and targeted intervention. Signal Transduct Targeted Ther 2021;6:263.