# Why does this system of PDEs take a long time to evaluate?

I am trying to simulate the growth of cell populations through a system of coupled PDEs. There are six equations in total but Mathematica takes a long time to solve the system (+6 hours running) and I do not know when it will stop or end.

The code of my system of equations is the following:

H[\[Psi]_, \[Sigma]_] := Which[\[Psi] <= 0, 1, \[Psi] > \[Sigma], 0, True, 1
- \[Psi]/\[Sigma]]

\[CapitalSigma][Y_, psi0_] := ((Y - psi0)/(1 - Y))*((2 - psi0)/(1 - psi0))

gaman = 0.746;  (*reproduccion celulas sanas*)
gamat = 0.97;    (*reproducion celulas tumorales*)
psi0 = 0.75;      (*razon de volumen libre de estres*)
psin = 0.6;        (*valor umbral que detiene el crecimiento c. sanas*)
psit = 0.8;        (*valor umbral que detiene el crecimiento c. \tumolaes*)
deltat = 0.03;  (*rapidez muerte celulas tumorales*)
deltan = 0.1;    (*rapidez muerte celulas sanas*)
mun = 0.1;           (*rapidez reproduccion MEC por celulas sanas*)
mut = 0.05;        (*rapidez reproduccion MEC por celulas tumorales*)
nu = 0.000016; (*coeficiente de degradacion de la MEC por enzimas*)
pin = 6000000;(*produccion de enzimas que degradan la MEC por c. \sanas*)
pit = 3000000;(*produccion de enzimas que degradan la MEC por c.\tumorales*)
tau = 0.005;         (*tiempo de vida medio de las enzimas*)
fo = 0.25;          (*suministro de oxigeno*)
fg = 0.16;          (*suministro de glucosa*)
betan = 1.2;      (*tasa de absorcion de nutrientes celulas \sanas*)
betat = 1.5;      (*tasa de absorcion de nutrientes celulas \tumorales*)

\[Sigma] = 0.2;
\[CapitalKappa]n = 0.1;
\[CapitalKappa]t = 0.3;
\[Kappa] = 0.00000734;
Diff = 1.0;

eqns = {
D[fin[t, x, y], t] == fin[t, x, y]*\[CapitalKappa]n*(D[fin[t, x,y]*\
[CapitalSigma][
fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0], x, x] +

D[fin[t, x, y]*\[CapitalSigma][fin[t, x, y] + fit[t, x, y] + mat[t,
x, y], psi0], y,y]) gaman*fin[t, x, y]*H[(fin[t, x, y] + fit[t, x, y]
+ mat[t, x, y]) - psin, \[Sigma]] - (1 - oxi[t, x, y])*deltan*fin[t,
x, y],

D[fit[t, x, y], t] == fit[t, x,y]*\[CapitalKappa]t*(D[fit[t, x, y]*\
[CapitalSigma][fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0], x,
x] + D[fit[t, x, y]*\[CapitalSigma][fin[t, x, y] + fit[t, x, y] +
mat[t, x, y], psi0], y, y]) + gamat*fit[t, x, y]*H[(fin[t, x, y] +
fit[t, x, y] + mat[t, x, y]) - psit, \[Sigma]] - (1 - glu[t, x,
y])*deltat*fit[t, x, y],

D[mat[t, x, y], t] == mun*fin[t, x, y]*\[CapitalSigma][fin[t, x, y] +
fit[t, x, y] + mat[t, x, y], psi0] + mut*fit[t, x, y]*\[CapitalSigma]
[fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0] - nu*enz[t, x,
y]*mat[t, x, y],

D[enz[t, x, y],t] == \[Kappa]*(D[enz[t, x, y], x, x] + D[enz[t, x,
y], y, y]) + pin*fin[t, x, y]*\[CapitalSigma][fin[t, x, y] + fit[t,
x, y] + mat[t, x, y], psi0] + pit*fit[t, x, y]*\[CapitalSigma][
fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0] - enz[t, x, y]/tau,

D[oxi[t, x, y], t] == Diff*(D[oxi[t, x, y], x, x] + D[oxi[t, x, y],y,
y]) - betan*fin[t, x, y]*oxi[t, x, y] + fo,

D[glu[t, x, y], t] == Diff*(D[glu[t, x, y], x, x] + D[glu[t, x, y],y,
y]) - betat*fit[t, x, y]*glu[t, x, y] + fg};

