# Modelling formation of stationary periodic patterns in growing population of
motile bacteria

Valentina Bucur and Bakhtier Vasiev  Department of Mathematical Sciences,
Mathematical Sciences Building, University of Liverpool, Liverpool L69 7ZL, UK

###### Abstract

Biological pattern formation is one of the most intriguing phenomena in
nature. Simplest examples of such patterns are represented by travelling waves
and stationary periodic patterns which occur during various biological
processes including morphogenesis and population dynamics. Formation of these
patterns in populations of motile microorganisms such as Dictiosteium
dictiostelium and E. coli have been shown in a number of experimental studies.
Conditions for formation of various types of patterns are commonly addressed
in mathematical studies of dynamical systems containing diffusive and
advection terms. In this work, we do mathematical study of spatio-temporal
patterns forming in growing population of chemotactically active bacteria. In
particular, we perform linear analysis to find conditions for formation of
stationary periodic patterns, and nonlinear (Fourier) analysis to find
characteristics, such as amplitude and wavelength, of these patterns. We
verify our analytical results by means of numerical simulations.

##  1 Introduction

Pattern formation is one of the outcomes of biological self-organization and
manifested by well known patterns, such as animal markings, and processes,
such as embryogenesis. Growing bacterial communities form patterns which, for
example, have an impact on biofilm formation [16]. Pattern formation is a
phenomenon which is widely observed in nature. Examples of non-biological
pattern formation include formation of dunes in deserts, temperature hot spots
in autocatalytic chemical systems [12] and many other processes in various
physical [27], chemical [33] and social [24] systems.

Turing, in his fundamental work [28], has pointed to two most representative
examples of patterns, namely travelling waves and stationary periodic
patterns. Formation of travelling waves in biological populations has been
studied for many decades by means of Kolmogorov-Fisher equation [9, 20] and
various models based on this equation [1, 8]. One of the main futures of such
models is an involvement of advection term which takes into account the
directional motion of biological species in response to chemical signals
(chemotaxis) [3, 10, 15, 17, 31].

The most important outcome of Turing’s work, presented in [28], was finding
mechanism of stationary periodic pattern formation in chemical systems.
Formation of this pattern was explained by so-called diffusion driven
instability which nowadays is commonly referred as ”Turing instability”.
Turing has looked at a system of differential equations describing the
dynamics of two interacting substances, or morphogens. Commonly, in a well-
mixed system (or in the absence of diffusion) the system has a stable
homogeneous solution. However, for certain types of interactions between the
morphogens, i.e. in so-called ”activator-inhibitor” system, homogeneous state
becomes unstable when the diffusion of one of morphogens (inhibitor) is larger
as compared with the diffusion of the other (activator). In this case,
formation of non-homogeneous patterns such as stripes or spots takes place [5,
18, 23]. This mechanism of pattern formation in the case of activator-
inhibitor systems is commonly phrased as ’local self-enhancement and long-
range inhibition’ [11, 32].

Stationary periodic patterns are often observed in populations of
chemotactically active microorganisms, which aggregate and form regularly
spaced clusters. Chemotactic activity is expressed as motion along or against
gradient of chemicals, produced externally or by microorganisms themselves.
Many theoretical studies have addressed the mechanisms of pattern formation in
growing populations of such microorganisms represented by bacteria or amoeba
[18, 13]. The general form of such models is given as

$\begin{cases}\displaystyle\frac{\partial n}{\partial
t}=D_{n}\nabla^{2}n-\nabla\left(\chi(n,c)\nabla c\right)+f(n,c),\\\\[10.0pt]
\displaystyle\frac{\partial c}{\partial t}=D_{c}\nabla^{2}c+g(n,c),\\\
\end{cases}$ (1)

where the variable $n$ represents the density of microorganisms, and the
variable $c$ \- the concentration of chemotactic agent. Parameter $D_{n}$ is a
diffusion coefficient describing the random motion of cells, $D_{c}$ \- the
diffusion constant of chemotactic agent, function $\chi(n,c)$ describes the
chemotactic response of the microorganisms, and functions $f(n,c)$ and
$g(n,c)$ describe the microbial and chemical kinetics, respectively. If
$\chi(n,c)=0$ (no chemotaxis) the system (1) represents a reaction-diffusion
model of the type of classical Lotka-Volterra model [22, 23]. If
$\chi(n,c)\neq 0$ then we have an advection term in the first equation and the
system (1) is a reaction-diffusion-advection model. Furthermore, if
$f(n,c)=0$, the system transfers into the well-known Keller-Segel model [18].

Keller-Segel model was presented in 1970 and designed to describe the
aggregation on Dictyosteium discoideum amoeba mediated by chemotactic agent
(cAMP, which at that time was known as acrasin) produced by the amoeba
themselves [18]. The original model consisted of four equation, but using
quasi-steady-state assumption (i.e. by replacing ’fast’ differential equations
by algebraic relationships) it was reduced to a system of two reaction-
diffusion-advection equations. The analysis was performed for the simplest
version of the model where it was assumed that the chemotactic sensitivity is
constant ($\chi(n,c)=\chi_{0}$), amoeba do not die and reproduce,
($f(n,c)=0$), and the chemical agent is produced by the amoeba and degrade at
a constant rate, $g(n,c)=hn-pc$, where $h$ is the production rate per amoebae
and $p$ is the degradation rate. Under these assumptions, the model (1) is
given by the following system:

$\begin{cases}\displaystyle\frac{\partial n}{\partial
t}=D\nabla^{2}n-\chi_{0}\nabla^{2}c,\\\\[5.0pt] \displaystyle\frac{\partial
c}{\partial t}=\nabla^{2}c+hn-pc,\\\ \end{cases}$

where $D=D_{n}/D_{c}$ is the ratio of diffusion constants of the amoeba and
chemical and the spatial variable is rescaled in a way that the diffusion of
the chemical $c$ is set to one. In the absence of chemotaxis ($\chi=0$) the
system evolves towards the stable steady state $(n,c)=(n_{0},c_{0}$) defined
by the initial amount of amoeba, so that $n_{0}$ is an average initial density
of amoeba and $c_{0}=n_{0}*h/p$. It was shown in [18] that in the presence of
chemotaxis the homogeneous state ($n_{0},c_{0}$) becomes unstable if

$\displaystyle\frac{\chi_{0}p}{Dh}+\frac{n_{0}h}{p}>1.$

This condition implies that an increase in the chemotactic activity $\chi$
breaks down the stability of homogeneous state and the amoeba aggregate into
clusters forming stationary periodic patterns. An assumption that the
chemotactic response doesn’t depend on the density of cells and on the
concentration of chemotactic agent (that is, $\chi(n,c)=\chi_{0}$) is
oversimplification which was made in [18] to make the mathematical analysis
easier. It is reasonable to assume that the chemotactic response is
proportional to the density of cells, and there is experimental evidence
indicating that it is reduced with an increase of the concentration of
chemotactic agent. Our analysis of the modification of model, where the
chemotactic term was assumed having the form
$\displaystyle\chi(n,c)=\chi_{0}\frac{n}{c}$, have shown that the aggregation
takes place if  $\displaystyle\frac{\chi_{0}}{D}+\frac{n_{0}h}{p}>1$. This
result doesn’t alter the main conclusion drawn from the previous version of
the model that the chemotactic response should be large enough to allow
formation of stationary periodic patterns. Other modifications of Keller-Segel
model (involving different forms of chemotactic term and/or bacterial and
chemical kinetics) have been considered by a number of researchers and some on
them are listed in [13].

E. coli bacteria is another biological species which exhibits chemotactic
activity. Patterns formed in E. coli colonies have been studied experimentally
in great details [4, 5]. It was shown that the waves of expansion of growing
E. coli colonies transform into stationary spots formed by clusters of
dead/quiescent bacteria. It was also shown that these patterns form because E.
coli respond to chemoattractant which they secrete themselves. It was
established that the observed patterns of cell aggregation depends on nutrient
consumption, cell proliferation, excretion of attractant and chemotactic
motility [4]. A number of models have been developed to describe patterns
formed in growing E. coli colonies [2, 25]. The earliest models reproducing
the formation of spots in E. coli colonies were presented in [29]. One of
these models (so-called ’liquid model’ which is designed to reproduce patterns
formed by bacteria in liquid environment) is very similar to Keller-Segel
model and given by the following two equations:

