2
$\begingroup$

I'm trying to model the diffusion of a ligand that is binding to a receptor (similar to the enzyme substrate model) where the reaction is

enter image description here

I tried using this question as a reference. I'm mainly interested in the rate of production of $ L $ as $ L[x]=(k_1*L[x]*r[x])-(d_1*[R]) $

This is my code so far

ClearAll[r, L, R, x];
reactions = {r + L -> R, R -> r + L};
vars = {r, L, R};
rates = {0.00193, 0.007};
init = <|r -> 1.7, L -> 1, R -> 0|>;
n = (200*^9)/(1.3*^-9);
NA = 6.02*^23;
NT = ((n/NA)*1930)*1.7;

det = NDSolveValue[{r'[x] == 0.007 R[x] - 0.2 r[x] L[x],
    L'[x] == 0.007 R[x] - 0.2 r[x] L[x],
    R'[x] == 0.2 r[x] L[x] - 0.007 R[x],
    r[0] == NT, L[0] == 1, R[0] == 0}, vars, {x, 0, 100}];

Plot[Evaluate@Through@det@x, {x, 0, 100}]

(*----------------------------------------------------------*)

diffCo = 0.0001;

bc = {DirichletCondition[u[x] == 1, x == 0], 
   DirichletCondition[u[x] == 0, x == 25]};

eqn = diffCo*u''[x] + L'[x] == 0;

(*Analytic Solution*)
solDSolve = u[x] /. First@DSolve[{eqn, bc}, u[x], {x, 0, 25}]

(*Numeric Solution*)
solNDSolve = 
 NDSolve[{eqn, bc}, u, {x, 0, 25}, 
  Method -> {"FiniteElement", MeshOptions -> MaxCellMeasure -> 0.001}]

Grid[{{Plot[solDSolve, {x, 0, 1}, PlotRange -> Full, ImageSize -> 300,
     PlotLabel -> "DSolve - Analytic Solution"]}, {Plot[
    Evaluate[u[x] /. solNDSolve], {x, 0, 1}, PlotRange -> Full, 
    ImageSize -> 300, PlotLabel -> "NDSolve - Numeric Solution"]}}, 
 Frame -> All, Spacings -> {1, 1}, Alignment -> Center]

I'm not too sure how to go about combining the solution for L[x] that I get in the top to the diffusion equation that I have later on.

If anyone could help me, I would greatly appreciate it.

$\endgroup$
  • 2
    $\begingroup$ I would like to shortly mention that while your equation, indeed, describes a version of the ligand-receptor interaction, the one for the Michaelis-Menten reaction is usually written in an essentially different form. For the Michaelis-Menten one can have a look here: J. D. Murray, Lectures on nonlinear differential equation model in biology. (Clarendon Press, Oxford, 1977), while for various forms of the ligand-receptor - here L. A. Segel, Biological kinetics. (Cambridge University Press, Cambridge, 1991). $\endgroup$ – Alexei Boulbitch Jun 12 '18 at 8:15
5
$\begingroup$

