Numerical Solution of diffusion PDE with varying boundary conditions.

First of all got to clarify that I have a biology background, so some things might seem obvious to the average math person. I have searched a lot about a solution to this and I only come across simple, constant boundary conditions which is not the case in my example. Also, there is no analytical solution possibility. I am thoroughly explaining the problem in purpose, from a biologist's point of view.

Background: I am trying to simulate a system of coupled ODEs and PDEs concerning the diffusion of two molecules between two different cells. For the sake of simplicity, let's just be concerned with the diffusion PDE of a single molecule from one cell to another.

Assume that the concentration of this molecule inside the cells and in the extracellular material is expressed as $$u_a$$, $$u_b$$ and $$u_e$$ respectively. Cells are assumed to be just points of production and they are separated by a distance $$L$$. Diffusion in this monodimensional distance is governed by a diffusion constant $$\Theta$$. There is also a rate of release/absorption from the inside of the cell to the outside, due to cell membrane separation, which is expressed by $$\eta$$.

The rate of change for quantities $$u_a$$ and $$u_b$$ is known, in terms of ordinary differential equations, and depends on the synthesis(if present), degradation, and absorption/release of this quantity to the extracellular material. Specifically:

$$\frac{du_a}{dt}(t) = K_Q s(t) + \eta \Big(u_e(t) - u_a(t) \Big) - \gamma_i u_a(t) \quad : \textrm{a cell}$$ $$\frac{du_b}{dt}(t) = \eta \Big(u_e(t) - u_b(t) \Big) - \gamma_i u_b(t) \quad : \textrm{b cell}$$

where $$K_Q$$ is a constant, $$s(t)$$ is just the concentration of another activating species (which itself is governed by another ODE), $$\gamma_i$$ the intracellular degradation constant and $$\eta$$ the cell wall diffusion constant ($$\eta = 2\ min^{-1}$$). Cell wall is assumed as point "boundary wall", that just allows molecules to pass with this $$\eta$$ coefficient, depending on their inside/outside differences. In cell $$b$$, there isn't any production of the molecule - just in/out flow and internal degradation.

Problem: I want to model the diffusion of this substance throughout the monodimensional domain (line) between the two cells and I have a problem applying the correct boundary conditions. I consider the boundary conditions as the points right outside the cells. Initial conditions are known (assumed to be zero everywhere: $$u_a(0)=0$$, $$u_b(0)=0$$ and $$u_e(x_i,0) = 0$$.).

As known, diffusion of a substance in one dimension is governed by the general equation:

$$\frac{\partial{u_e}}{\partial{t}} = \Theta \frac{\partial^2{u_e}}{\partial{x^2}} + g(x,t)$$

where $$\Theta$$ is the constant diffusion coefficient and $$g(x,t)$$ a term of addition or removal of concentration in each point. It differs for the boundary points.