$\begin{cases}\displaystyle\frac{\partial n}{\partial
t}=D\nabla^{2}n-\chi_{0}\nabla\left(\frac{n}{(1+c)^{2}}\nabla
c\right),\\\\[7.0pt] \displaystyle\frac{\partial c}{\partial
t}=\nabla^{2}c+\omega\frac{n^{2}}{\nu+n^{2}}.\end{cases}$

In this model the death and proliferation of cells is again neglected
($f(n,c)=0$). The chemotactic sensitivity is given as
$\displaystyle\chi(n,c)=\chi_{0}\frac{n}{1+c^{2}}$, that is, it is assumed to
be proportional to the density of cells and the dependence on concentration of
the chemotactic agent is modified in a way that the absence of chemical
doesn’t imply the division by zero. The decay of the chemical is neglected
while its production is given by the Hill function $\displaystyle
g(n,c)=\omega\frac{n^{2}}{\nu+n^{2}}$. Analysis of this model performed in
[29] has shown that for sufficiently large $\chi_{0}$ the homogeneous state
becomes unstable with formation of periodic patterns for which, however, the
spacial periodicity changes over on time. Thus, the patterns obtained using
models presented in [29] are not stationary. Formation of stationary periodic
patterns was shown later by other researchers in a framework of models for
dynamics of E. coli population represented by larger systems of differential
equations [2, 25].

While mechanisms and conditions of periodic pattern formation have been
addressed in many studies [28, 18, 22] their properties, such as wavelength
and amplitude, where studied in a much lesser extend. The estimation of
wavelength of patterns is commonly performed on the basis of linear
instability analysis and considered to be given by the wave number of most
unstable mode (see Chapter 2 of [23]). Most recently properties of Turing
patterns formed in a class of two-variable reaction-diffusion models, with
kinetics terms given by polynomial functions, were addressed in [6]. Amplitude
of the patterns were estimated in assumption that they have a sinusoidal form,
that is, that the solution of (1) (in one-dimensional case) is given as:

$\begin{cases}n(x,t)=n_{0}+\alpha\cos(kx),\\\
c(x,t)=c_{0}+\beta\cos(kx),\end{cases}$

where ($n_{0}$, $c_{0}$) represents homogeneous steady state of the system,
$k$ is the wave number of most unstable mode found from linear stability
analysis and $\alpha$ and $\beta$ give amplitude of Turing pattern emerging in
the model. The estimate of the amplitude, obtained this way, has matched
reasonably well numerical results presented in this research. However, it is
evident that more accurate estimates of the wavelength and amplitude of Turing
patterns should be based on nonlinear analysis of these patterns.

In this paper we study properties of periodic patterns forming in a
representative version of model (1) which is relatively simple, justifiable
biologically, and allows such pattern formation. The choice of the model was
made on the basis of analysis of patterns forming in a range of models of type
(1) with non-zero advection term. We use Fourier series for analysis of
patterns obtained in numerical simulations as well as for analytical
estimation of the amplitude and wavelength of Turing patterns emerging in the
model. For greater accuracy we perform our analytical studies for a small
domain which can contain only half (or one) spike. We also study how the
amplitude and wavelength of Turing patterns depend on the values of model
parameters.

##  2 Model

In this section, we introduce the model of growing population of motile
bacteria which we will analyse throughout this paper. This model is
represented by a system of two partial-differential equations which include
diffusion, reaction and advection terms. It is known that such systems allow
formation of travelling wavefronts [10, 21] as well as Turing patterns [13,
22, 23].

We look for Turing patterns forming in one-dimensional domain (of size $L$)
under no-flux boundary conditions. We have performed preliminary numerical
simulations using various models of type (1) and based on the results of these
simulations we have chosen the following system:

$\begin{cases}\displaystyle\frac{\partial n}{\partial
t}=D_{n}\frac{\partial^{2}n}{\partial
x^{2}}-\tilde{\chi}\frac{\partial}{\partial x}\left(n\frac{\partial
c}{\partial x}\right)+r_{0}n\left(1-\frac{n}{k}\right),\\\\[10.0pt]
\displaystyle\frac{\partial c}{\partial t}=D_{c}\frac{\partial^{2}c}{\partial
x^{2}}+hn-pc,\\\ \end{cases}$

This system is a version of the system (1) and describes the dynamics of
population of bacteria, whose density is given by $n$, which respond
chemotactically to a chemical agent of concentration $c$, produced by bacteria
themselves. Here the chemotactic sensitivity of cells is assumed to be
proportional to $n$ and takes the form $\chi(n,c)=\tilde{\chi}n$. Also, cells
reproduce according to logistic law with growth rate $r_{0}$ and carrying
capacity $k$: $f(n,c)=r(1-r/k)$. The kinetics of chemical is assumed to be
linear , $g(n,c)=hn-pc$ (like in Keller-Segel model [18]), that is, it is
being produced by cells with rate $h$ and decays with rate $p$. This system
can be nondimensionalised using the substitutions $\tilde{t}=tp$,
$\displaystyle\tilde{x}=x\sqrt{\frac{p}{D_{c}}}$,
$\displaystyle\tilde{n}=\frac{n}{k}$, $\displaystyle\tilde{c}=\frac{cp}{hk}$.
This gives (after dropping the tildes):

$\begin{cases}\displaystyle\frac{\partial n}{\partial
t}=D\frac{\partial^{2}n}{\partial x^{2}}-\chi_{0}\frac{\partial}{\partial
x}\left(n\frac{\partial c}{\partial x}\right)+rn(1-n),\\\\[10.0pt]
\displaystyle\frac{\partial c}{\partial t}=\frac{\partial^{2}c}{\partial
x^{2}}+n-c,\\\ \end{cases}$ (2)

with new constants $\displaystyle D=\frac{D_{n}}{D_{c}}$,
$\displaystyle\chi_{0}=\frac{\tilde{\chi}hk}{D_{c}p}$ and $\displaystyle
r=\frac{r_{0}k}{p}$ While the parameters $D$ and $r$ are essentially non-
negative, $\chi$ can be positive (chemoattraction) or negative
(chemorepulsion).

In the absence of diffusion and chemotaxis, the system is said to be well-
mixed, meaning that both variables, $n$ and $c$, don’t depend on spatial
variable $x$. In this case, the first equation is detached from the second and
describes the logistic growth of population of bacteria. The system has two
steady states $(n,c)=(0,0)$ and $(1,1)$, which are unstable and stable
respectively [22]. After putting back the diffusion term ($D\neq 0$,
$\chi_{0}=0$) the first equation transfers into Kolmogorov-Fisher equation
which is known for travelling front solutions [20, 9]. The system now allows
formation of travelling wavefronts between the steady states $(n,c)=(0,0)$ and
$(n,c)=(1,1)$, which are now stable and unstable, respectively [22].

When the system includes the flow of cells due to chemotaxis ($\chi_{0}\neq
0$) then the two equations in (2) are coupled and the known solutions of this
system include stationary periodic, or Turing, patterns [22]. When the
chemotaxis is weak ($\chi_{0}$ is small), we have a model, which, similarly to
the Fisher-Kolmogorov equation, allows only formation of travelling wavefront
(see Fig.1A). However, if the chemotactic flow is strong enough, that is, the
value of $\chi_{0}$ is above certain value, then the formation of stationary
periodic patterns given by periodic functions, $n(x)$ and $c(x)$, can be
observed (see Fig.1B,C).

Figure 1: Numerical simulation of travelling front and Turing patterns forming
in the system (2). Profiles of cell density, $n$, (blue line) and
concentration of chemotactic agent, $c$, (red line) are shown. Solutions are
obtained in the domain $0\leq x\leq 200$, $t\geq 0$. Initial conditions:
$n(x,0)=0$ if $x>2$ and $n(x,0)=1$ if $0\leq x\leq 2$; $c(x)=0$ $\forall x$.
A: formation of travelling front in a system when the chemotaxis is weak
($\chi_{0}=1.7$). B, C: formation of Turing pattern in a system when the
chemotaxis is strong ($\chi_{0}=1.9$). Other model parameters: $D=1$, $r=0.1$.
Profiles are shown at time $T=200$ in panels A and B, and $T=800$ in panel C.

##  3 Conditions for formation of Turing patterns

Our numerical simulations, illustrated in Fig.1, have indicated that Turing
patterns can appear only if the chemotactic sensitivity, $\chi_{0}$, is above
some threshold value. This threshold value can be found analytically by
applying the standard linear analysis technique, which is well known from the
literature [18, 23]. It is known that Turing patterns form due to so-called
Turing instability when a steady state, which is stable in the well-mixed
system, becomes unstable in the full reaction-diffusion-advection system. This
means that the equilibrium point $(n_{0},c_{0})=(1,1)$ in system (2), which is
stable if $D=\chi_{0}=0$, should become unstable for some non-zero values of
$D$ and $\chi_{0}$.

