Help plotting interacting Turing modes

[jssdenton72]

Hello, everybody.

I'm trying to replicate plots from

Yang et al. 2002. Spatial resonances and superposition patterns in a reaction-diffusion model with interacting Turing modes. Phys. Rev. 88(20)

using Mathematica, to explore specifically under what parameter values the pattern of Figure 3c is stable.

I am a biologist by training, and the notation of PDEs is a bit beyond me (I am also not sure how to define the initial conditions [perturbations of steady state]), but what I have so far to replicate the plots is (in Mathematica 7):

(* Fixed parameter values used to generate Figure 3c *)

a = 3;
Dx1 = 1.31;
Dy1 = 9.87;
Dx2 = 34;
Dy2 = 344.8999;
b = 6;
\[Alpha] = 1;

(* from equation 1, equations Structured as: Laplacian + interaction + Brusselator *)

e1 = D[c1[x1, y1, x2, y2, t], t] ==
Dx1*(D[c1[x1, y1, t], x1, x1] +
D[c2[x1, y1], y1, y1]) + \[Alpha]*(c3[x2, y2, t] -
c1[x1, y1, t]) +
a - (b + 1)*c1[x1, y1, t] + c1[x1, y1, t]^2*c2[x1, y1, t];

e2 = D[c2[x1, y1, x2, y2, t]] ==
Dy1*(D[c1[x1, y1, t], x1, x1] + D[c2[x1, y1, t], y1, y1]) +
\[Alpha]*(c4[x2, y2, t] - c2[x1, y1, t]) +
b*c1[x1, y1, t] - c1[x1, y1, t]^2*c2[x1, y1, t];

e3 = D[c3[x1, y1, x2, y2, t], t] ==
Dx2*(D[c3[x2, y2, t], x2, x2] +
D[c4[x2, y2, t], y2, y2]) + \[Alpha]*(c1[x1, y1, t] -
c3[x2, y2, t]) +
a - (b + 1)*c3[x2, y2, t] + c3[x2, y2, t]^2*c4[x2, y2, t];

e4 = D[c4[x1, y1, x2, y2, t], t] ==
Dy2*(D[c3[x2, y2, t], x2, x2] +
D[c2[x1, y1, t], y1, y1]) + \[Alpha]*(c2[x1, y1, t] -
c4[x2, y2, t]) + b*c3[x2, y2, t] -
c3[x2, y2, t]^2*c4[x2, y2, t];

InteractingBrusselator = NDSolve[{e1, e2, e3, e4},

{c1, c2, c3, c4}, {x1, x2,y1,y2,t}]

I get several errors, including one that says the equation c1[x1,y1,t] is not a function of all the variables. However, c1 and c2 define the interactions in one layer, and c3 and c4 define interactions in another layer. It is not clear to me how to proceed.

I don't need time dynamics for the equations, just the stable pattern, which should look something like the attached images from the software Ready.

Any help on setting up the system for plotting would be greatly appreciated.

John

 

[Bill Simpson]

Usually a function you define in Mathematica has a fixed number of arguments. But you have examples of c1, c2, c3 and c4 with different numbers of arguments.

{c1[x1, y1, t], c1[x1, y1, x2, y2, t],
 c2[x1, y1], c2[x1, y1, t], c2[x1, y1, x2, y2, t],
 c3[x2, y2, t], c3[x1, y1, x2, y2, t],
 c4[x2, y2, t], c4[x1, y1, x2, y2, t]}

Is there any chance you might be mixing up the numbers of variables that you want to differentiate with the arguments to your function? I'm thinking that based on the warning/error you are getting.

If I'm confused and your code is exactly what you meant then I apologize for my misunderstanding.

 

[jssdenton72]

Hi Bill,

I think my issue comes from not understanding the notation of PDEs in Mathematica, in general. The code was as close as I could get, but there is something wrong. As you point out, the c1, c2, etc. do not match up, but how do I define the elements properly for Mathematica to solve it as a PDE?

For example, in equation 1 of the paper attached to my initial post, the equations for the components in the first layer are:

∂x_i / ∂t = D_{x_{i}}∇^2 x_i + \alpha (x_j - x_i) + f(x_i, y_i)
∂y_i / ∂t = D_{y_{i}}∇^2 y_i + \beta (y_j - y_i) + g(x_i, y_i)

where f() and g() are the Brusselator functions

f(x_i, y_i) = a - (1 - b)x + x^2y
g(x_i, y_i) = bx - x^2y

These components interact with one another in one layer, according to the first (Laplacian) and third (Brusselator) elements in the RHS, and interact with the second layer according to the second element (diffusion) on the RHS.

So I am defining c1[], c2[], etc as the functions for x1, x2, etc. It looks to me like the requirement in Mathematica is to specify any entry c1[], c2[] as a function of *all* the variables, but based on the equations, that does not make sense to me...for example, f() and g() are only functions of the variables in the local layer, so this is why I am confused.