Untangling Complex Systems, page 61
t
∂ 2
r
∂ y
r
2 y
1
= β
α
y 1+
xy
xy +
∂
δ
∂
+
2
t
β
∂
r
All the four conditions for having a Turing pattern, considering a diffusion-driven instabil-
ity from the stationary state ( x
) 0 0)
s , ys = ( ,
, are verified. In fact,
I.
α + β = 0 899
.
− 0 9
. 1 < 0;
II.
αβ − γ = ( .
0
)(
899 − .
0 9 )
1 − (− .
0
)
899 > 0;
III.
δ Dβ +αδ = ( .
0
)(
0021 .
0
)
516 (− .
0 9 )
1 + ( .
0
)
899 ( .
0
)
0021 > 0;
IV.
δ Dβ +αδ ≅ × −4 >
δ
D 2 (αβ −γ ) ≅
× −
9 10
2
8 6 10 4
.
The critical wavenumber is K = [
D +
]1 2 ≅
−
(
)((
) (
)) /
1 2
14
1
α δ
β δ
cm . The wavelength is
λ = ( π
2 K ) ≅ 0 449
.
cm. Since the spatial domain is L = [−1, +1], its total length is 2. It
follows that the number of waves in the domain is (2 / λ) ≈ .
4 5. This can be confirmed
by solving numerically our system of partial differential equations, using the following
MATLAB file (adapted from Schneider 2012):
%Grid size
Tf=3000;
a=-1;
% Lower boundary
b=1;
% Upper boundary
M=100;
% M is the number of spaces between points a and b.
dr=(b-a)/M;
% dr is delta r
r=linspace(a,b,M+1); % M+1 equally spaced r vectors including a and b.
%Time stepping
dt=0.04;
% dt is delta t the time step
N=Tf/dt;
% N is the number of time steps
300
Untangling Complex Systems
%Constant Values
D=0.516;
% D is the Diffusion coefficient Du/Dv
delta=0.0021;
% sizes the domain for particular wavelengths
alpha=0.899;
% a is alpha, a coefficient in f and g
beta=-0.91;
% b is beta, another coefficient in f and g
r1=3.5;
% r1 is the cubic term
r2=0;
% r2 is the quadratic term
gamma=-alpha;
% g is for gamma
%pre-allocation
xnp1=zeros(M+3,1);
ynp1=zeros(M+3,1);
%Initial Conditions
xn=-0.5+rand(M+3,1); %Begin with a random point between [-0.5,0.5]
yn=-0.5+rand(M+3,1);
for n=1:N
xn(1)=xn(3);
%Boundary conditions on left flux is zero
xn(M+3)=xn(M+1);
%Boundary conditions on right
yn(1)=yn(3);
yn(M+3)=yn(M+1);
for j=2:M+2
%Source function for x and y
srcx=alpha*xn(j)*(1-r1*yn(j)^2)+yn(j)*(1-r2*xn(j));
srcy=beta*yn(j)*(1+(alpha*r1/beta)*xn(j)*yn(j))+xn(j)*(gamma+r2*yn(j));
Lapx=(xn(j-1)-2*xn(j)+xn(j+1))/dr^2; %Laplacian x
Lapy=(yn(j-1)-2*yn(j)+yn(j+1))/dr^2; %Laplacian y
xnp1(j)=xn(j)+dt*(D*delta*Lapx+srcx);
ynp1(j)=yn(j)+dt*(delta*Lapy+srcy);
end
xn=xnp1;
yn=ynp1;
% Graphing
if mod(n,25000)==0
%subplot(2,1,2)
plot(r, xn(2:M+2),r, yn(2:M+2))
axis([-1,1,-1,1]);
fprintf(‘Time t = %fn’,n*dt);
input(‘Hit enter to continue:’)
end
end
In such MATLAB file, we apply the Euler’s method (see Appendix A) to a semi-discretized
Reaction-Diffusion system. Moreover, the second derivative is expressed as finite differences:
d 2 x
x( r
r
∆ ) 2 x( r) x( r
r
=
∆ )
Lapx =
−
−
+
+
= ( xn( j −1) − 2* xn( j) + xn( j + )
1 ) / dr 2
dr 2
( r 2
∆ )
The initial conditions are fixed as a randomly generated set of values between [−0.5,
+0.5].
The result of the calculation is shown in Figure 9.33.
We notice that the x and y species are spatially distributed in anti-phase conditions:
in the point where x is at its maximum, y is at its minimum and vice versa. This result
verifies that in the places where the activator prevails, the inhibitor succumbs and
vice versa. Moreover, we confirm that the number of waves in the spatial domain is
about 4.5.
The Emergence of Order in Space
301
0.2
xy
sie 0.1
ec
0.0
−0.1
Abundance of sp
−0.2−1.0
−0.5
0.0
0.5
1.0
r
FIGURE 9.33 Turing pattern involving x and y in one-dimensional space of coordinate r.
9.7. If you have MATLAB at your disposal, an example of .m file that can be used to find the
solution of the exercise is the following one (adapted from Schneider 2012):
%Grid size
Tf=100000;
a=-1;
% Lower boundary
b=1;
% Upper boundary
M=100;
% M is the number of spaces between points a
and b.
dx=0.04;
%(b-a)/M; % dx is delta x
dy=0.04;
%(b-a)/M;
x=linspace(a, b,M+1);
% M+1 equally spaced x vectors including a and b.
y=linspace(a, b,M+1);
%Time stepping
dt=0.08;
%100*(dx^2)/2; % dt is delta t the time step
N=Tf/dt;
% N is the number of time steps in the interval
[0,1]
%Constant Values
D=0.516;
% D is the Diffusion coefficient Du/Dv
delta=0.0021;
% sizes the domain for particular wavelengths
alpha=0.899;
% a is alpha, a coefficient in f and g (-a is
gamma)
beta=-0.91;
% b is beta, another coefficient in f and g
r1=3.5;
% r1 is the cubic term
r2=0;
% r2 is the quadratic term
gamma=-alpha;
% g is for gamma
%pre-allocation
unp1=zeros(M+3,M+3);
vnp1=zeros(M+3,M+3);
%Initial Conditions
un=-0.5+rand(M+3,M+3); %Begin with a random point between [-0.5,0.5]
vn=-0.5+rand(M+3,M+3);
for n=1:N
for i=2:M+2
un(i,1)=un(i,3);
%Boundary conditions on left flux is zero
un(i,M+3)=un(i,M+1);
%Boundary conditions on right
vn(i,1)=vn(i,3);
vn(i,M+3)=vn(i,M+1);
end
302
Untangling Complex Systems
for j=2:M+2
un(1,j)=un(3,j);
%Boundary conditions on left
un(M+3,j)=un(M+1,j);
%Boundary conditions on right
vn(1,j)=vn(3,j);
vn(M+3,j)=vn(M+1,j);
end
for i=2:M+2
for j=2:M+2
%Source function for u and v
srcu=alpha*un(i, j)*(1-r1*vn(i, j)^2)+vn(i, j)*(1-r2*un(i, j));
srcv=beta*vn(i, j)*(1+(alpha*r1/beta)*un(i, j)*vn(i, j))
+un(i, j)*(gamma+r2*vn(i, j));
uxx=(un(i-1,j)-2*un(i, j)+un(i+1,j))/dx^2; %Laplacian u
vxx=(vn(i-1,j)-2*vn(i, j)+vn(i+1,j))/dx^2; %Laplacian v
uyy=(un(i, j-1)-2*un(i, j)+un(i, j+1))/dy^2; %Laplacian u
vyy=(vn(i, j-1)-2*vn(i, j)+vn(i, j+1))/dy^2; %Laplacian v
Lapu=uxx+uyy;
Lapv=vxx+vyy;
unp1(i, j)=un(i, j)+dt*(D*delta*Lapu+srcu);
vnp1(i, j)=vn(i, j)+dt*(delta*Lapv+srcv);
end
end
un=unp1;
vn=vnp1;
% Graphing
if mod(n,6250)==0
%subplot(2,1,2)
hdl = surf(x,y,un(2:M+2,2:M+2));
set(hdl,’edgecolor’,’none’);
axis([-1, 1,-1,1]);
%caxis([-10,15]);
view(2);
colorbar;
fprintf(‘Time t = %fn’,n*dt);
ch = input(‘Hit enter to continue:’,’s’);
if (strcmp(ch,’k’) == 1)
keyboard;
end
end
end
The patterns achieved when r
3.5,
0, and there are not quadratic terms in the PDEs,
1 =
r 2 =
are shown in Figure 9.34. They are made of stripes.
The patterns obtained when r
0.02,
0.2, and the contribution of the cubic terms
1 =
r 2 =
in the PDEs is pretty small, are shown in Figure 9.35. They are made of spots.
The patterns obtained when r
3.5,
0.2 and both the cubic and the quadratic terms
1 =
r 2 =
are relevant, are shown in Figure 9.36. They contain both stripes and spots.
In synthesis, when there is only the contribution of the cubic terms, we observe stripes; when
the quadratic term is dominant, we have spots. Finally, when both the quadratic and the cubic
terms contribute appreciably, we observe patterns with spots and stripes (Barrio et al. 1999).
The Emergence of Order in Space
303
u
v
−
−
0.1750
0.2250
−
0.8
0.1344
−0.1666
0.8
−0.09375
−0.1063
−0.05313
−0.04688
−0.01250
0.01250
0.4
0.02813
0.07187
0.4
0.06875
0.1312
0.1094
0.1905
y 0.0
0.2500
0.1500
y 0.0
−0.4
−0.4
−0.8
−0.8
−0.8
−0.4
0.0
0.4
0.8
−0.8
−0.4
0.0
0.4
0.8
x
x
FIGURE 9.34 Profiles of u (on the left) and v (on the right) obtained when r
3.5,
0, and within the
1 =
r 2 =
square box having x and y as spatial coordinates.
u
v
−5.960
−11.23
−
−
3.344
9.547
0.8
−0.7375
0.8
−7.899
1.869
−6.191
4.475
−4.513
0.4
7.081
−2.834
0.4
9.688
−1.156
12.29
0.5219
y
14.90
2.200
0.0
y 0.0
−0.4
−0.4
−0.8
−0.8
−0.8
−0.4
0.0
0.4
0.8
−0.8
−0.4
0.0
0.4
0.8
x
x
FIGURE 9.35 Profiles of u (on the left) and v (on the right) obtained when r
0.02,
0.2, and within the
1 =
r 2 =
square box having x and y as spatial coordinates.
u
v
−0.1740
−0.1660
−
−
0.1218
0.1308
0.8
−0.09950
0.8
−0.09550
−0.01725
−0.06025
0.03500
−0.02500
0.4
0.08725
0.4
0.01025
0.1395
0.04550
0.1918
0.08075
y
0.2440
0.1160
0.0
y 0.0
−0.4
−0.4
−0.8
−0.8
−0.8
−0.4
0.0
0.4
0.8
−0.8
−0.4
0.0
0.4
0.8
x
x
FIGURE 9.36 Profiles of u (on the left) and v (on the right) obtained when r = 3.5, = 0.2, and within the 1
r 2
square box having x and y as spatial coordinates.
304
Untangling Complex Systems
TABLE 9.4
Six Combinations of the Parameters of the Turing Model and Pictures of the Patterns They
Generate. In the Pictures, the Activator is Gray, Whereas the Inhibitor is Black
Name
au
bu
cu
du
av
bv
cv
dv
Du
Dv
Pattern
Spot
0.08
−0.08
0.005
0.03
0.1
0
−0.15
0.06
0.02
0.5
Net
0.08
−0.08
0.2
0.03
0.1
0
−0.15
0.06
0.02
0.5
Stripe2
0.08
−0.08
0.04
0.03
0.1
0
−0.15
0.08
0.02
0.5
Stripe3
0.08
−0.08
0
0.03
0.1
0
−0.15
0.08
0.02
0.5
Dapple
0.15
−0.08
0.025
0.03
0.15
0
−0.15
0.06
0.02
0.5
Spot2
0.1
−0.08
0.025
0.03
0.15
0
−0.1
0.06
0.02
0.5
9.8. The combinations of the parameters available in the RD simulator by Kondo and his team
are listed in Table 9.4. In the last column of Table 9.4, there are the pictures of the patterns generated by numerical solution of the partial differential equations.
A new combination of parameters values is the following one:
au = 0 0
. 8, bu = −0 0
. 8, cu = 0 2
. , du = 0 0
. 3, Du = 0 0
. 2, av = 0 1
. , bv = 0, cv = −0 1
. 5,
dv = 0 0
. 7, Dv = 0 5
. . This combination differs from that named as “Net” in Table 9.4 only in
dv. The new combination verifies the four conditions required to have a Turing pattern. In fact,
The Emergence of Order in Space
305
I.
tr( J ) = au − du + bv − dv = .
0 08 − .
0 03 + 0 − .
0 07 < 0
II.
det( J ) = ( a
)(
)
u − du
bv − dv − buav = ( .
0 0 )(
5 −0. )
07 − (−0.0 )(
8
.
0 )
1 = .
0 00445 > 0
III.
D (
)
(
)
u bv − dv + Dv au − du = (0.0 )(
2 −0. )
07 + (0. )(
5
.
0 0 )
5 = .
0 0236 > 0
IV.
D (
)
(
) 0 0236
.
2
det( )]1/2
/
1 2
[
]
u bv − dv + Dv au − du =
> [ DuDv
J
=2 (0 0
. 2)(0 5
. )(0.
)
0045
= .
0 013
If we start from a randomized condition, the final stationary pattern is shown in
Figure 9.37 as graph a. If we start from configuration a and we reduce the values of the
diffusion coefficient to Du = 0 002
.
and Dv = 0 0
. 5, we obtain the pattern shown in b of
Figure 9.37. On the other hand, if we use the values Du = 0 002
.
and Dv = 0 0
. 5 from an
initial randomized distribution of u and v into space, we obtain the pattern labeled as c in
Figure 9.37.
If we take the pattern obtained by the combination named as “Stripe2” in Table 9.4 as