Investigation of stability of steady state $(n_{0},c_{0})$ involves
linearisation of the system (2) at this state and checking the dynamics of
perturbations, which are given in the form $e^{\lambda t}\cos(kx)$. Parameter
$k$ represents a perturbation mode and can have any value from the set, given
as $k=i\pi/L$, with $i$ taking any integer value so that the perturbation
satisfies no-flux boundary conditions. It appears that these perturbations
represent solutions of the linearised system provided that $\lambda$ is an
eigenvalue of the characteristic matrix, which is, for the system (2) at the
steady state $(1,1)$, given as

$M=\begin{pmatrix}-Dk^{2}-r&\chi_{0}k^{2}\\\ 1&-k^{2}-1\end{pmatrix}.$ (3)

The steady state $(1,1)$ is stable if both eigenvalues of this matrix are
negative, which is known to require $\det M>0$ and $\operatorname{Tr}M<0$
[23]. We can see that $\operatorname{Tr}M=-Dk^{2}-r-k^{2}-1$ is always
negative, so for the system to be driven unstable by perturbation, we require
$\det M<0$, which leads to the following condition for Turing instability,

$\det M=Dk^{4}+k^{2}(D+r-\chi_{0})+r<0.$

From this, a necessary but not sufficient condition is

$D+r-\chi_{0}<0$ (4)

since if $D+r-\chi_{0}>0$ then $\det M>0$ for any $k$.The second necessary
condition for initiation of Turing patterns comes from the fact that the
equation $\det M=0$ should have real roots, so that $\det M<0$ between these
two roots. Clearly, the roots are given as

$\displaystyle
k^{2}_{1,2}=\displaystyle\frac{-(D+r-\chi_{0})\pm\sqrt{(D+r-\chi_{0})^{2}-4Dr}}{2D},$
(5)

and therefore we require

$(D+r-\chi_{0})^{2}-4Dr>0.$

Taking into account (4) we get

$\chi_{0}-D-r>2\sqrt{Dr}.$

After rearranging terms in this inequality and, by introducing new notation
$R_{T}$, we can state the condition for Turing instability in the form:

$R_{T}=\frac{\chi_{0}}{D+r+2\sqrt{Dr}}>1.$ (6)

Looking again at Figure 1, we see that in panel A, the factor of instability
is $R_{T}=0.98<1$, while in panels B and C, $R_{T}=1.1>1$. Thus, results shown
in Fig.1 confirm that if $R_{T}<1$ there observed travelling waves rather than
Turing pattern, while if $R_{T}>1$ formation of Turing patterns is observed.

Using the above technique we have performed stability analysis for a number of
models, which are represented by system (1) but differ by functions
$\chi(n,c)$ and $f(n,c)$. Following [18] we assumed that the chemical kinetics
is linear and didn’t alter it across models, i.e. $g(n,c)=n-c$ in all models.
List of the analysed models together with the obtained formulas for
instability factor, $R_{T}$, is shown in Table 1. Note that travelling wave
solutions form if $R_{T}<1$ (and $f(n,c)\neq 0$), while if $R_{T}>1$ then
formation of Turing patterns should be observed.

| $\chi(n,c)$ | $f(n,c)$ | $R_{T}$
---|---|---|---
$M_{1}$ | $\chi_{0}$ | 0 | $\displaystyle\frac{\chi_{0}}{D}$
$M_{2}$ | $\chi_{0}n$ | 0 | $\displaystyle\frac{\chi_{0}n_{0}}{D}$
$M_{3}$ | $\chi_{0}n$ | $rn(1-n)$ | $\displaystyle\frac{\chi_{0}}{D+r+2\sqrt{Dr}}$
$M_{4}$ | $\displaystyle\chi_{0}\frac{1}{c}$ | 0 | $\displaystyle\frac{\chi_{0}}{Dn_{0}}$
$M_{5}$ | $\displaystyle\chi_{0}\frac{1}{c}$ | $rn(1-n)$ | $\displaystyle\frac{\chi_{0}}{D+r+2\sqrt{Dr}}$
$M_{6}$ | $\displaystyle\chi_{0}\frac{n}{c}$ | $rn(1-n)$ | $\displaystyle\frac{\chi_{0}}{D+r+2\sqrt{Dr}}$
$M_{7}$ | $\displaystyle\chi_{0}\frac{n}{c+\nu}$ | $rn(1-n)$ | $\displaystyle\frac{\chi_{0}}{(\nu+1)(D+r+2\sqrt{Dr})}$
$M_{8}$ | $\displaystyle\chi_{0}\frac{n}{(c+\nu)^{2}}$ | $rn(1-n)$ | $\displaystyle\frac{\chi_{0}}{(\nu+1)^{2}(D+r+2\sqrt{Dr})}$
$M_{9}$ | $\displaystyle\chi_{0}\frac{n+\nu_{1}}{c+\nu_{2}}$ | $rn(1-n)$ | $\displaystyle\frac{\chi_{0}(1+\nu_{1})}{(\nu_{2}+1)(D+r+2\sqrt{Dr})}$
Table 1: Results of Turing instability analysis for nine models represented by
system (1). Model notations are given in the first column. Functions
$\chi(n,c)$ and $f(n,c)$ corresponding to each model are given in columns 2
and 3. Calculated formula for instability factor, $R_{T}$, is given in column
4. In models M2 and M4, where $f(n,c)=0$, the instability factor, $R_{T}$,
depends on $n_{0}$, which is an average value of $n$ in the domain and
determined by the initial conditions.

Most of the models listed in Table 1 have been introduced earlier (see [13]),
however, the conditions for Turing instability for most of these models have
not been stated so far. These conditions are qualitatively similar for all
considered models and this is illustrated by domains where these conditions
satisfied (in parameter planes for models $M3$ and $M7$) in Fig. 2. We see
from this figure that for Turing instability to take place the chemotactic
sensitivity, $\chi_{0}$ should be strong enough, while the diffusion of
chemotactic agent, $D$ and reproduction rate $r$ are small enough.

Figure 2:  Domains $R_{T}>1$ (where Turing instability takes place) in
parameter planes for models $\mathbf{M3}$ (solid) and $\mathbf{M7}$ (dashed).
Plots $R_{T}=1$ are drown according to the formulas for $R_{T}$ given in Table
1 in an assumption that one of parameters has a fixed (default) value, which
are $D=1$, $\chi_{0}=1.9$, $r=0.1$. Dashed line is plotted for $\nu=0.2$
(model $M7$). The dot indicates position of the point with coordinates given
by default values of parameters.

We have reproduced numerically Turing patterns forming in all models given in
Table 1 aiming to find the simplest model, which allows formation of Turing
patterns, such that its shape is within physical constrains. Here is a summary
of our numerical results:

  * •

$M_{1}$ is the simplest model where the chemotactic sensitivity is assumed to
be constant, and no cell death/proliferation takes place [18]. Numerical
simulations show that Turing patterns in this model evolve towards infinite
amplitude and involve negative values for density of cells and concentration
of chemical (see Fig.3A), which is physically impossible. So, we can conclude,
that $\chi(n,c)$ has to depend on at least one of the two variables, $n$ and
$c$.

  * •

In $M_{2}$ the chemotactic sensitivity is assumed to be proportional to the
concentration of cells, $\chi(n,c)=\chi_{0}n$, but there is still no cell
death/proliferation are taking place. When the instability condition is met we
observe only one sharp spike forming in response to stimulus (see Fig. 3B).
Thus, no Turing pattern can be observed in this model.

  * •

$M_{3}$ is $M_{2}$ with added logistic growth of cells, $f(n,c)=rn(1-n)$. In
this model Turing patterns of finite amplitude and non-negative values of $n$
and $c$ are observed. Thus, this is a simplest model which allows reproduction
of Turing pattern of physically justifiable shape.

  * •

In $M_{4}$ the chemotactic sensitivity is inversely proportional to the
concentration of cells and there is no cell proliferation. Here numerical
simulations end up with division by zero and abruptly stop after short period
of time.

  * •

$M_{5}$ is $M_{4}$ with added logistic growth of cells. Simulations can be ran
for longer periods of time, however division by zero is still taking place.

  * •

In $M_{6}$ $\chi(n,c)$ is proportional to $n$ and inversely proportional to
$c$. Eventually, division by zero is not taking place here and we obtain
appropriately shaped Turing patterns.

  * •