I deduce from the few that I know about chemistry that the first ODE (the one whose solution is stored in det) describes the reaction under the assumption of infinite diffusion, i.e., under the assumption that the concentrations at each time instance are constant in the whole medium (I interpret x as time variable here.

If you are going to couple that to a diffusion, then the concentrations have to be functions of both time and space and you will get a system of reaction-diffusion equations, so the equations should look more like this:

xx = {x, y, z};
sx = Sequence @@ xx;
{
 D[r[sx, t], t] == dr Laplacian[r[sx, t], xx] + 0.007 R[sx, t] - 0.2 r[sx, t] L[sx, t],
  D[L[sx, t], t] == dL Laplacian[L[sx, t], xx] + 0.007 R[sx, t] - 0.2 r[sx, t] L[sx, t],
  D[R[sx, t], t] == dR Laplacian[R[sx, t], xx] + 0.2 r[sx, t] L[sx, t] - 0.007 R[sx, t]
 }

with appropriate boundary conditions: an initial condition for each concentration and probably homogeneous Neumann conditions (if you have a closed system). Here, dr, dL and dR are the diffusivities of the reagents, t is the time variable and {x,y,z} is a three-dimensional space variable. For the one-dimensional case, just use xx = {x}.

If you are looking for the steady state/equilibrium solution, you can drop the dependencies on t and set the left hand sides to 0:

xx = {x, y, z};
sx = Sequence @@ xx;
{
 0 == dr Laplacian[r[sx], xx] + 0.007 R[sx] - 0.2 r[sx] L[sx], 
 0 == dL Laplacian[L[sx], xx] + 0.007 R[sx] - 0.2 r[sx] L[sx], 
 0 == dR Laplacian[R[sx], xx] + 0.2 r[sx] L[sx] - 0.007 R[sx]
 }

along with boundary conditions of your choice.

So, a complete system with Dirichlet boundary conditions could look like this:

dr = 0.001; dL = 0.001; dR = 0.001; 
bc1 = {
   DirichletCondition[r[x] == 0.8, x == 0], 
   DirichletCondition[L[x] == 1, x == 0], 
   DirichletCondition[R[x] == 0, x == 0], 
   DirichletCondition[r[x] == 0, x == 100],
   DirichletCondition[L[x] == 0, x == 100],
   DirichletCondition[R[x] == 0.8, x == 100]
   };
xx = {x};
sx = Sequence @@ xx;
eqn = {
   0 == dr Laplacian[r[sx], xx] + 0.007 R[sx] - 0.2 r[sx] L[sx], 
   0 == dL Laplacian[L[sx], xx] + 0.007 R[sx] - 0.2 r[sx] L[sx], 
   0 == dR Laplacian[R[sx], xx] + 0.2 r[sx] L[sx] - 0.007 R[sx]
   } ;
sol = NDSolve[Join[eqn, bc1], {r[sx], L[sx], R[sx]}, {x, 0, 100}]

Unfortunately, Mathematica reminds us that NDSolve's FEM solver cannot solve nonlinear PDEs, yet.

A handwoven semi-implicit transient solver (experimental)

First some helper functions

(* Computes load vector for the reaction-diffusion system *)
computeLoad[pat_, celldata_, r_, R_, L_] := cAssembleDenseVector[
   pat,
   Flatten[Join[
     cLoad[celldata, 0.007 R, -0.2 r, L],
     cLoad[celldata, 0.007 R, -0.2 r, L],
     cLoad[celldata, -0.007 R, 0.2 r, L]
     ]],
   {3 Length[r]}
   ];

(* Computes load vector for function (f1 + f2 * f3) *)
cLoad = Block[{PP, P, f1, f2, f3, UU, VV, WW, u, v, w, load},
   PP = Table[Compile`GetElement[P, i, 1], {i, 1, 2}];
   UU = Table[Compile`GetElement[f1, i], {i, 1, 2}];
   VV = Table[Compile`GetElement[f2, i], {i, 1, 2}];
   WW = Table[Compile`GetElement[f3, i], {i, 1, 2}];
   u = t \[Function] (1 - t) UU[[1]] + t UU[[2]];
   v = t \[Function] (1 - t) VV[[1]] + t VV[[2]];
   w = t \[Function] (1 - t) WW[[1]] + t WW[[2]];
   load = Abs[PP[[2]] - PP[[1]]] Integrate[(u[t] + v[t] w[t]) {(1 - t), t}, {t, 0, 1}];

   With[{code = N[load]},
    Compile[{{P, _Real, 2}, {f1, _Real, 1}, {f2, _Real, 1}, {f3, _Real, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
    ]
   ];

cAssembleDenseVector = 
  Compile[{{ilist, _Integer, 1}, {values, _Real, 1}, {dims, _Integer, 1}},
   Block[{A},
    A = Table[0., {i, 1, Compile`GetElement[dims, 1]}];
    Do[A[[Compile`GetElement[ilist, i]]] += Compile`GetElement[values, i], {i, 1, Length[values]}];
    A
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Next, we intitialize the finite element method and try to use as much built-in capability as possible.

Needs["NDSolve`FEM`"]

ν = 0.01;
dr = ν; dL = ν; dR = ν;
xx = {x};
sx = Sequence @@ xx;

bc = {
   DirichletCondition[r[x] == 0.8, x == 0],
   DirichletCondition[L[x] == 1, x == 0],
   DirichletCondition[R[x] == 0, x == 0],
   DirichletCondition[r[x] == 0, x == 100],
   DirichletCondition[L[x] == 0, x == 100],
   DirichletCondition[R[x] == 0.8, x == 100]
   };

(*Initialization of Finite Element Method*)
reg = ToElementMesh[
   DiscretizeRegion[Line[{{0}, {100}}]],
   "MeshOrder" -> 1,
   MaxCellMeasure -> 2
   ];
vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{r, R, L}, {x}}];
sd = NDSolve`SolutionData[{"Space"} -> {reg}];
cdata = InitializePDECoefficients[vd, sd,
   "DiffusionCoefficients" -> -DiagonalMatrix[{dr, dR, dL}],
   "MassCoefficients" -> IdentityMatrix[3]
   (*"LoadCoefficients"\[Rule]{{Sin[2 x]+Cos[x+3 y]}}*)
   ];
bcdata = InitializeBoundaryConditions[vd, sd, bc];
mdata = InitializePDEMethodData[vd, sd];

(*Discretization*)
dpde = DiscretizePDE[cdata, mdata, sd];
dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd];
{load, stiffness, damping, mass} = dpde["All"];
DeployBoundaryConditions[{load, stiffness}, dbc];

cells = reg["MeshElements"][[1, 1]];
n = Length[reg["Coordinates"]];
m = Dimensions[cells][[2]];
celldata = Partition[reg["Coordinates"][[Flatten[cells]]], m];
pat = Flatten[Join[cells, cells + n, cells + 2 n]];

Setting up the time integrator. We use $\theta$-method with θ = 0.8 and time stepsize τ = 0.15.

θ = .8;
τ = 0.15;
aplus = mass + (τ θ) stiffness;
DeployBoundaryConditions[{load, aplus}, dbc];
sol = LinearSolve[aplus, Method -> "Banded"];
aminus = mass + τ (1. - θ) stiffness;
DeployBoundaryConditions[{load, aminus}, dbc];
step[u_] := Module[{load},
  load = SparseArray[
    Partition[
     aminus.u + 
      computeLoad[pat, celldata, Sequence @@ Partition[u, n]], 1]];
  DeployBoundaryConditions[{load, stiffness}, dbc];
  sol[Flatten[load]]
  ]

Setting up the initial conditions and computing 10000 time iterations.

r0 = 0.2 + 0.1 RandomReal[{-1, 1}, n];
R0 = 0.2 + 0.1 RandomReal[{-1, 1}, n];
L0 = 0.2 + 0.1 RandomReal[{-1, 1}, n];
data = NestList[step, Join[r0, R0, L0], 10000];
{rlist, Rlist, Llist} = Partition[data, {Length[data], n}][[1]];

A plot of the solutions

Manipulate[
 ListLinePlot[{rlist[[i]], Rlist[[i]], Llist[[i]]},
  PlotRange -> {-1, 1} 4,
  PlotLegends -> {"r", "R", "L"},
  AxesLabel -> {"x", "Concentration"}
  ],
 {i, 1, Length[rlist], 1}
 ]

enter image description here

The result is somewhat odd. I guess something must be wrong about the reaction part of the equations... (Of course, bugs can be also contained in other parts.)

$\endgroup$
  • $\begingroup$ Thank you so much for your response! So, for this situation, I would need a boundary conditions that accounts for r, L, and R correct? For the steady-state case, I know for the boundary conditions, I was thinking at x=0, r=0.8, L=1, R=0 while at x=100, r=0, L=0, R=0.8. Could the Dirichlet Condition be used in this situation? Or I was also thinking about the Neumann Condition where at x=100, the flux/derivative of everything is 0. I'm a bit lost in how to go about setting up on how to set up DSolve/NDSolve when there's multiple equations. $\endgroup$ – AhWee Jun 11 '18 at 8:28
  • $\begingroup$ I tried setting it up, but I don't believe I went about it the right way since I got back multiple errors. dr = 0.001 dL = 0.001 dR = 0.001 bc1 = {DirichletCondition[r[0] = 0.8, L[0] = 1, R[0] = 0], DirichletCondition[r[100] = 0, L[100] = 0, R[100] = 0.8]} xx = {x, y, z}; sx = Sequence @@ xx; eqn = {0 == dr Laplacian[r[sx], xx] + 0.007 R[sx] - 0.2 r[sx] L[sx], 0 == dL Laplacian[r[sx], xx] + 0.007 R[sx] - 0.2 r[sx] L[sx], 0 == dR Laplacian[R[sx], xx] + 0.2 r[sx] L[sx] - 0.007 R[sx]} sol = NDSolveValue[{eqn, bc1}, {r, L, R}, {x, 0, 100}] $\endgroup$ – AhWee Jun 11 '18 at 8:36
  • $\begingroup$ So, you are interested in the steady-state solution? In that case, homogeneous Neumann conditions would imply that r, L, R are constant over the domain. Dirichlet conditions are somewhat artificial in this context since it is in general not possible to fix the concentration at the boundary. One could use inhomogenous Neumann conditions for modelling inflow and outflow of reagents or Robin boundary conditions for controlled in/outflow depending on concentration. $\endgroup$ – Henrik Schumacher Jun 11 '18 at 8:36
  • $\begingroup$ You need to use == (a.k.a Equal) instead if = (a.k.a. Set) for specifying equations. $\endgroup$ – Henrik Schumacher Jun 11 '18 at 8:37
  • $\begingroup$ I wanted to try getting steady state to work first and then maybe afterwards see if I can get the time dependent examples to work. My overall interest is to see how far L will travel, given that it's diffusing in 1D with the reaction rate consuming it. My initial thought is to set in the beginning (x=0) an initial concentration and setting the endpoint to be a significant distance, where L has depleted to 0 and to pinpoint at what distance this occurs. So, Dirichlet conditions wouldn't work in this case? (Ohh. == in the portion where I had my boundary conditions correct?) $\endgroup$ – AhWee Jun 11 '18 at 8:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.