If we discretize the line of length $$L$$ between the cells in $$N$$ points, then we get three cases: a set of $$N-2$$ ODEs for the intermediate points $$\frac{du_e^i}{dt}$$, $$i \in \{1,...,N-1\}$$, and two boundary conditions for points $$i=0$$ and $$i=N$$ (right outside the cells a and b. My problem is that I don't understand what form should these boundaries take. They should take into account the intracellular and extracellular concentration of the molecule ($$u_a$$, $$u_b$$) in order to calculate how much substance is diffused outwards or inwards the cells. I understand that I cannot use Dirichlet boundaries since this concentration is varying through time. I tried using Neumann boundaries (i.e. setting $$\frac{du_e^0}{dt}(x_0,t)$$ and $$\frac{du_e^N}{dt}(x_N,t)$$ instead of $$u(x_0,t)$$ and $$u(x_N,t)$$). Should I include diffusion to the extracellular material and degradation of the substance in these boundaries? If I do this, I obtain the following system:

$$\begin{cases} \begin{gather*} \frac{du_e^i}{dt} = \frac{\Theta}{\Delta x^2} \Big(u_e^{i+1}(t) - 2u_e^i(t) + u_e^{i-1}(t) \Big) - \gamma_e u_e^{i}(t) \textrm{ , for } i \in \{1,...,N-1\} \\[1.5em] \frac{du_e^0}{dt} = \eta \Big( u_a (t) - u_e^0(t) \Big) + \frac{\Theta}{\Delta x^2} \Big(2u_e^1(t) - 2u_e^0(t) \Big) - \gamma_e u_e^0(t) \\[1.5em] \frac{du_e^N}{dt} = \eta \Big( u_b (t) - u_e^N (t) \Big) + \frac{\Theta}{\Delta x^2} \Big(2u_e^{N-1}(t) - 2u_e^N(t) \Big) - \gamma_e u_e^N(t) \end{gather*} \end{cases}$$

which doesn't converge - concentrations in the extracellular material tend to $$\pm \infty$$.

The alternative would be to ommit the extracellular diffusion and degradation terms for $$i=0$$ and $$i=N$$, but would this be physically accurate?

For solving this system, I use a low-level Forward-Euler solver I coded into python. I always follow the necessary $$dt \leq \frac{\Delta x^2}{2 \Theta}$$ for the timestep.

Goals: Would like to get some insight as to why my extracellular concentrations do not converge: is it a problem of the finite differences method, even when following the restriction, or of the boundary conditions I'm using? I would prefer to solve the system using the simple methods I have already coded, rather than trying more elaborate techniques (implicit solution or Crank-Nicolson scheme).

P.S. My discretization is based on this life-saving manual: http://hplgit.github.io/prog4comp/doc/pub/p4c-sphinx-Python/._pylight006.html

EDIT: for clarification, I have added the exact equations for cells $$a$$ and $$b$$ concerning the quantity $$u$$. The question then boils down to how to connect extracellular $$u$$ with the boundary point $$u$$? For example, at point $$i=N$$, the diffused $$u_e$$ all the way from cell $$a$$ needs to be passed on to $$\frac{du_e^N}{dt}$$ and this would involve the inclusion of a term considering $$u_e^{N-1}$$ - if not, the rate of change will always be zero and the amount constant. In what way to include this concentration/point is what I cannot find out.

EDIT2 Adding the full set of equations after request. Diffused quantities are $$Q_j^1$$ and $$Q_j^2$$ molecules (correspond to $$u$$), where index $$j$$ defines the locality of these quantities: $$c$$ for controller cell, $$p$$ for producer cell and $$e$$ for extracellular material. All values in $$[x]$$ refer to concentrations of species $$x$$, and correspond to functions of time: $$[x](t)$$. $$[Ref](t)$$ is assumed constant.

$$\begin{gather} \textrm{Cell C}=\begin{cases} \begin{split} \frac{d[A\!:\!Q^2]}{dt} & = \Big( \chi_{A:Q_2,r,0} + \chi_{A:Q_2,r} \frac{K_{r}^{n_r}}{K_{r}^{n_r} + [Ref]^{n_r}} \Big) \cdot \\ & \Big( \chi_{A:Q_2,a,0} + \chi_{A:Q_2,a} \frac{[Q_c^2]^{n_q}}{K_{q}^{n_q} + [Q_c^2]^{n_q}} \Big) - \gamma_{A:Q_2}[A\!:\!Q^2] \end{split} \\[3ex] \frac{d[B]}{dt} = \chi_{B_0} + \chi_B \frac{[A\!:\!Q^2]^{n_b}}{K_b^{n_b} + [A\!:\!Q^2]^{n_b}} - \gamma_B [B] \end{cases}\\ \textrm{Cell P}=\begin{cases} \frac{d[C]}{dt} = \chi_{C_0} + \chi_C \frac{[Q_p^1]^{n_c}}{K_c^{n_c} + [Q_p^1]^{n_c}} - \gamma_C [C] \\ \frac{d[D]}{dt} = \chi_{D_0} + \chi_D \frac{K_d^{n_d}}{K_d^{n_d}+[C]^{n_d}} - \gamma_D [D] \end{cases} \\ \frac{d[Q_c^1]}{dt} = K_{Q^1} [B] + \eta \Big( [Q_e^1] - [Q_c^1] \Big) - \gamma_i [Q_c^1]\\ \frac{d[Q_p^1]}{dt} = \eta \Big( [Q_e^1] - [Q_p^1] \Big) - \gamma_i [Q_p^1] \\ \frac{\partial [Q_e^1]}{\partial t} = \eta \Big( [Q_c^1] - [Q_e^1] \Big) + \eta \Big( [Q_p^1] - [Q_e^1] \Big) - \gamma_e [Q_e^1] + \Theta \nabla^2[Q_e^1] \\ \frac{d[Q_c^2]}{dt} = \eta \Big([Q_e^2] - [Q_c^2] \Big) - \gamma_i [Q_c^2] \\ \frac{d[Q_p^2]}{dt} = K_{Q^2} [D] + \eta \Big([Q_e^2] - [Q_p^2] \Big) - \gamma_i [Q_p^2] \\ \frac{\partial [Q_e^2]}{\partial t} = \eta \Big([Q_c^2] - [Q_e^2] \Big) + \eta \Big([Q_p^2] - [Q_e^2] \Big) - \gamma_e [Q_e^2] + \Theta \nabla^2[Q_e^2] \\ \end{gather}$$

• You might be able to solve this analytically using the Laplace Transform. See the book The Mathematics of Diffusion, by John Crank. smile.amazon.com/… Typically, you get an infinite series solution. You can fit that solution to experimental data to get more information about the system than you would if you just hand the PDE over to a solver. – Adrian Keister Aug 6 '19 at 16:34
• For a similar problem, e.g., see this paper co-authored by my father, J. C. Keister: heidmann.com/eyedrop/index.htm. – Adrian Keister Aug 6 '19 at 16:39
• How about you first make one step in the ode's for the endpoints, using the extracellular concentration from the previous step. And then solve the 1d diffusion equation with Dirichlet BC, taking the newly computed (through ode's) concentrations at the next time step? It's merely a suggestion. – VorKir Aug 9 '19 at 16:41