IC = {fin[0, x, y] == 0.45, fit[0, x, y] == 0.005,
mat[0, x, y] == 0.2, enz[0, x, y] == 0.3, oxi[0, x, y] == 0.25,
glu[0, x, y] == 0.16};
BC = {
Derivative[0, 1, 0][fin][t, 0, y] == 0,
Derivative[0, 1, 0][fin][t, 1, y] == 0,
Derivative[0, 0, 1][fin][t, x, 0] == 0,
Derivative[0, 0, 1][fin][t, x, 1] == 0,

Derivative[0, 1, 0][fit][t, 0, y] == 0,
Derivative[0, 1, 0][fit][t, 1, y] == 0,
Derivative[0, 0, 1][fit][t, x, 0] == 0,
Derivative[0, 0, 1][fit][t, x, 1] == 0,

Derivative[0, 1, 0][mat][t, 0, y] == 0,
Derivative[0, 1, 0][mat][t, 1, y] == 0,
Derivative[0, 0, 1][mat][t, x, 0] == 0,
Derivative[0, 0, 1][mat][t, x, 1] == 0,

Derivative[0, 1, 0][enz][t, 0, y] == 0,
Derivative[0, 1, 0][enz][t, 1, y] == 0,
Derivative[0, 0, 1][enz][t, x, 0] == 0,
Derivative[0, 0, 1][enz][t, x, 1] == 0,

Derivative[0, 1, 0][oxi][t, 0, y] == 0,
Derivative[0, 1, 0][oxi][t, 1, y] == 0,
Derivative[0, 0, 1][oxi][t, x, 0] == 0,
Derivative[0, 0, 1][oxi][t, x, 1] == 0,

Derivative[0, 1, 0][glu][t, 0, y] == 0,
Derivative[0, 1, 0][glu][t, 1, y] == 0,
Derivative[0, 0, 1][glu][t, x, 0] == 0,
Derivative[0, 0, 1][glu][t, x, 1] == 0};

solu = NDSolve[{eqns, IC, BC}, {fin, fit, mat, enz, oxi, glu}, {x, 0,
1}, {y, 0, 1}, {t, 0, 80}];


I want to see how cell populations grow in 2D but this it take long time to evaluate the system. My boundary conditions are arbitrary so i don't know if they are correct.

I attach an image of my system of equations that I try to solve:

UPDATE I change my BC and my Mathematica finally solve the system quickly but when using DensityPlot to see how behaves in spacetime is not what I really expected.

Really what i hope is see how evolve the volume ratios occupied by two types cells represented by $\phi_n$ and $\phi_t$ (in the image, in the code are "fin" and "fit") into a mesh. In the mesh i want to represent a blood vessel (a circle in a middle of the mesh, for example) and from there to see how the initial group of cells evolve or expand in the confined space of the mesh.

I update my code and I put it below:

H[\[Psi]_, \[Sigma]_] := Which[\[Psi] <= 0, 1, \[Psi] > \[Sigma], 0, True, 1
- \[Psi]/\[Sigma]]

\[CapitalSigma][Y_, psi0_] := ((Y - psi0)/(1 - Y))*((2 - psi0)/(1 - psi0))

gaman = 0.746;
gamat = 0.97;
psi0 = 0.75;
psin = 0.6;
psit = 0.8;
deltat = 0.03;
deltan = 0.1;
mun = 0.1;
mut = 0.05;
nu = 0.000016;
pin = 6000000;
pit = 3000000;
tau = 0.005;
fo = 0.25;
fg = 0.16;
betan = 1.2;
betat = 1.5;
\[Sigma] = 0.2;
\[CapitalKappa]n = 0.1;
\[CapitalKappa]t = 0.3;
\[Kappa] = 0.00000734;
Diff = 1;
L = 5;
T = 80;

