# NonLinear system for chemotaxis

I want to solve the chemotaxis mode, given by the next non-linear system:

It is taken from Murray's book: equation (11.30) at pag. 408

$$\frac{\partial n}{\partial t} = D \frac{\partial^2 n}{\partial x^2} -\xi_0 \partial_x \Bigl( n \frac{\partial a}{\partial x} \Bigr)$$

$$\frac{\partial a}{\partial t} = hn - ka + D_a \frac{\partial^2 a}{\partial x^2}$$

where $$h,k,D_a,D$$ are just parameters, and $$D_a>D$$ and the domain is $$x \in [-6,6]$$

I decided to take as no flux boundary conditions, i.e. $$\partial_x(n(-6,t))=\partial_x (a(-6,t))=0$$ $$\partial_x(n(6,t))=\partial_x (a(6,t))=0$$

and as initial conditions $$n(0,x)=e^{-x^2}$$ $$a(0,x)=\cos( \pi x)$$

Note that numerically the conditions are compatbile since the exponential is "flat". I know that analytically it's not true.

I integrated up to time $$T=0.1$$ with my own FEM solver (with linear finite elements) and obtain the following, using the parameters $$D = 2 \quad D_a = 5.5 \quad h = 0.5 \quad k = 0.5 \quad \xi_0 = 0.2$$

I'd like to use Mathematica to check my results and to try what comes out by changing some parameters, but I can't understand how to solve a non-linear system like the one above. Could someone show the plot I should obtain with Mathematica, and, if possible, the right code-snippet?

EDIT:

Here is what I obtain, which has the shape of Daniel answer's, which seems to be similar to his one

EDIT:

The pysical principle behind the model is:

The amoebae of the slime mould Dictyostelium discoideum, with density n(x, t), secrete a chemical attractant, cyclic-AMP, and spatial aggregations of amoebae start to form. THe book says ti use zero-flux boundary conditions, and that's fine. But what initial conditions could I use for $$n(x,t)$$ and $$a(x,t)$$ that are physically relevant?

If you use the Finite Element Method, no flux is the default boundary condition, so there is no need to specify. An alternative to Daniel's answer would be:

(* Define parameters *)
l = 6;
tend = 0.1;
parms = {d -> 2, da -> 5.5, h -> 0.5, k -> 0.5, x0 -> 0.2};
(* Create Parametric PDE operators for n and a *)
parmnop =
D[n[t, x], t] - d D[n[t, x], x, x] + x0 D[n[t, x] D[a[t, x], x], x];
parmaop = D[a[t, x], t] - da D[a[t, x], x, x] + k a[t, x] - h n[t, x];
(* Setup PDE System *)
pden = (parmnop == 0) /. parms;
pdea = (parmaop == 0) /. parms;
icn = n[0, x] == Exp[-x^2];
ica = a[0, x] == Cos[π x];
(* Solve System *)
{nif, aif} =
NDSolveValue[{pden, pdea, icn, ica}, {n, a}, {t, 0, tend}, {x, -l,
l}, Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> MaxCellMeasure -> 0.1}}];

(* Display results *)
Manipulate[
Plot[{nif[t, x], aif[t, x]}, {x, -l, l}, PlotRange -> All], {t, 0,
tend}, ControlPlacement -> Top]


• Well, by setting Method -> {"FiniteElement"}, NDSolve uses pure FEM for solving the system, which isn't quite good for IBVP. Just check Plot[D[aif[t,x],x]/.x->l,{t,0,tend}]. Method -> {MethodOfLines,SpatialDiscretization->{"FiniteElement",MeshOptions->MaxCellMeasure->0.1}} would be a better choice. Sep 29, 2020 at 2:50
• I took the liberty to add @xzczd suggestion. It's better to use the method of lines for time dependent PDEs. Sep 29, 2020 at 6:35
• @xzczd Thank you very much! Actually, when I was working the answer, I was using your recommended approach. I got over zealous trying to tighten the code. I just eye-balled the results without a thorough check. Sep 29, 2020 at 16:18
• @user21 Cheers for the edit! Sep 29, 2020 at 16:18
• @TimLaska Tanks for your answer. Just one last point, about the physics of the system (see my last EDIT). I don't know what are some suitable initial conditions for $n(x,t)$ and $a(x,t)$. Do you have any idea? Sep 30, 2020 at 9:40

Here is my code. Unfortunately, at t==0.1, it does not duplicate your result. I hope I did not make a mistake.

eq = {D[n[x, t], t] ==
d  D[n[x, t], {x, 2}] - c0 D[n[x, t] D[a[x, t], x], x],
D[a[x, t], t] == h  n[x, t ] - k a[x, t] + da  D[a[x, t], {x, 2}],
(D[n[x, t], x] /. x -> -6) == 0, (D[a[x, t], x] /. x -> -6) ==
0, (D[n[x, t], x] /. x ->   6) ==
0, (D[a[x, t], x] /. x ->   6) == 0,
n[x, 0] == Exp[-x^2], a[x, 0] == Cos[Pi x]} /. {d -> 2, da -> 5.5,
h -> 0.5, k -> 0.5, c0 -> 0.2};
sol[x_] = {n[x, 0.1], a[x, 0.1]} /.
NDSolve[eq, {n, a}, {t, 0, 0.1}, {x, -6, 6}][[1]]
Plot[sol[x], {x, -6, 6}, PlotRange -> All]


• Thanks for your answer @DanielHuber. I edited my question with the actual plot I obtain. How did you impose noflux boundary conditions in your code? Sep 28, 2020 at 22:14
• @The b.c.s are (D[n[x, t], x] /. x -> -6) == 0, (D[a[x, t], x] /. x -> -6) == 0, (D[n[x, t], x] /. x -> 6) == 0, (D[a[x, t], x] /. x -> 6) == 0, which are straightforward translations of those b.c.s in traditional math notation. If you're having difficulty in understanding them, please check document of D, ReplaceAll (/.), Rule (->) by pressing the F1 key. Sep 29, 2020 at 1:43
• You may specify no flux b.c. in different way: I used e.g: D[f[x,t],x] /. x->6 . But you can also say: Derivative[1,0][f][6,t ] . What is clearer is up to your taste. Note also for function with only one variable you can abbreviate the first derivative by: f'[x] Sep 29, 2020 at 8:52
• @DanielHuber and xzczd, thanks for you answer/comments. I edited my question with a last point, about the physics of the system. I don't know what are some suitable initial conditions for $n(x,t)$ and $a(x,t)$. Sep 30, 2020 at 9:39
• @Vefhug Hi, I do not know enough about the real problem to give a real good answer. However, no flux boundary means, when there are no walls, that the functions at xmin and xmax stays constant. On the other hand, if you introduce position depended parameters, you can model something of a soft or hard wall and drop the condition that derivatives at the ends are zero. But I think this is more complicated than your teacher intend it to be. Sep 30, 2020 at 10:44