The main feature of models $M_{7,8,9}$ is that $\chi(n,c)$ is proportional to
linear function of $n$ and inverse proportional to linear function, or a
square of linear function of $c$. The latter property lets us to safely avoid
division by zero (see [13, 29]). In all these cases, we get reasonably shaped
Turing patterns, which can be sharper or smoother for particular models (i.e.
$M_{7}$, $M_{8}$ or $M_{9}$).

Figure 3: Numerical simulations of patterns forming in models $M_{1}$ (panel
A), $M_{2}$ (panel B) and $M_{3}$ (panel C). Profiles of cell density, $n$,
(blue line) and concentration of chemotactic agent, $c$, (red line) at time
$T=500$ in the domain of size $L=50$ shown in all three panels. Initial
condition $n=c=0$ in (A) and (C), and $n=1$ and $c=0$ in (B). Patterns were
initiated using stimulus applied to the centre of the domain: $n=1.1$ for
$24<x<26$ at time $t=0$. Model parameters: $D=1$, $\chi_{0}=1.8$ and $r=0.1$
(in $M_{3}$).

Our main conclusion after analysis and numerical simulations presented in this
chapter that model $M_{3}$ is the simplest model where periodic stationary
patterns, which do not violate any physical constrains, are observed. Model
$M_{3}$ (which is given by the system (2)) will be used in the rest of this
paper for analysis of properties of Turing patterns such as their amplitude
and wavelength.

##  4 Estimation of the wavelength of Turing patterns

In the previous sections, we have found conditions for formation of Turing
patterns (given, for example, by Eq. (4) for the model (2)) and have presented
numerical simulations (see Fig.1) to illustrate the validity of this
estimation. Our next task is to analyse properties of Turing patterns, which
are given by their wavelength and amplitude. The amplitude is defined as the
difference between maximal and minimal values of cell density, $n$, in the
pattern, while the wavelength - as a distance representing spatial periodicity
of Turing pattern. As an alternative to the wavelength we will also consider a
characteristic length which is a distance between the consecutive maxima and
minima and twice smaller than the wavelength.

In this section we consider the wavelength of Turing pattern and will estimate
it on the basis of linear stability analysis presented in the previous
section. Turing pattern is a periodic pattern, and as a such can be
represented as a sum of cosinusoidal profiles, each given as $a(k)\cos(kx)$,
where $k=i\pi/L$, $i$ is an integer number and $L$ is the domain size. Integer
$i$ should be such that $k$ belongs to the interval $k\in(k_{1},k_{2})$ where
$k_{1}$ and $k_{2}$ are given by Eq. (5). Wavenumbers $k$ in this interval
correspond to unstable modes, each having the wavelength $\Lambda=\pi/k$. It
is reasonable to assume that Turing pattern is predominantly defined by the
most unstable mode. Popular way to choose this mode is to take the one
corresponding to the average of $k_{1}^{2}$ and $k_{2}^{2}$ given by Eq.(5)
[23, 6], which is:

$\displaystyle k^{2}_{av}=\displaystyle\frac{\chi_{0}-D-r}{2D}.$

The wavelength corresponding to this mode is given as:

$\displaystyle\Lambda=\pi\sqrt{\frac{2D}{\chi_{0}-D-r}}.$ (7)

More accurate way to estimate the most unstable mode is to find a mode for
which the characteristic matrix (3) has the maximal positive eigenvalue. In
terms of trace and determinant, the two eigenvalues of the characteristic
matrix, $M$, are

$\displaystyle\lambda_{1,2}=\frac{trM\pm\sqrt{trM^{2}-4\det M}}{2},$

which gives

$\displaystyle\lambda_{1,2}=\frac{-Dk^{2}-k^{2}-r-1}{2}\pm\frac{\sqrt{k^{4}(D^{2}-2D+1)+k^{2}(2D(r-1)+4\chi_{0}-2r+2)+1-2r}}{2}.$

Clearly, $\lambda_{1}<0$ always, so we are interested in the value of $k^{2}$
which corresponds to the largest value of $\lambda_{2}$. This means that we
are looking for the value of $k^{2}$, at which the derivative of $\lambda_{2}$
with respect to $k^{2}$ is zero:

$\displaystyle\frac{d\lambda_{2}}{dk^{2}}=\frac{-D-1}{2}+\frac{\tau(1-D)^{2}+(2\chi_{0}+Dr-
D-r+1)}{2\sqrt{\tau^{2}(1-D)^{2}+2\tau(2\chi_{0}+Dr-D-r+1)-2r+1}}=0.$ (8)

Using Matlab code, we have found the value of $k^{2}$ (for fixed values of
parameters $D$, $\chi_{0}$ and $r$) satisfying the above equation, and hence
calculated the corresponding wavelength, $\Lambda=\pi/k$.

Figure 4:  Dependence of the wavelength of Turing pattern, $\Lambda$, on
parameters of model (2) found according to the equation (7) (dashed) and the
implicit formula (8) (solid). Default values of model parameters: $D=1$,
$\chi_{0}=1.9$ and $r=0.1$.

To complete the presented linear analysis we have checked the dependence of
the wavelength, $\Lambda$, on the values of model parameters for the system
(2) according to formula (7) or implicit formula (8). The corresponding plots
are shown in Fig. 4. One can see that according to both formulas the
wavelength is monotonically increasing function of diffusion, $D$ (panel A),
and decreasing function of chemotactic sensitivity, $\chi_{0}$ (panel B).
Vertical asymptotes in panels A and B confirm that Turing patterns can’t be
observed, when the diffusion is too fast, or the chemotactic sensitivity is
too weak (see the condition given by Eq.(6) and domains for Turing instability
in Fig. 2). While we observe certain quantitative differences between plots in
Panels A and B, plots on Panel C differ qualitatively: the wavelength of most
unstable mode in Turing pattern is a decreasing function of the growth rate,
$r$, according to formula (8) and an increasing function according to the
formula (7). Our numerical simulations (shown in Section 7) indicate that the
dependence given by the formula (8) is not only correct qualitatively but also
quite accurate quantitatively.

We note that Equations (7) and (8) only give approximations of the wavelength
for Turing pattern based on the linear analysis of the system (2). These
approximations will later be compared with results obtained by means of
nonlinear analysis and numerical simulations.

##  5 Fourier series for numerically simulated profiles

Turing patterns are stationary periodic patterns and as such can naturally be
represented as Fourier series. Pattern forming in a domain of size $L$ under
no-flux boundary conditions is given as a superposition of cosines, that is in
the case of the system (2) we have:

$\begin{cases}\displaystyle n(x)=\sum_{i=0}^{M}\alpha_{i}\cos\frac{i\pi
x}{L},\\\\[10.0pt] \displaystyle c(x)=\sum_{i=0}^{M}\beta_{i}\cos\frac{i\pi
x}{L}.\end{cases}$ (9)

The coefficients $\alpha_{i}$ and $\beta_{i}$ define the amplitudes of mode
$i$ for the variables $n$ and $c$ (the amplitudes are given as $2\alpha_{i}$
and $2\beta_{i}$). For smooth profiles these coefficients quickly tend to zero
as $i$ increases and therefore we can truncate the series by taking the first
$M$ terms, with $M$ to be carefully defined. For known profile, $n(x)$, the
coefficients are found using the formulas

$\alpha_{0}=\frac{1}{L}\int_{0}^{L}n(x)dx,$

and for $i>0$

$\alpha_{i}=\frac{2}{L}\int_{0}^{L}n(x)\cos\frac{i\pi x}{L}dx.$ (10)

In the rest of this chapter we will focus on profiles $n(x)$ (and coefficients
$\alpha_{i}$) having in mind that the analysis of profile $c(x)$ is done in
the same way.

Formulas (10) can be used for spectral analysis of Turing patterns obtained
numerically. A typical stationary profile for cell density, $n$, which is
obtained in numerical simulations of the system (2) is shown in Fig.5A.
Spectral decomposition of this profile shows that only three modes have
reasonably high coefficients, namely, $\alpha_{0}=0.8979$,
$\alpha_{19}=0.3945$ and $\alpha_{38}=0.1293$, while all other coefficients
are considerably smaller (lesser than 0.1). $\alpha_{0}$ gives the average
level of $n$ for the entire pattern, $\alpha_{19}$ defines the amplitude of
the mode whose characteristic length is 1/19th of the domain size, that is the
one represented by 9.5 spikes. This is exact number of spikes we see in
Fig.5A. Finally, $\alpha_{38}$ corresponds to the second harmonic of the main
harmonic given by $\alpha_{19}$. Thus, the amplitude of the pattern shown in
Fig.5A can be estimated as $2\alpha_{19}$. Varying the size of the domain will
result to a change in the number of observed spikes, however their amplitude
and spatial periodicity wouldn’t change. This is illustrated by patterns shown
in panels B and C of Fig.5, which are obtained in simulations of small domains
where only half (B) and one (C) spike can fit in.

