Untangling Complex Systems, page 62
the initial condition and we reduce the values of the diffusion coefficients to Du = 0 002
.
and Dv = 0 0
. 5, we achieve pattern D of Figure 9.38. If we do the same from the output of
the combination named as “Stripe3,” we achieve pattern E of Figure 9.38.
(a)
(b)
(c)
FIGURE 9.37 Three Turing patterns obtained for different values of the parameters and distinct initial
conditions. Read the text for details.
Stripe2
D
Du = 0.002
Dv = 0.05
Stripe3
E
Du = 0.002
Dv = 0.05
FIGURE 9.38 Effects of the diffusion coefficients values on the Turing patterns.
306
Untangling Complex Systems
These modifications demonstrate the plasticity of the Turing patterns. Such plasticity has
been experimentally demonstrated by Watanabe and Kondo (2012) in transgenic zebrafish.
9.9. The steady state solution for the homogeneous case can be determined by fixing the two
differential equations, devoid of the diffusion terms, equal to zero:
dX *
* *
*
X Y
= a − X − 4
2 = f ( X *, Y * ) = 0
dτ
1+
X *
dY *
* *
*
bX Y
= bX −
= g ( X * Y *
,
) = 0
*2
dτ
1+
X
The steady state solution is when X *
*
2
s = ( a
)
5 , Ys = (1+ ( a 2 ))
5
This steady state is stable if tr( J ) < 0, det( J ) > 0. The Jacobian is
2
a
1−
25
−
20 a
1−
4
−
f
∂
f
∂
2
a
(25 + 2
a )
2
3 a − 5(25)
20
1+
−
*
*
X
∂
Y
∂
25
2
(25 + a )
(225
2
+ a )
J
s
s
=
=
=
g
∂
g
∂
a 2
2 2
a b
5
ba
−
*
*
1−
X
∂
Y
∂
2
2
( + a )
(
+
s
s
25
25 a )
25
ba
5
b −
b
−
a 2
(25 + a 2)
1+
25
The trace of the Jacobian J is negative if b > (3 )
5 a − (25 a); the determinant of J is positive
if a > 5.
9.10. The values of the kinetic constants and the concentrations presented in the text of
the exercise allow us calculating the rate constants that appear in [9.32]. We find that
′ = 7 1
. 4 × −6 −
k
1
1
6
1
1
10 s M, ′ =
−
k 2 3 6
. s , ′ = 2 6
. 5× − −
k 3
10 s M. It derives that the parameters
a = ( k′
( ′3/ ′2 α )
1 / k′2
α ) and b = k k
are equal to 19.8 and 7.36, respectively. Note that a and
b are dimensionless. Now, we determine the values of the terms of the Jacobian J.
2.68
−0.96
J =
14.14
−
1.75
Note that two of the four terms of the Jacobian are positive: they regard the activator I −.
The other two terms are negative, and they regard the inhibitor ClO−2. det( J ) = 8.88 > 0.
The determinant J’ becomes negative if the relationships [9.37] and [9.38] are both
verified. The first relation imposes that d > 0.65. The second one imposes that d > 6.
9.11. The PDEs are:
∂[ X ]
2
2
2
∂ [ X ]
∂ [ X ]
k 1 k 2 X
k 3 X
Y
DX
2
2
D
∂ = ′ − [ ] + [ ] [ ] +
X
t
∂
+
r 1
∂
r 2
∂[ Y ]
2
2 Y
2
∂ [ Y ]
∂ [ ]
= k′4 − k 3[ X ] [ Y ] + DY
DY
2
2
+
∂
t
∂
r 1
∂
r 2
The Emergence of Order in Space
307
where
k′1 = k 1[ A] and k′4 = k 4[ B]. Note that in the original PDEs we have five variables and six parameters. To obtain the nondimensionalization of the PDEs, we define the following
dimensionless variables: u = X [ ]
[ ]
C X , v = YC Y , τ = tCt, x = xC r 1, y = yC r 2.
Introducing these new variables in the PDEs, we get:
∂ u
X
2
D 2
2
2
2
X x
C
u DX y
C
u
C
u
u v
∂
=
k′1
k 2
k 3
2 +
2
∂τ
−
+
+
∂
tC
tC
XCtC YC
tC
∂
x
tC
∂
y
∂ v
Y
u 2 v D x 2 2 v D y 2 ∂2 v
= k′
C
Y C
Y C
4
k
−
3
+
∂
2 +
2
∂τ
tC
XC tC
tC x
∂
tC y
∂
The next step requires the definition of nondimensionalizing constants. Since L is the
length of the sides of the square, it derives that xC = yC = (1 L). Introducing these defini-
tions of x
2
C and yC in the PDEs, it is spontaneous to fix tC = ( DX / L ). If we also impose d = ( D / )
Y DX , we rewrite the PDEs as
∂ u k
2
2 u 2 u
2 k′1
k
3
u v
=
XC − u +
+ ∂
2 + ∂
2
∂τ
tC k 2
k 2 XC YC ∂
x
∂
y
∂ v k
2
2 v
2 v
2 k′
k
u
=
4
YC − 3
v +
∂
d
d
2 +
∂
2
∂τ
t
C
k
k X
∂ x
∂ y
2
2
C
If we fix X
1 2
1 2
3/
/
2 )
3/
/
2 )
C = ( k k
and YC = ( k k
, we have:
∂ u
k
1/2
2
2 u
2 k′1 k 3
2
u
=
u u v + ∂
+ ∂
2
∂τ
− +
t
2
C
k 2 k 2
∂
x
∂
y
∂ v
1/2
k
2 v
2 v
2 k′4 k 3
=
2
u v +
∂
d 2 +
∂
d 2
∂
−
τ
tC k k
∂ x
∂ y
2 2
Finally, if we fix γ = ( k
1/2
( ′
1/2
4
2 )( 3
2 )
2 / tC ), a = (( k′
, and b = k k
k k
, the PDEs
1 k 2 )( k 3 k 2 ))
assume their final nondimensionalized form:
∂ u
2
2
2
u ∂ u
= γ { a − u + u }
v + ∂
2 + 2
∂τ
∂
x
∂
y
∂ v
2
2
2
v
v
= γ b
{ − u v}+ ∂
d 2 +
∂
d 2
∂τ
x
∂
y
∂
wherein we still have five variables, but now only four parameters.
9.12. The nondimensionalized ordinary differential equations of the Schnackenberg model are:
∂ u
= γ { a − u + u 2 }
v
∂τ
∂ v
= γ { b − u 2 }
v
∂τ
308
Untangling Complex Systems
If we put the ODEs equal to zero, we find that the steady-state solution is
b
us = a + b, vs =
+ > 0.
( a + b)2 with b > 0 and a b
The Jacobian for the Schnackenberg model is:
( b − a)
( a + b 2
)
u
2
1
2 ( +
svs −
us
a b)
J =
=
−
u
2
2
−
svs
− us
b
2
−( a + b)2
( a + b)
Since
[ R
]
[
]
u U
( ) and Rv V
( ) must have opposite signs, we must have b > a. The conditions
s
s
required to have diffusion-driven instability are:
( b a)
I.
tr( J ) =
−
− ( a + b)2 < 0, i.e., ( b − a) < ( a + b)3;
( a + b)
II.
det( J ) = ( a + b)2 > 0;
III.
[ R
]
[
]
3
v V
( ) + d R U
( ) > 0, i.e., d( b − a) > ( a + b) ;
s
u
s
2
2
IV.
{ R
3
2
( )
} 4 ) ( )
(
) (
)
4 (
)
v V
d R U
> ( d det J
+ ( )
, i.e., d b − a > a + b
d a b
s
u
s
>
+
The relations (I)–(IV) define a domain in ( a, b, d) parameter space that is favorable to
observe Turing structures.
9.13. I report an example of .m file that is useful to solve the exercise by MATLAB as follows.
function []=RDzeroflux(N,a,b,d,gamma)
% defining the mesh
N = 100;
h = 1/N; % step size in x and y
x = h*(0:N); % x coordinates of grid
y = h*(0:N); % y coordinates of grid
[xx, yy] = meshgrid(x, y); % 2D x and y mesh coordinates
dt =.01*h^2; % time step - usually small
% parameter values
a =.1;
b = 1.35;
d = 20;
gamma = 10000;
% Initial data at t=0:
us = (a+b); % u steady state
vs = b/us^2; % v steady state
u = (a+b)*ones(size(xx)); % u steady state
v = (b./u.^2); % v steady state
u = u +.1*randn(size(xx)); % add small perturbations about the steady
state
v = v +.1*randn(size(xx)); % add small perturbations about the steady
state
% Time-stepping:
t = 0;
tmax = 10;
nsteps = round(tmax/dt); % number of time steps
for n = 1:nsteps % main time-stepping loop
t = t+dt;
uE = u(:,[2:N+1 N]);
The Emergence of Order in Space
309
uW = u(:,[2 1:N]);
uN = u([2 1:N],:);
uS = u([2:N+1 N],:);
vE = v(:,[2:N+1 N]);
vW = v(:,[2 1:N]);
vN = v([2 1:N],:);
vS = v([2:N+1 N],:);
% finite difference formula
u2v = u.^2.*v;
u = u + gamma*dt*(a-u+u2v)+dt*(uE+uW+uN+uS-4*u)/h^2;
v = v + gamma*dt*(b-u2v)+d*dt*(vE+vW+vN+vS-4*v)/h^2;
end;
% plot solution:
colormap(‘grey’)
surf(x, y,u)
title([‘u at t = ‘num2str(t)],’fontsize’,16)
zlim([0 4])
shg
end