eqns = {
D[fin[t, x, y], t] ==
fin[t, x,
y]*\[CapitalKappa]n*(D[
fin[t, x, y]*\[CapitalSigma][
fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0], x, x] +
D[fin[t, x, y]*\[CapitalSigma][
fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0], y, y]) +
fin[t, x, y]*gaman*fin[t, x, y]*
H[(fin[t, x, y] + fit[t, x, y] + mat[t, x, y]) -
psin, \[Sigma]] - deltan*fin[t, x, y],

D[fit[t, x, y], t] ==
fit[t, x,
y]*\[CapitalKappa]t*(D[
fit[t, x, y]*\[CapitalSigma][
fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0], x, x] +
D[fit[t, x, y]*\[CapitalSigma][
fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0], y, y]) +
fit[t, x, y]*gamat*fit[t, x, y]*
H[(fin[t, x, y] + fit[t, x, y] + mat[t, x, y]) -
psit, \[Sigma]] - deltat*fit[t, x, y],

D[mat[t, x, y], t] ==
mun*fin[t, x, y]*\[CapitalSigma][
fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0] +
mut*fit[t, x, y]*\[CapitalSigma][
fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0] -
nu*enz[t, x, y]*mat[t, x, y],

D[enz[t, x, y],
t] == \[Kappa]*(D[enz[t, x, y], x, x] + D[enz[t, x, y], y, y]) +
pin*fin[t, x, y]*\[CapitalSigma][
fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0] +
pit*fit[t, x, y]*\[CapitalSigma][
fin[t, x, y] + fit[t, x, y] + mat[t, x, y], psi0] -
enz[t, x, y]/tau,

D[oxi[t, x, y], t] ==
Diff*(D[oxi[t, x, y], x, x] + D[oxi[t, x, y], y, y]) -
betan*fin[t, x, y]*oxi[t, x, y] + fo,

D[glu[t, x, y], t] ==
Diff*(D[glu[t, x, y], x, x] + D[glu[t, x, y], y, y]) -
betat*fit[t, x, y]*glu[t, x, y] + fg};

IC = {
fin[0, x, y] == 0.45,
fit[0, x, y] == 0.05,
mat[0, x, y] == 0.2,
enz[0, x, y] == 0.3,
oxi[0, x, y] == 0.25,
glu[0, x, y] == 0.16};
BC = {
fin[t, -L, y] == fin[t, L, y],
fin[t, x, -L] == fin[t, x, L],
fit[t, -L, y] == fit[t, L, y],
fit[t, x, -L] == fit[t, x, L],
mat[t, -L, y] == mat[t, L, y],
mat[t, x, -L] == mat[t, x, L],
enz[t, -L, y] == enz[t, L, y],
enz[t, x, -L] == enz[t, x, L],
oxi[t, -L, y] == oxi[t, L, y],
oxi[t, x, -L] == oxi[t, x, L],
glu[t, -L, y] == glu[t, L, y],
glu[t, x, -L] == glu[t, x, L]
};

pdes = Flatten@{eqns, IC, BC};

sol1 = NDSolve[
pdes, {fin, fit, mat, enz, oxi, glu}, {t, 0, T}, {x, -L, L}, {y, -L,
L}, PrecisionGoal -> 2,
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 100}}]

Plot[Evaluate[fin[t, 1, 1] /. sol1], {t, 0, 80}, PlotRange -> All]

Plot[Evaluate[mat[t, 1, 1] /. sol1], {t, 0, 80}, PlotRange -> All]

Table[DensityPlot[
Evaluate[glu[t, x, y] /. sol1], {x, -L, L}, {y, -L, L},
ColorFunction -> "SunsetColors", PlotLegends -> Automatic], {t, 0,
80, 10}]


Any suggestions are very appreciated.

• Have you tried solving it in 1D first, to see if that even works? Commented Dec 7, 2017 at 2:49
• I changed the boundary conditions and could already solve the system. It takes approximately five minutes to solve. Thanks you all. Commented Dec 8, 2017 at 4:49
• I was looking at your profile and I realized that working with reaction-diffusion equations and biological populations. I working with cell populations and i have some questions about my system. I already solved the problem with BC but now I need to solve how those cell populations evolve over time arround a blood vessel (represented graphically as a circle) into a mesh. I used DensityPlot for to seehow the system of equations behaves in space and time but it is not what I expected and i don't know how represent a blood vessel (a circle, for example) into a mesh or in DensityPlot. Commented Dec 11, 2017 at 22:34
• I update my post Commented Dec 11, 2017 at 23:24
• Hi Rodrigo - most of my problems don't have a complex geometry underneath but I'll look at your edit. Commented Dec 11, 2017 at 23:54

NDSolve[{eqns, IC, BC}, {fin, fit, mat, enz, oxi, glu}, {x, 0, 1}, {y,