Figure 5: Profiles of cell density, $n(x)$, for Turing patterns in system (2)
obtained numerically. Domain size: $L=100$ in A, $L=5.5$ in B, and $L=11$ in
C. Values of model parameters $D=1$, $\chi_{0}=1.9$ and $r=0.1$.

Patterns occurring in a small domain (like those shown in Fig.5B and C) will
be of particular interest for us, as corresponding Fourier series can be
truncated at reasonably law value of $M$ (see eq. (9)). Using Fourier
decomposition of numerically obtained patterns we have investigated how the
coefficients of Fourier modes depend on the size of domain. In Fig.6 we have
shown plots of $\alpha_{i}$, for $i=0,1,...,10$, against the size of domain,
$L$, which varied from $0$ to $25$. We can see that while the value of
$\alpha_{0}$ does not change vary much (stays in the range between $0.9$ and
$1$), and $\alpha_{10}$ is always negligibly small, other coefficients vary
significantly. For $L<4$ $\alpha_{0}=1$ and all other coefficients are zero,
which means that the system is in homogeneous state and no spikes have formed.
For $4<L<7.5$ the first coefficient, $\alpha_{1}$, is larger than any
subsequent coefficient ($\alpha_{1}>\alpha_{i}$ for any $i>1$), and this
reflects the fact that only half of the spike can fit in. Within this range of
domain sizes the value of $\alpha_{1}$ increases from zero and then decreases
to zero reaching its maximum value $\alpha_{1}(=\alpha_{max})=0.385$ at
$L=5.5$. This domain size, which corresponds to the maximal value of the
coefficient $\alpha_{1}$ will be considered as a proper characteristic length
(or half-wavelength) of Turing pattern and will be denoted as $\Lambda_{0}$,
i.e. $\Lambda_{0}=5.5$ for the system (2) with model parameters in simulations
used to produce Fig. 6.

Figure 6: Dependence of coefficients $\alpha_{i}$, $i=0,1,...,10$ on domain
size, $L$. Plots are obtained by applying formula (10) to $n-$profiles from
numerical simulations. Model parameters: $D=1$, $\chi_{0}=1.9$ and $r=0.1$.
Vertical dashed lines indicate domain sizes given in whole number of
characteristic length $\Lambda_{0}$.

We see in Fig.6 that with the further increase of medium size ($L>7.5$) the
coefficient $\alpha_{2}$ increases and then decreases, and then all the same
with the coefficients $\alpha_{3}$ and $\alpha_{4}$, each reaching the maximum
of $\alpha_{max}=0.385$ when the size of the domain $L=i*\Lambda_{0}$ for
$\alpha_{i}$ (where $i=2$, $3$ and $4$), which are indicated by vertical
dashed lines. From these observations we conclude that Turing patterns
observed in the system (2) with values of model parameters used to produce
plots in Fig.6 have the characteristic length $\Lambda_{0}=5.5$ (size of half-
spike) and the amplitude $A=2\alpha_{max}=0.77$ (compare with the amplitude of
all three patterns shown in Fig. 5).

In this section we have performed Fourier decomposition of patterns obtained
in numerical simulations. Our next task is to find the Fourier coefficients
describing Turing pattern analytically. This would let us to compare
analytical results with numerical ones and besides predict the dependence of
the wavelength, $\Lambda_{0}$ , and amplitude, $A$, of Turing patterns on
model parameters in (2).

##  6 Nonlinear analysis of Turing patterns

Stationary solutions of the system (2) represented by their Fourier series,
satisfy the system (2) where partial derivatives over time are set to zero and
variables $n$ and $c$ are replaced by their Fourier series (9). This gives the
following equations for the coefficients of Fourier series, $\alpha_{i}$ and
$\beta_{i}$:

$\begin{cases}\displaystyle-D\sum_{i=0}^{M}(i\kappa)^{2}\alpha_{i}\cos(i\kappa
x)+\chi_{0}\frac{\partial}{\partial
x}\left(\sum_{i=0}^{M}\alpha_{i}\cos(i\kappa
x)\sum_{i=0}^{M}i\kappa\beta_{i}\sin(i\kappa x)\right)+\\\ \hfill
r\sum_{i=0}^{M}\alpha_{i}\cos(i\kappa
x)\left(1-\sum_{i=0}^{M}\alpha_{i}\cos(i\kappa x)\right)=0,\\\\[10.0pt]
\displaystyle-\sum_{i=0}^{M}(i\kappa)^{2}\beta_{i}\cos(i\kappa
x)+\sum_{i=0}^{M}\alpha_{i}\cos(i\kappa x)-\sum_{i=0}^{M}\beta_{i}\cos(i\kappa
x)=0,\end{cases}$ (11)

where $\kappa=\pi/L$. We will start analysis of this system by finding the
homogeneous solutions, or the solutions for which $\alpha_{i}=\beta_{i}=0$
when $i>0$. This means that in order to find the homogeneous solution we
truncate the above system at $M=0$ and get:

$\begin{cases}\alpha_{0}-\alpha_{0}^{2}=0,\\\\[7.0pt]
\alpha_{0}-\beta_{0}=0.\end{cases}$

Roots of this system correspond to two homogeneous solutions:
$\alpha_{0}=\beta_{0}=1$ (meaning that $u(x,t)=c(x,t)=1$, $\forall x,t$) and
$\alpha_{0}=\beta_{0}=0$ (meaning that $u(x,t)=c(x,t)=0$, $\forall x,t$).

To derive equations for $\alpha_{i}$ and $\beta_{i}$ when we truncate the
system (11) at arbitrary $M$ we use trigonometric identities (to replace
product of cosines by cosines of sums). By equating terms containing
$\cos(ikx)$ for each $i$ we break the first equation in the system (11) into
$M+1$ equations:

$\begin{cases}\displaystyle\alpha_{0}-\alpha_{0}^{2}-\sum_{i=1}^{M}\frac{\alpha_{i}^{2}}{2}=0,\\\\[12.0pt]
\displaystyle-D\alpha_{1}k^{2}+\chi_{0}k^{2}\left(\alpha_{0}\beta_{1}+\sum_{i\geq
1}^{M}\frac{i+1}{2}\alpha_{i}\beta_{i+1}-\sum_{i\geq
2}^{M}\frac{i-1}{2}\alpha_{i}\beta_{i-1}\right)+\\\ \displaystyle\hfill
r\left(\alpha_{1}-2\alpha_{0}\alpha_{1}-\sum_{i\geq
1}^{M}\alpha_{i}\alpha_{i+1}\right)=0,\\\\[12.0pt]
\displaystyle-4D\alpha_{2}k^{2}+\chi_{0}k^{2}\left(4\alpha_{0}\beta_{2}+\alpha_{1}\beta_{1}+\sum_{i\geq
1}^{M}(i+2)\alpha_{i}\beta_{i+2}-\sum_{i\geq
3}^{M}(i-2)\alpha_{i}\beta_{i-2}\right)+\\\ \displaystyle\hfill
r\left(-\frac{\alpha_{1}^{2}}{2}+\alpha_{2}-2\alpha_{0}\alpha_{2}-\sum_{i\geq
1}^{M}\alpha_{i}\alpha_{i+2}\right)=0,\\\\[12.0pt]
\displaystyle-9D\alpha_{3}k^{2}+\chi_{0}k^{2}\left(9\alpha_{0}\beta_{3}+3\alpha_{1}\beta_{2}+\frac{3}{2}\alpha_{2}\beta_{1}+\sum_{i\geq
1}^{M}\frac{3(i+3)}{2}\alpha_{i}\beta_{i+3}-\sum_{i\geq
4}^{M}\frac{3(i-3)}{2}\alpha_{i}\beta_{i-3}\right)+\\\ \displaystyle\hfill
r\left(\alpha_{3}-2\alpha_{0}\alpha_{3}-\alpha_{1}\alpha_{2}-\sum_{i\geq
1}^{M}\alpha_{i}\alpha_{i+3}\right)=0,\\\\[12.0pt] ....\\\ \end{cases}$ (12)