My previous answer was incorrect; here's a more careful derivation.

In the inter-cellular space we have as a starting point the law of diffusion $$J = -\Theta \frac{\partial u_e}{\partial x}$$ for flux $$J$$, from which it follows that $$\frac{\partial u_e}{\partial t} = \Theta \Delta u_e - \gamma_e u_e$$ on $$0 < x < L$$ and \begin{align*} - \Theta \frac{\partial u_e}{\partial x}(0) &= g(0)\\ \Theta \frac{\partial u_e}{\partial x}(L) &= g(L) \end{align*} at the boundary. Here $$g(x)$$ is the rate of concentration injection at the boundary, due to diffusion through the cell wall: $$g(0) = \eta \left[u_a - u_e(0)\right]$$ $$g(L) = \eta \left[u_b - u_e(L)\right]$$

The boundary condition on the left can then be discretized as $$\Theta \frac{u_e^1(t+\Delta t) - u_e^0(t + \Delta t)}{\Delta x} = - \eta\left[u_a(t) - u_e^0(t)\right]$$ and likewise for the right boundary.

Here's a C++ implementation, and an animation of what the concentration looks like over time (left-most and right-most points are $$u_a$$ and $$u_b$$), for some arbitrarily-chosen values of the physical constants. Note that you need a quite small time step for the simulation to be stable. Switching to an implicit discretization will dramatically improve stability.

#include <iostream>
#include <iomanip>

int main()
{
const int N = 10;
double L = 1;
double Theta = 1;
double eta = 1;
double KQ = 1;
double gammai = 0.1;
double gammae = 0.1;

// initial conditions
double ua = 0;
double ub = 0;
double ue[N];
for(int i=0; i<N; i++)
ue[i] = 0;

double dx = L / double(N-1);
double dt = 0.001;

int nsteps = 1000;
for(int step=0; step<nsteps; step++)
{
double newua = ua + dt * (KQ + eta * (ue[0] - ua) - gammai * ua);
double newub = ub + dt * (eta * (ue[N-1]-ub) - gammai * ub);
double newue[10];
newue[0] = newue[1] + dx / Theta * eta * (ua-ue[0]);
newue[N-1] = newue[N-2] + dx / Theta * eta * (ub - ue[N-1]);
for(int i=1; i<N-1; i++)
{
newue[i] = ue[i] + dt * (Theta * (ue[i-1] - 2.0 * ue[i] + ue[i+1]) / dx / dx - gammae * ue[i]);
}
ua = newua;
ub = newub;
for(int i=0; i<N; i++)
ue[i] = newue[i];
std::cout << ua << ' ';
for(int i=0; i<N; i++)
std::cout << ue[i] << ' ';
std::cout << ub << std::endl;

}
}


• Thanks so much for the quick response. I see exactly what you mean, though I have an objection and a clarification to make (added to the original question as well), and I think really this is the boiled down question concerning this system. If I use these boundaries, I don't take into account the diffused molecule at the N point: $u$ is not passed from point N-1 to point N with any term, thus yielding a constant 0 concentration at that point (u molecule is not produced into cell b, only diffused in or out from it). – Soter Aug 6 '19 at 18:43
• @Soter You're right. I've revised my answer. – user7530 Aug 6 '19 at 21:25
• Thanks for the edit. I still meet some problems with this solution: (1) in the boundary conditions discretization you propose the term $\frac{\Theta}{\Delta x}$, but this results in ambiguities concerning the units: $\Theta$ is epxressed on $\mu m^2$ /sec. (2) Even with this discretization, the system never evolves into a steady state as it should - all the values keep increasing exponentially. I'm not so sure about the second point though, because I haven't present the complete system, and maybe its due to a different reason. – Soter Aug 7 '19 at 15:40
• @Soter You are right about (1); the time derivative and decay terms vanish when you integrate over the vanishing neighborhood around the boundary points. I don't know about (2); I'm not seeing any issues in my simulations (see code I posted). Remember that the choice to explicitly discretize the Laplacian imposes a severe timestep restriction ($\Delta t \approx \frac{\Delta x^2}{\Theta}$) and if you are seeing instability my first instinct would be to see if the problem persists with a much smaller timestep. – user7530 Aug 7 '19 at 17:00
• For (2), I have implemented the restriction on $\Delta t$ inside the solver, using always a value of $0.9 \cdot \frac{\Delta x^2}{2 \Theta}$. I have tried to lower the step even more, but I'm still getting the same results. I have also tried to increase the discretization points, but after a limit it becomes impossible to compute: my theta value is huge ($\Theta = 800 \mu m^2\ sec^{-1}$) - relating to the length $L=20\mu m$, and this imposes many problems... Concerning (1) I don't see exactly what you mean, but I guess what I thought of is indeed the case? Trying with odespy right now. – Soter Aug 7 '19 at 17:26