Untangling Complex Systems, page 49
9.2 THE REACTION-DIFFUSION MODEL
The time evolutions of the many systems we studied in the previous chapters give rise to diversified
dynamical landscapes. Such landscapes may be imagined as mountain ranges. Locally, there are
many shapes, like those sketched in Figure 9.2. Each one represents a particular dynamical situ-
ation. We may stumble across a stable periodic orbit in a groove (like in (a)); an unstable periodic
orbit on a brim of a crater (like in (b)); an unstable fixed point on top of a bowl (like in (c)); a stable
fixed point at the bottom of a bowl (like in (d)); a saddle point (like in (e)). Such shapes may be
mapped by the linear stability analysis applied locally to each fixed point of the dynamical land-
scape. The eigenvalues and the eigenvectors that we determine tell us how the system responds to
local perturbations.
241
242
Untangling Complex Systems
(a)
(c)
(e)
(b)
(d)
(f)
FIGURE 9.1 Examples of chemical waves in (a) and (b), Turing patterns in (c) and (d), and periodic precipita-
tion in (d) and (e) observed in a laboratory (top row) and nature (bottom row).
(a)
(b)
(c)
(d)
(e)
FIGURE 9.2 Possible shapes in a dynamical landscape: a groove in (a), a crater in (b), a bowl in (c) and (d), and a saddle in (e).
If we focus on non-linear chemical reactions enriched by diffusion, we expect to see more com-
plicated, odd and bizarre dynamical landscapes. To unveil them, we use our traditional approach:
we linearize the non-linear system around the steady-state solutions, and we determine the eigen-
vectors and the eigenvalues. The eigenvalues tell us how fast the system rolls along the eigenvectors.
Let us imagine having non-linear chemical reactions taking place inside an open and unstirred
system of volume V and area A (see Figure 9.3). Let us indicate with m the concentration of M spe-
cies within V. The change in the amount of m depends on the flow of M in ( J
)
(
)
o i
→ and out Ji→ o
of V, and on the rates of all the chemical reactions involving M within V. If we indicate with J
( →
→ )
M = Ji o − J o i , the net flow through the infinitesimal area dA, and with R( M) the net rate of chemical transformations involving M within the infinitesimal volume element dV (see Figure 9.3), the balance equation, based on the Conservation Law of Mass, is
The Emergence of Order in Space
243
nˆ
J→ M
dA
M
dV
R( M)
V
FIGURE 9.3 Terms of the balance equation for M: J M represents the net flow of M out of V in a small area dA having n as its normal vector; R( M) represents the net rate of M’s production within the volume V.
m
∂
dV = R( M ) dV − JM dA
∫
[9.1]
∂
∫
∫
t
V
V
A
According to the divergence theorem (see Box 9.1 of this chapter)
J
( )
M dA =
∇ JM dV
∫
∫
[9.2]
A
V
BOX 9.1 THE DIVERGENCE THEOREM
For the demonstration of the divergence theorem, let us consider the flow of M out of a tiny
volume having the shape of a parallelepiped with dx, dy, and dz as sides lengths (see the following figure, within this box). Let us start to calculate the flow of M through the two sides that are
orthogonal to the x-axis. The flow through the ABCD side (having A
as area) is given by:
ABCD
dΦ
, ( , ,
)
ABCD = J M dAABCD = − J M x x y z dydz
[B1.1]
where x represents the position of the ABCD face with respect to the x-axis; y and z are the average values of the coordinates y and z for the ABCD face. Since the normal vector to the
ABCD face has opposite orientation with respect to the x-axis, we need to introduce the minus
sign in equation [B1.1].
y
C
G
F
B
D
H
A
E
x
x + dx
x
z
( Continued)
244
Untangling Complex Systems
BOX 9.1 (Continued) THE DIVERGENCE THEOREM
Similarly, the flow of M through the EFGH side, having x + dx as its own coordinate along the x-axis, is
dΦ
, (
, , )
EFGH = J M dAEFGH = + J M x x + dx y z dydz
[B1.2]
The equation [B1.2] can be transformed into the equation [B1.3] by developing
JM, x( x + dx, y, z ) in Taylor series with ( x, y, z ) as its accumulation point, and truncating the development at its first term, i.e., at its first derivative:
J
∂ ,
d
M x
Φ
, ( , , )
EFGH = J M x x y z dydz +
dxdydz
[B1.3]
∂
x ( x, y, z)
If we sum the contributions dΦ ABCD and dΦ EFGH, we obtain
J
∂ ,
d
M x
Φ ABCD + dΦ EFGH =
dV
[B1.4]
∂
x ( x, y, z)
Analogously, we may determine the contributions of J M through the pairs of faces that are
orthogonal to the y- and z- axis, respectively. Finally, the total flux will be
J ,
,
,
M x
J
∂ M y
J
d
M z
Φ tot = ∂
+
+ ∂
dV
[B1.5]
x
∂
y
∂
z
∂
The term ∂ JM, x
∂ JM , y
∂
( + + JM, z is the divergence of J
( ). The flux of M can also be
x
y
z )
∂
∂
∂
M : div J M
defined by equation [B1.6]:
dΦ tot = JMdAtot
[B1.6]
The term dAtot represents the sum of the area of the six faces of the parallelepiped. By
integration, we conclude that the flow through the macroscopic area A, encompassing the
volume V, is
J
( )
M dA =
div JM dV
∫
∫
[B1.7]
A
V
Note that the divergence is a differential operator. When it is applied to the vector J M, it gives
the scalar function: div
( J )
M . We can use the nabla ∇ operator, which is a vector differential
operator
∂
∇ = i + ∂
j + ∂
k , to formulate the divergence theorem. With ∇, the divergence
∂ x
∂ y
∂ z
theorem is reported in equation [9.2]. The divergence of a vector, like J M, can be expressed as
a scalar product between the nabla operator ∇ and the vector itself.
The Emergence of Order in Space
245
Introducing equation [9.2] into [9.1], we may rewrite the balance equation in the following form:
m
∂
dV = R( M ) dV − (∇ J )
M dV
∫
[9.3]
∂
∫
∫
t
V
V
V
Since equation [9.3] must be valid for any volume, we can equate the integrands. This step gives us
the differential form of the balance equation:
∂ m
= R( M ) − ∇ JM
[9.4]
∂ t
Reminding Fick’s law [3.66] ( J M = − DM∇ m), and introducing it in [9.4], we obtain
∂ m
= R( M ) + D
2
M ∇ m
[9.5]
∂ t
The operator ∇2 = ( ∂2
2
2
is the Laplacian, which is the sum of the second derivatives with
2 + ∂ 2 + ∂
∂
∂
∂ 2
x
y
z )
respect to the spatial variables.
Now, we need to express the term R( M). We imagine having nonlinear chemical processes inside
the volume V. For instance, we may have an autocatalytic reaction of the type B + M → 2M, whose
rate is kbm (where k is the kinetic constant, and b is the concentration of B). For simplicity, we imagine having just one spatial coordinate, r. Hence, the balance equation [9.5] becomes
∂ m
2 m
= kbm +
∂
DM
[9.6]
2
∂ t
∂
r
To find the steady-state solution, we impose (∂ m ∂ t) = 0. It derives that
2 m
kbm = −
∂
DM
[9.7]
2
r
∂
The solution of equation [9.7] is a function whose second derivative with respect to r is equal to the
same function with the changed sign. A solution is
2π
m = m 0 sin
r +ϕ
[9.8]
λ
In fact, if we introduce [9.8] in [9.7], we obtain
∂2 m
2π 2
2π
kbm
2
0
m
sin
r ϕ
π
sin
r ϕ
[9.9]
2 = −
0
+
= −
+
∂
r
λ
λ
DM
λ
The meaning is that a stationary sinusoidal pattern emerges. Its wavelength is
λ = π D
2
M
[9.10]
kb
246
Untangling Complex Systems
If the mono-dimensional system extends in the range [− L,+ L] and the boundary conditions are
∂ m
∂
(
m
0. Then, it derives that
r )
= ( r ) =
∂
+
∂
L
− L
∂ m
∂ m
2π
2π
=
=
m 0
cos
(± L) +
ϕ = 0
[9.11]
∂ r
λ
λ
+ L
∂ r − L
The latter equation is verified when
2π ±
n
L +ϕ
π
(
)
with n 1,3,5,
= ±
=
…
[9.12]
λ
2
4 L
If
ϕ = 0, λ =
[9.13]
n
TRY EXERCISES 9.1 AND 9.2
9.3 TURING PATTERNS
The Reaction-Diffusion (RD) model we have developed in the previous paragraph is rather trivial.
We need more sophisticated models to interpret the emergence of spatially ordered structures, like
Turing patterns, chemical waves, and periodic precipitation.
For the interpretation of the Turing patterns (see, e.g. , Figure 9.1 c and d), we need a model that includes two species, X and Y, participating in non-linear kinetic laws and diffusing freely in space.
Let us indicate the concentrations of the two chemical variables as x and y, respectively. The Partial Differential Equations (PDE) of the type [9.5] become
∂ x
= R( X ) + D
2
X ∇ x
∂ t
[9.14]
∂ y
= R Y
( ) + D 2
Y ∇ y
∂ t
In [9.14], R( X ), DX , R Y
( ), and DY represent the net production rate and the diffusion coefficients
for X and Y, respectively. We simplify the mathematical treatment of the chemical system if we
consider only one spatial dimension that we label as r. Therefore, equation [9.14] transforms into
equation [9.15]:
∂ x
2 x
= R( X ) +
∂
DX 2
∂ t
∂
r
[9.15]
∂ y
2 y
= R Y
( ) +
∂
DY 2
∂ t
∂
r
Let us assume that, initially, the system is in a stable steady-state. The concentrations of X and Y
are ( x
)
s , ys , and the chemical compounds are homogeneously distributed in space. Since we are in a
steady-state, we have that [ (
R X )] s = [ (
R Y )] s = 0: the sums of the rates of the transformations involv-
ing X and Y, both determined at the steady-state ( x
)
s , ys , are equal to zero. Moreover, the steady-
state is stable. Therefore, tr( J ) < 0 and det( J ) > 0,
The Emergence of Order in Space
247
R
( )
x X
R X
( )
where J
s
y
s
=
R
( )
x Y
R Y
( )
s
y
s
is the Jacobian matrix, whose elements are the partial derivatives of the rates evaluated at the
steady-state ( x
)
s , ys .
Now, let us imagine slightly perturbing the system at one point. The terms δx and δy represent the
extents of the perturbations on x and y, respectively. In the absence of diffusion, when the chemical system is under stirring, the perturbations annihilate because the steady-state is stable. What happens in the presence of diffusion?
We expect that diffusion should act to lessen, and eventually annihilate, the inhomogeneities
induced by the perturbations, promoting the restoration of the original uniform distribution of the
chemical species into space. Does the system behave as we expect? The answer is: “it depends.”
Let us have a look at the possible solutions of the partial differential equations linearized about
the steady-state ( x
)
s , ys :
∂2
D
0
∂ δ x R
X
2
x ( X )
δ x
δ x
s
Ry ( X )
s
∂
r
∂
=
+
t
[9.16]
δ y
δ y
2 δ y
R
∂
x ( Y )