Here the first equation is build as an identity for coefficients of $\cos(0)$,
the second equation - for $\cos(kx)$, $i$-th equation - for $\cos(ikx)$. As
the second equation in (11) is linear it breaks into $M+1$ identical
equations:

$\displaystyle-\beta_{i}(ik)^{2}+\alpha_{i}-\beta_{i}=0.$

This gives $\displaystyle\beta_{i}=\frac{\alpha_{i}}{1+(ik)^{2}}$ so that the
coefficients $\beta_{i}$ can be taken out from the system (12), and we get a
system of $M+1$ equations for $M+1$ unknowns.

By truncating the system (12) at $M=1$ and using the formula
$\displaystyle\beta_{1}=\frac{\alpha_{1}}{1+k^{2}}$ we get the following
system for finding coefficients $\alpha_{0}$ and $\alpha_{1}$:

$\begin{cases}\displaystyle\alpha_{0}-\alpha_{0}^{2}-\frac{\alpha_{1}^{2}}{2}=0,\\\\[7.0pt]
\displaystyle-
Dk^{2}\alpha_{1}+k^{2}\chi_{0}\alpha_{0}\frac{\alpha_{1}}{1+k^{2}}+r(-2\alpha_{0}\alpha_{1}+\alpha_{1})=0.\end{cases}$
(13)

This system has four solutions of which two, namely $(0,0)$ and $(1,0)$,
correspond to homogeneous solutions of the system (2) and two others,

$\displaystyle\alpha_{0}=\frac{(Dk^{2}-r)(1+k^{2})}{\chi_{0}k^{2}-2r(1+k^{2})}\mbox{
}\mbox{and}\mbox{
}\displaystyle\alpha_{1}=\pm\sqrt{2\alpha_{0}(1-\alpha_{0})}$ (14)

\- to spatial periodic patterns forming in (2). The roots associated with
homogeneous solutions were found earlier after truncating the system (12) at
$M=0$ while the solutions (14) correspond to the first mode of the spectral
decomposition of Turing pattern. Thus, equations (14) give the dependence of
$\alpha_{0}$ and $\alpha_{1}$ on model parameters and can be compared with the
results obtained numerically and shown in Figure 6.

Truncation of the system (12) at $M=2$ gives three equations for three
unknowns, $\alpha_{0}$, $\alpha_{1}$ and $\alpha_{2}$:

$\begin{cases}\displaystyle\alpha_{0}-\alpha_{0}^{2}-\frac{\alpha_{1}^{2}}{2}-\frac{\alpha_{2}^{2}}{2}=0,\\\\[7.0pt]
\displaystyle-D\alpha_{1}k^{2}+\chi_{0}k^{2}\left(\alpha_{0}\beta_{1}+\alpha_{1}\beta_{2}-\frac{1}{2}\alpha_{2}\beta_{1}\right)+r(\alpha_{1}-2\alpha_{0}\alpha_{1}-\alpha_{1}\alpha_{2})=0,\\\\[7.0pt]
\displaystyle-4D\alpha_{2}k^{2}+\chi_{0}k^{2}(4\alpha_{0}\beta_{2}+\alpha_{1}\beta_{1})+r\left(\alpha_{2}-\frac{\alpha_{1}^{2}}{2}-2\alpha_{0}\alpha
2\right)=0.\end{cases}$

Similarly, systems of equations for 4 and 5 unknowns, when the system (12) is
truncated at $M=3$ or $M=4$, can be derived. Although roots of these system
can’t be expressed explicitly, their values (for a given set of model
parameters) can be found using Matlab code. Using such code we have find
solutions of the system (12) truncated up to $M=4$. These solutions for the
model (2) with values of model parameters $D=1$, $\chi_{0}=1.9$ and $r=0.1$
(see Figs 5 and 6) and domain size $L=5.5$ are shown in Table 2. For
comparision we have also shown in this table the coefficients $\alpha_{i}$
obtained by spectral decomposition of the numerically simulated pattern (shown
in Fig.5B). As we can see from this table, the obtained solutions are
consistent across different truncations and fairly comparable with the results
of spectral decomposition of simulated pattern. Also, the coefficients
$\alpha_{3}$ and $\alpha_{4}$ are considerably small and this indicates that
the truncation of system (12) at $M=4$ should give a fairly accurate solution
for the model (2) when the domain size is reasonably small (in the presented
case $L=\Lambda_{0}$).

| $\alpha_{0}$ | $\alpha_{1}$ | $\alpha_{2}$ | $\alpha_{3}$ | $\alpha_{4}$
---|---|---|---|---|---
Simulation | 0.9009 | 0.3849 | 0.1423 | 0.0428 | 0.0143
$M=1$ | 0.8462 | 0.5103 | 0 | 0 | 0
$M=2$ | 0.8868 | 0.4153 | 0.1680 | 0 | 0
$M=3$ | 0.8938 | 0.4035 | 0.1563 | 0.0508 | 0
$M=4$ | 0.8941 | 0.4033 | 0.1550 | 0.0496 | 0.0155
Table 2: Coefficients of Fourier series, $\alpha_{i}$, for the first five
modes of Turing pattern represented by half-spike. Coefficients in the first
row obtained using eq. (10) for the profile shown in Fig. 5. Values in the
following four rows obtained by truncation of the system (12) at $M=1$, $2$,
$3$ and $4$. Values of the model parameters are the same as in Fig. 5B.

To check the accuracy of truncated solutions presented in Table 2 we have
reproduced corresponding profiles (for $M=1$, $2$ and $3$) and compared them
with the profile obtained in numerical simulations. All these profiles are
shown in Fig. 7A. Visually it is evident that the solution of (12) truncated
at $M=2$ fits better the profiles obtained in numerical simulation than the
one truncated at $M=1$. Furthermore, the profile corresponding to the solution
truncated at $M=3$ is almost indistinguishable from the simulated profile. As
a numerical measure of discrepancy between profiles, $n_{i}(x)$ and
$n_{0}(x)$, one can use the integral:

$I=\int_{0}^{L}[n_{i}(x)-n_{0}(x)]^{2}dx,$

where $n_{i}(x)$ correspond to one of the solutions of (12) (truncated at
$M=i$, $i>0$), while $n_{0}(x)$ \- to the profile obtained in numerical
simulations. The value of this integral for the solutions truncated at $M=1$
is $I=0.1353$, at $M=2$ \- $I=0.0114$ and at $M=3$ \- $I=0.0026$. Thus, the
discrepancy between the solution, truncated at $M=3$, and the one obtained
numerically is lesser than 1%, and besides the numerical discrepancy for $M=4$
is $I=0.0020$ and fairly the same as for $M=3$.

The amplitude of the profile for variable $n$ in a small domain (when
$L\approx\Lambda_{0}$) is predominantly defined by the coefficient
$\alpha_{1}$, while its average value by $\alpha_{0}$. Figure 7B shows the
dependence of these two coefficients on the size of the domain, $L$, obtained
from the spectral decomposition of the simulated profile and according to the
equation (11) truncated at $M=1$ and $M=3$. We see that while all three plots
are qualitatively similar, the plots for simulated profile and for the
solution of eq. (11) truncated at $M=3$ also show a good quantitative
agreement.

Figure 7: Comparison of analytical and numerical solutions. A: Profiles of
solutions to (12) truncated at $M=1$ (red) and $M=3$ (cyan) as compared to
numerically simulated profile (blue, same as in Fig. 5B). B: Dependence of the
coefficients $\alpha_{0}$ (solid lines) and $\alpha_{1}$ (dashed lines) on the
domain size $L$ found by spectral decomposition of profiles obtained in
simulations (blue) and from (12) truncated at $M=1$ (red) and $M=3$ (cyan).
Values of model parameters are the same as in Fig.5.

In this section we have focused on finding the analytical solutions of the
system (2) represented by truncated Fourier series. Using Matlab code we were
able to find solutions truncated up to $M=4$, and have shown that these
solutions (at $M\geq 3$) are fairly accurate if the modelled domain is small.
Throughout this section we have fixed the model parameters: $D=1$,
$\chi_{0}=1.9$ and $r=0.1$. In the next section we will check how the
characteristics of Turing patterns obtained numerically and given by the
analytical solutions of (11) depend on the model parameters.

##  7 Impact of model parameters on characteristics of Turing pattern

