01. Theory

Mathematical Background

In applied mathematics, partial differential equations (PDEs) are used to describe the development of physical phenomena. Frequently these PDEs may incorporate a small parameter, that scale some of the terms in the equation. If these terms are important then the presence of this small parameter may be critical; setting the parameter to zero can fundamentally alter the nature and behavior of the equations. A notable example is the Reynolds number in the Navier-Stokes equations, which balances viscous and momentum terms. When the Reynolds number is small, viscous forces predominate, allowing the momentum terms to be neglected. Conversely, when the Reynolds number is large, momentum forces dominate. However, it is not possible to completely neglect the viscous terms, as high Reynolds number flows near boundaries create very thin, highly viscous regions known as boundary layers. Such equations are called singularly perturbed equations. For a simple example of this phenomena applied to an ordinary differential equation (ODE), see example 1 in Section 9.

In this exemplar, we will develop a partial differential equation (PDE) solver to address a class of singularly perturbed equations commonly encountered in mathematical biology, specifically reaction-diffusion systems. These systems describe interactions between two competing agents, such as predator and prey or different chemical substances. When one agent diffuses more slowly than the other, the system exemplifies a singularly perturbed equation.

Reaction-diffusion systems are renowned for their ability to produce a wide variety of spatial patterns, known as Turing patterns. These patterns are hypothesized to be responsible for natural phenomena such as the stripes on zebras, the spots on leopards, and the complex skin patterns of giant pufferfish.

General Model:

Let

$$u := u(t,x,y) \quad \text{and}\quad v := v(t,x,y)$$

represent agent densities which vary in both time $t$ and spatial coordinates $x$ and $y$. Reaction-diffusion systems can then be expressed in the form:

$$\frac{\partial u }{\partial t} = \epsilon \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial x^2} \right) + F(u,v), \quad \frac{\partial v }{\partial t}= \left(\frac{\partial^2 v}{\partial x^2}+ \frac{\partial^2 v}{\partial y^2} \right) + G(u,v),$$

where $\epsilon >0$ is some small parameter and $F(u,v)$ and $G(u,v)$ are given functions. $F$ and $G$ represent reaction terms that may be non-linear.

For the spatial boundary conditions, we will consider either Dirichlet boundary conditions or Neumann boundary conditions.

We will additionally specify initial conditions of the system at the start of the simulation, to be evolved over time.

Predator-Prey Model:

In a predator-prey models, $u$ represents the prey density and $v$ the predator density. We can imagine the prey $u$ as a quantity of deer wandering around an enclosed area, while the predators $v$ are wolves that limit the deer population. Then $F$ and $G$ may be of the form

$$F(u,v) = u(1-u) - v\frac{u}{1+u},$$

$$G(u,v) = v\frac{u}{1+u}.$$

In this context, the first term of $F$, $u(1-u)$, represents the prey density increasing with logistic growth. In contrast, the second term in $F$ (and $G$ itself) represent the predator's response to the prey population, with the understanding that the predator can only consume a limited amount of prey (see type II response). The small parameter $\epsilon$ represents a slowly dispersing prey relative to the predators.

Here we show a solution to this equation in one dimension where we have the initial conditions

$$u(t=0,x) = e^{-x^2}\cos^2(4\pi x),$$

$$v(t=0,x) = e^{(x+0.8)^2}(0.2 + 0.1\sin(10\pi x))^2.$$

We additionally applied Dirichlet boundary conditions at $x=0$ and $x=1$ and considered a time domain of $[0,0.5]$. With $\epsilon = 0.05$, we discretised the system with $6000$ points in time and $400$ points in space.

Notice how the shape of the prey density is retained for a much longer time than that of the predator density. This is due to the small parameter $\epsilon$.

Activator-Substrate Depletion Model:

The Activator-Substrate Depletion Model is a framework used to describe a specific type of oscillatory behavior observed in biochemical systems. - Activator $u$: A molecule or substance that promotes or accelerates a particular biochemical reaction or process. - Substrate $v$: The substance upon which the activator acts, typically being converted into a product by the action of the activator. Then $F$ and $G$ may be of the form

$$F(u,v) = u^2v - u,$$

$$G(u,v) = \mu(1 - u^2v).$$

The constant $\mu$ qualitatively describes the relative rates of chemical reactions for substrate compared to activator. Small $\epsilon$ models a scenario when the activator disperses slowly relative to the substrate's dispersion. When $\epsilon$ is small and $\mu$ greater than $1$, Turing patterns may emerge in the solution.

Here we have an activator - substrate example where the initial conditions have been given as $u(0,x,y)=1+r$, where $r$ is a small amount of random noise and $v(0,x,y) = 1$. Dirichlet conditions were used at the boundaries. Notice how, despite looking similar, the activator $u$ has more defined structures than $v$. $v$ has more diffusion between the structures. This is due to $\epsilon$ being small.

See Example 4 in Section 9 for more details.

Problem Complexities

This system of PDEs is fairly complicated:

• Coupled Equations: This system of coupled equations involves both $u$ and $v$, with each equation's solution depending on the other. This interdependence complicates the solution process, as the equations cannot be solved in isolation.
• Nonlinearity: The functions $F(u,v)$ and $G(u,v)$ are typically non-linear, i.e. of the form $u^2+uv$ or $v\frac{u}{1+u}$. Non-linear PDEs are much harder to solve both analytically and numerically. We will tackle nonlinearities numerically with Newton iteration.
• Spatial and Temporal Dependencies: The variables $u$ and $v$ depend on both time $t$ and spacial coordinates $(x,y)$. The equations are first order in time $\frac{\partial u }{\partial t}$, which we can handle numerically with temporal marching.
• Order of Equations: The equations are also second order in space $\frac{\partial^2 u }{\partial x^2}+\frac{\partial^2 u }{\partial y^2}$. This second order derivative means we must solve these equations at every point in space concurrently. This second derivative quantity, called the Laplacian of $u$, measures diffusion, or how $u$ changes around a point relative to the points neighbours.
• Solution Method: To solve this equation numerically we employ a spatial discretisation scheme such as the finite difference scheme. See this exemplar for a introduction to finite differences.
• Small parameter $\epsilon$: As discussed above, the small parameter $\epsilon$ indicates that the diffusion of $u$ is much smaller than the diffusion of $v$. When $\epsilon$ is smaller, we will see less diffusion in $u$: any structures of $u$ will appear sharper than in $v$.
• Boundary Layers: If the behavior of the equation results in the formation of a boundary layer, this is handled numerically using non-uniform grids. These grids are denser in regions where rapid changes occur, such as near the boundary layer, and coarser elsewhere. To implement this numerically, the non-uniform grids are mapped to a uniform computational grid, which is then utilised for computing finite differences. Note that in both the examples above, stretched grids are not necessary as the sharpness of the diffusion of $u$ appears uniformly throughout the domain. See example 1 in Section 9 to see the effect of the stretched grid.
• Boundary conditions: At a particular boundary $(x_0,y)$ of the domain (this corresponds to an $x$ boundary which is defined for all points in $y$) a boundary condition must be specified for both $u$ and $v$. Two common types of boundary condition are:
  • A Dirichlet boundary condition is one of the form $u(t,x=x_0,y) = L$, where $L$ is a constant.
  • A Neumann boundary condition is one of the form $\frac{\partial u}{\partial x} (t,x=x_0,y) = L$, where $L$ is a constant.
• Initial Conditions: We must specify initial condition of the form $u(t=t_0,x,y) = f(x,y)$, where $f(x,y)$ is a given function. If we are solving for both $u$ and $v$ then initial conditions will be given for both.