Based on the linear analysis of pattern formation presented in Section 3 we
were able to predict dependence of the pattern’s wavelength on model
parameters (see eqs (7), (8) and Fig. 4). However, linear analysis didn’t let
us to make any predictions concerning the amplitude of Turing pattern. In this
section we will make such predictions using nonlinear analysis, presented in
the previous section, and compare them with results of numerical simulations.
The outcome of this study is presented in Fig. 8. Since $D$ represents the
ratio of the rate of random motion of cells to the diffusion coefficient of
chemotactic agent it is essentially smaller than one and therefore $D<1$ in
our plots. Also, ranges of values for $\chi_{0}$ and $r$ were chosen in a way
that the condition for Turing instability given by the equation (6) (and
illustrated in Fig. 2) is satisfied, that is $\chi_{0}>1.9$ or $r<0.14$ for
the default set of parameters.

Figure 8: Dependence of the characteristics of Turing pattern on model
parameters. A-C Dependence of the wavelength on the diffusion coefficient,
$D$, (panel A), chemotactic coefficient, $\chi_{0}$, (panel B), and
reproduction rate, $r$, (panel C): black lines - from linear analysis (same as
solid line in Fig. 4), green lines - from Fourier series truncated at $M=3$,
and blue lines - from simulations. D-F Coefficients of Fourier series versus
the diffusion coefficient, $D$, (panel D), chemotactic coefficient,
$\chi_{0}$, (panel E), and reproduction rate, $r$, (panel F) found after
spectral decomposition of numerically simulated profile (solid lines) and
found analytically from eq. (12) with $M=3$ (dashed lines). G-I Profiles of
cell density obtained numerically for different values of the diffusion
coefficient, $D$, (panel G), chemotactic coefficient, $\chi_{0}$, (panel H),
and reproduction rate, $r$, (panel I). Each profile obtained in the domain of
size of characteristic length, $L=\Lambda_{0}$ (for given set of parameter
values) and presented versus scaled spatial variable $\xi=x/\Lambda_{0}$.

The dependence of pattern’s wavelength on model parameters as predicted on the
basis of linear analysis (eq. (8)), nonlinear analysis (eq. (12)) and found
numerically is shown on panels A-C of Fig. 8. The wavelength in numerical
simulations was found as a length of the domain for which Fourier
decomposition of obtained profile gave maximal value of $\alpha_{2}$ (or
$2\Lambda_{0}$ according to Fig. 6). Similarly, using equation (12) truncated
at $M=4$ we found the value of $k$ for which the value of $\alpha_{2}$ is
maximal, and the wavelength was found as a length of domain corresponding to
this value of $k$, namely, $\Lambda=\pi/k$. As we can see from Fig. 8A-C the
dependences found this way indicate that the wavelength is increasing
(linearly) with the diffusion coefficient, $D$, and decreasing with
chemotactic coefficient, $\chi_{0}$, (can be approximated by the hyperbola)
and reproduction rate, $r$, (appears as a linear) as it was predicted by
linear analysis in Section 3. Concerning the dependence on the reproduction
rate, $r$, the prediction made by the implicit formula (8) rather than by the
explicit (7) (see Fig. 4) was confirmed both by simulations and analytically.
It is evident that the predictions made using nonlinear analysis (on the basis
of equation (12)) are far more accurate as compared to those made using linear
analysis (on the basis of equation (8)) in reproducing numerical results.

The dependence of the amplitude of Turing pattern on model parameters is
addressed on panels D-I of Fig. 8. The amplitude is predominantly defined by
the maximal value of $\alpha_{i}$ ($i>0$) when the size of the domain $L$ is
varied (see Fig. 6). As it was found in Section 5, for the default values of
model parameters, the maximal value of $\alpha_{i}$ is $\alpha_{max}=0.385$
for $i=1$, $2$, $3$, $4$ and evidently for all other modes. So, the main
question we address now is how $\alpha_{max}$ depends on the model parameters.
Panels D-F show how $\alpha_{max}$ which is given here by the value of
$\alpha_{1}$ depends on model parameters. In addition we have presented here
the values of $\alpha_{0}$, $\alpha_{2}$ and $\alpha_{3}$ which contribute to
the shape of the $n-$profile. We do not consider contributions of $\alpha_{2}$
and $\alpha_{3}$ to the amplitude of Turing pattern for the following reasons:
(1) due to symmetry in Fourier expansion, even modes do not contribute to the
amplitude at all and therefore $\alpha_{2}$ has no impact to the amplitude of
pattern; (2) according to Table 2 $\alpha_{3}$ by order of magnitude smaller
than $\alpha_{1}$ which means that it contributes only about 10% to the
amplitude of patterns. Besides, the plots in panels D-F indicate that this
contribution is roughly the same for any set of model parameters and thus we
can consider $\alpha_{1}$ as an indicator for the behaviour of the pattern’s
amplitude.

The main conclusion drawn from the plots presented on panels D-F of Fig. 8 is
that the amplitude of the Turing pattern (given by the value of $\alpha_{1}$)
saturates at low values of the diffusion, $D$, and reproduction rate, $r$, and
at high values of the chemotactic sensitivity, $\chi_{0}$, while it increases
with the increase of $D$ and $r$ and decreases with the increase of
$\chi_{0}$. Another notable conclusion, one can draw from the plots presented
on these panels, is that the average value of cell density, $n$, given by the
value of $\alpha_{0}$, is increasing with the diffusion coefficient $D$ and
reproduction rate, $r$, while decreasing with the chemotactic sensitivity,
$\chi_{0}$. It saturates at low values of $D$ and high values of $\chi_{0}$
while depends linearly on $r$.

Panels G-I of Fig. 8 show how the shape of $n-$profile depends on model
parameters. Plots on panel G confirm that the amplitude of the profile
decreases and the average value of the density of cells, $n$, increases with
the increase of diffusion coefficient, $D$. Plots on panel H confirm that the
amplitude of the profile increases and the average value of the density of
cells, $n$, decreases with the increase of chemotactic sensitivity,
$\chi_{0}$. And finally, plots on panel I confirm that the amplitude of the
profile decreases and the average value of the density of cells, $n$,
increases with the increase of reproduction rate, $r$. Dependence of the
wavelength of patterns on model parameters can’t be seen on panels G-I as the
profiles there, although obtained for media of different sizes, are spatially
scaled and presented as $n(x)=n(\xi)$ in scaled variable $\xi=x/L$.

##  8 Discussion

In this paper, we have analysed periodic stationary patterns forming in a
model of growing population of motile bacteria. The growth of modelled
population is following the logistic law while migration of cells is due to
chemoattraction by the chemical produced by the cells themselves. In terms of
mathematics, we have studied stationary solutions of nonlinear two variable
reaction-diffusion-advection system. It is known that nonlinear system without
advective term should have at least cubic terms to produce periodic patterns
of finite amplitude [30, 26]. It is also known that for formation of periodic
patterns the advection term should represent chemoattraction [18], while
chemorepulsion only causes formation of moving bands of cells [17, 14,
Vasiev2016].

In our analysis of the models of type (2) we have identified conditions for
formation of Turing patterns (or for Turing instability) for a number of
commonly used models [13]. This conditions where presented as a certain
expression, which we denoted as $R_{T}$, given in terms of model parameters
and such that $R_{T}>1$ for Turing pattern to emerge. The expressions defining
$R_{T}$ for a number of models are presented in Table 1. It appears that the
domains in parameter space defined by the condition $R_{T}>1$ for different
models are very similar, which is illustrated in Fig. 2 for models $M3$ and
$M7$. The main conclusion of this analysis is that for Turing patterns to
form, the chemoattraction should be strong enough, while the diffusion of
chemical and the reproduction of cells are slower than certain threshold rate
defined by the condition $R_{T}>1$.

We have also presented our analysis of spatial characteristics of Turing
patterns. Commonly, the wavelength of Turing pattern is estimated on the basis
of linear analysis using formula (7) [23]. We have suggested a correction to
this estimate, which is also based on linear analysis and given by the
implicit formula (8). Comparison with the wavelengths of simulated patterns
has indicated that the formula (8) gives better approximation to the
wavelength of Turing pattern not only quantitatively but also qualitatively
when it comes to checking the dependence of the wavelength on model parameters
(see Fig. 4 and Fig. 8, panel C). The wavelength of Turing pattern was also
estimated on the basis of nonlinear analysis presented in Section 6. The first
order approximation can be made analytically by finding the value of $k$
corresponding to the maximal value of $\alpha_{1}$ in equations (14). The link
between this value of $k$ with the wavelength is given by formula
$k=\pi/\Lambda_{0}$ where $\Lambda_{0}$ is a characteristic length, which is
twice smaller than the wavelength. Dependence of $\alpha_{1}$ on the domain
size, $L$ found on the basis of this approximation is shown by dashed-red plot
in Fig. 7, panel B. It is evident from this figure that the length of the
medium for which $\alpha_{1}$ reaches its maximal value (and which gives the
value of $\Lambda_{0}$) is much higher that the one found numerically
($\Lambda_{0}=5.5$, see Fig. 6). However when we find the characteristic
length using the model (12) truncated at $M=3$ we find that it is very close
to the one obtained numerically (compare locations of maxima for dashed blue
and cyan in Fig. 7, panel B).

One of the most important problems addressed in this work is about analytical
estimation of the amplitude of Turing pattern. This problem gets recently
attention of many researchers. Using multi-scale analysis by means of
amplitude equation for analysis of Turing patterns turned to be technically
sophisticated and factually not convincing [7, 31], although further
development of this method may turn to be more useful. The linear approach for
finding the amplitude of Turing pattern suggested in [6] looks somehow
oversimplified. In terms of the system (2) it brings to the equation for
$\alpha_{1}$ given by the second equation in system (13) where the value of
$\alpha_{0}$ is taken from the homogeneous solution, $\alpha_{0}=1$. We have
found that the value of $\alpha_{1}$ found this way fits numerical results
worth than the one found as a solution of the system (13). Furthermore, the
amplitude found from system (12) truncated at $M=3$ is far better even when
compared with the one found using system (13) and the Turing pattern
reproduced using this solution is almost indistinguishable from the one
obtained in numerical simulations (see Fig. 7, panel A). The only concern
about finding solutions of algebraic system (12) is that this solution can be
expressed explicitly only when this system truncated at $M=1$ (see (14)),
while when $M>1$ solutions can be found only numerically using Matlab (or
Maple) code. This is especially true about the dependence of wavelength and
amplitude of Turing patterns found for the system (12) and presented in Fig.
8. All corresponding plots are drawn on the basis of sets of solutions of
system (12) truncated at $M=3$, which where found numerically.

In this paper we have presented the analysis of spatial characteristics and
amplitude of Turing patterns forming in the model given by system (2). As this
model is a representative for a wide class of models (see Table 1), presented
results, at least, qualitatively are applicable to many other models. The
methods used for the presented analysis, particularly, representation of
Turing patterns by means of Fourier series, appeared to be very helpful. It
can be used for the analysis of other classes of models where formation of
Turing patterns can be observed. Particularly, it can be extended for the
analysis of patterns forming in a system of two competing species of which one
is producing chemical agent and the other responds to it chemotactically. We
are planning to present such study in the follow up publication.

## Acknowledgment

This work has been supported by the EPSRC (Engineering and Physical Sciences
Research Council) and NBIC (National Biofilms Innovation Centre) scholarships
to Valentina Bucur. Authors are grateful to Prof Rachel Bearon for helpful
discussions.

## References

  * [1] M.J. Ablowitz and A. Zeppetella.  Explicit solutions of fishers equation for a special wave speed.  Bulletin of Mathematical Biology, 41(6):835–840, 1979.
  * [2] A. Aotani, M. Mimura, and T. Mollee.  A model aided unedrstanding of spot pattern formation in chemotactic E. coli colonies.  Japan Journal of Industrial and Applied Mathematics, 27:5 – 22, 2010.
  * [3] J.J. Bramburger.  Exact minimum speed of travelling waves in a keller-segel model.  Applied Mathematics Letters, 111:106594, 2021.
  * [4] E.O. Budrene and H.C. Berg.  Complex patterns formed by motile cells of Escherichia coli.  Letters to Nature, 349:630–633, 1991.
  * [5] E.O. Budrene and H.C. Berg.  Dynamics of formation of symmetrical patterns by chemotactic bacteria.  Letters to Nature, 376:49–53, 1995.
  * [6] Y. Chen and J. Buceta.  A non-linear analysis of turing pattern formation.  Plos One, 14(8), AUG 9 2019.
  * [7] A.K. Dutt.  Amplitude equation for a diffusion-reaction system: The reversible Sel’kov model.  AIP Advances, 2(4), 2012.
  * [8] P.C. FIFE.  Asymptotic states for equations of reaction and diffusion.  Bulletin of the American Mathematical Society, 84(5):693–726, 1978\.
  * [9] R.A. Fisher.  The wave of advance of advantageous genes.  Annalf of Eugenics, 7(4):355–369, JUN 1937.
  * [10] B. Franz, C. Xue, K.J. Painter, and R. Erban.  Travelling waves in hybrid chemotaxis models.  Bulletin of Mathematical Biology, 2013.
  * [11] A. Gierer and H. Meinhardt.  Theory of biological pattern formation.  Kybernetik, 12(1):30–39, 1972.
  * [12] A. Gupta and S. Chakraborty.  Linear stability analysis of high- and low-dimensional models for describing mixing-limited pattern formation in homogeneous auticatalytic reactors.  Chemical Engineering Journal, 145:399 – 411, 2009.
  * [13] T. Hillen and K.J. Painter.  A user’s guide to pde models for chemotaxis.  Journal of Mathematical Biology, 58:183 – 217, 2009.
  * [14] N.C. Harrison and R. Diez del Corral and B. Vasiev.  Coordination of Cell Differentiation and Migration in Mathematical Models of Caudal Embryonic Axis Extension.  Plos One, 6(7), 2011.
  * [15] Y. Hosono.  Traveling waves for the lotka-volterra predator-prey system without diffusion of the predator.  Discrete and Continuous Dynamical Systems, 20:161 – 171, 2015.
  * [16] G.A. James, L. Beaudette, and J.W. Costerton.  Interspecies bacterial interactions in biofilms.  Journal of Industrial Microbiology, 15:257 – 262, 1995.
  * [17] E.F Keller.  Traveling bands of chemotactic bacteria: A theoretical analysis.  Journal of Theoretical Biology, 30:235 – 248, 1971.
  * [18] E.F. Keller and L.A. Segel.  Initiation of slime mold aggregation viewed as an instability.  Journal of Theoretical Biology, 26:399 – 415, 1970.
  * [19] E.F. Keller and L.A. Segel.  Model for chemotaxis.  Journal of Theoretical Biology, 30:225 – 234, 1971.
  * [20] A. Kolmogorov, I. Petrovskii, and N. Piskunov.  A study of the diffusion equation with increase in the amount of substance.  Bulletin of Moscow University, 1:1–25, 1993.
  * [21] T. Li and J. Park.  Travelling waves in a chemotaxis model wtih logistic growth.  Discrete and Continuous Dynamical Systems Series B, 24:6465 – 6480, 2019.
  * [22] J.D. Murray.  Mathematical Biology: I. An Introduction.  Springer, third edition edition, 2003.
  * [23] J.D. Murray.  Mathematical Biology: II. Spatial Models and Biomedical Applications.  Springer, third edition edition, 2003.
  * [24] A. Piunovskiy and B. Vasiev.  Modelling ethnogenesis.  Biosystems, 219, 2022.
  * [25] A.A. Polezhaev, R.A. Pashkov, A.I. Lobanov, and I.B. Petrov.  Spatial patterns formed by chemotactic bacteria Escherichia coli.  International Journal of Developmental Biology, 50(2-3, SI):309–314, 2006.
  * [26] M. Rasolonjanahary and B. Vasiev.  Scaling of morphogenetic patterns in reaction-diffusion systems.  Journal of Theoretical Biology, 404, 2016
  * [27] L. Smolin.  Galactic disks as reaction-diffusion systems.  arXiv:astro-ph/9612033.
  * [28] A. Turing.  The chemical basis of morphogenesis.  Royal Society of London, 237:37 – 72, 1952.
  * [29] R. Tyson.  Pattern formation by E. coli; Mathematical and numerical investigation of a biological phenomenon.  PhD thesis, University of Washington, 1996.
  * [30] O. Vasieva and M.I. Rasolonjanahary and B. Vasiev  Mathematical modelling in developmental biology.  Reproduction, 145, 2013.
  * [31] Y. Wang and H. Li.  Travelling wave fronts of lotka-volterra reaction-diffusion system in the weak competition case.  The Royal Society of Edinburgh Section A, 152:912 – 938, 2022.
  * [32] M. Yamaguchi, E. Yoshimoyo, and S. Kondo.  Pattern regulation in the stripe of zebrafish suggests an underlying dynamic and autonomous mechanism.  Applied Physical Sciences, 104:4790 – 4793, 2007.
  * [33] A.N. Zaikin and A.M. Zhabotinsky.  Concentration wave propagation in 2-dimensional liquid-phase self-oscillating system.  Nature, 225(5232):535+, 1970.