# Evolutionary dynamics in vascularised tumors, PDE-ODE coupled system

I have to solve the following PDE-ODE system

$$\displaystyle{\partial_{t} n = \bigl[a(s) - b(s)(y - h(s))^{2} - d\int_{\mathbb{R}} n \, dy \; + \; \beta \, \partial^2_{yy} n \;} \\\\ \displaystyle{s' = -\int_{\mathbb{R}} \bigl[\eta_{s}(g(y,s))_{+}n \bigl] \, dy - \lambda_{s}s + S0 \;}$$ where $$n = n(t,y)$$ is the distribution of tumor cells depending on time and phenotype state $$y$$, $$s = s(t)$$ the oxygen concentration with intial oxygen concentration $$S_{0} = 6.3996 \times 10^{-7}$$ and $$(.)_{+}$$ the positive part. The phenotypic state is $$y \in \mathbf{R}$$, the time $$t \in [0,+\infty]$$ . The function $$a(s),b(s),h(s)$$ depending on the oxygen concentration are so defined: $$a(s) := \Bigg( \gamma_{s} \frac{s}{\alpha_{s} + s} + \frac{\varphi^{2}}{\varphi + \gamma_{s}\dfrac{s}{\alpha_{s} + s}} \Bigg)$$ $$b(s) := \bigg(\varphi + \gamma_{s}\dfrac{s}{\alpha_{s} + s} \bigg)$$ $$h(s) := \frac{\varphi}{\varphi + \gamma_{s}\dfrac{s}{\alpha_{s} + s}}$$ the function $$g(y,s)$$ depending on phenotype and oxygen concentration as $$g(y,s) := \gamma_{s} \frac{s}{\alpha_{s} + s}(1 - y^{2})$$

The parameters are the following $$\alpha_{s} \ \text{Michaelis-Menten constant of oxygen} \ 1.5 \times 10^{-7}$$ $$d \ \text{Rate of death due to competition for space} \ 2 \times 10^{-13}$$ $$\varphi_{} \ \text{Maximum rate of cell division via anaerobi pathway} \ 1 \times 10^{-5}$$ $$\gamma_{s} \ \text{Maximum rate of cell division via aerobic pathways} \ 1 \times 10^{-4} \$$ $$\eta_{s} \ \text{Conversion factor for cell consumption of oxygen} \ 2 \times 10^{-11} \$$

My question is how to solve this problem numerically (I'm using Matlab). To start, I used a uniform grid for the domain phenotype $$y∈[−2L,2L]$$ with $$L= 2,Ny= 1200$$ points (https://biologydirect.biomedcentral.com/articles/10.1186/s13062-016-0143-4). After that I also discretized time $$t$$ between the starting time $$t=0$$ and the final time $$t_{f} = 5$$ in units of $$10^{4}s$$ with a time step $$dt= 1s$$ (and $$N_{t}$$time step), according with the CFL condiction.

The method for calculating numerical solutions is based on a time splitting scheme between the conservative part and the reaction term. I approximate the diffusion term through finite difference approximation (three point stencil) and we used an explicit finite difference scheme for the reaction term. The scheme that I implemented is the following $$n^{k+1}_{i} = n^{k}_{i} + dt \, \bigl[ a(s^{k}) - b(s^{k})(y_{i} - h(s^{k}))^{2} - d \, dy \, \sum_{j=1}^{N_{y}}(n^{k}_{j})\bigl] \,n^{k}_{i} + \frac{\beta dt}{dy^{2}} \biggl[ n^{k}_{i+1} - 2n^{k}_{i} + n^{k}_{i-1}\biggl]$$ $$s^{k+1} = s^{k} - dt \, dy \, \sum_{j=1}^{Ny}(\bigl[\eta_{s}(g(y_{j},s^{k}))_{+}n^{k}_{j} \bigl]) - dt \, \lambda_{s} \, s^{k} + dt \, S0$$ $$k = 1,...,N_{t}-1, \ i = 2,...,N_{y}-1.$$ But unfortunately the solution isn't correct and oxygen concentration has very abnormal behavior

Here is the documentation on which this work is based https://arxiv.org/abs/1910.08566 This is an easy case in which there is not a spatial dependence.

    clc
clear all
close all

set(0,'DefaultAxesFontName', 'Times New Roman')
set(0,'DefaultAxesFontSize', 18)
set(0,'defaultaxeslinewidth',1)
set(0,'defaultpatchlinewidth',1)
set(0,'defaultlinelinewidth',4)
set(0,'defaultTextInterpreter','latex')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% To solve
%%%%% n_t = beta*n_yy + (a(s) - b(s)*(y - h(s))^2 - chi*rho)*n
%%%%% s_t = ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters
betay = 1e-6; % rate of spontaneous phenot. variations
chi = 2*1e-13; d = chi; % rate of death for competition
phi = 1e-5;
alfac = 2*1e-6;
alfas = 1.5*1e-7;
muc = 4.4*1e-6;
mus = 2*1e-5;
gammac = 1.8*1e-4;
gammas = 1e-4;
etac = 4*1e-11;
etas = 2*1e-11;
lambdac = 2.3*1e-4;
lambdas = 2.78*1e-6;

%% Time, Space and Phenotype
Ly = 4;
Ny = 1200;
t0 = 0;
tf = 5*1e4;

%% Discretization
ysp = linspace(-Ly,Ly,Ny)';
dy = ysp(end-1) - ysp(end-2);
dt = 1;
tsp = t0:dt:tf;
Nt = length(tsp);

%% Oxygen concentrarion at time 0, Oxyeng into system
S0 = 6.3996*1e-7;

%% Def : a,b,h;
a = @(s) gammas*(s./(alfas + s))  + + (phi).^2./(phi + gammas*s./(alfas + s));
b = @(s) phi + gammas*(s./(alfas + s));
h = @(s) (phi)./(phi + gammas*s./(alfas + s));
g = @(s) (gammas*s./(alfas + s)).*(1 - ysp.^2);

%% Numerical approach, Initialization
Rho0 = 1e8; Mu0 = 0.5; Vu0 = 1; % Data from reference
n0 = Rho0*sqrt(Vu0/(2*pi))*exp(-(Vu0/2)*(ysp - Mu0).^2).*ones(1,1)';
R = zeros(Ny,1);    % initialization
nOld = n0;          % initialization
nNew = n0;          % initialization
sOld = S0;
sNew = S0;
sMemory = S0*ones(Nt,1); % To memorize s(t)
%% I-M-V
rho = ones(Nt,1).*Rho0;
mean = ones(Nt,1).*Mu0;
var = ones(Nt,1).*Vu0;
initial = [1/Vu0,Mu0,Rho0];

for kk = 1:Nt-1 %%% TIME
tk = kk*dt; % time counter
%% Explicit-Eu
%%%% Reaction and Integral for n
R = a(sOld) - b(sOld).*(ysp - h(sOld)).^2 - d*dy*sum(nOld);
%%%% New n
nNew(2:Ny-1) = nOld(2:Ny-1) + dt*R(2:Ny-1).*nOld(2:Ny-1) + ...
+ (betay*dt/(dy^2))*(nOld(3:Ny) - 2*nOld(2:Ny-1) + nOld(1:Ny-2));
%%%% Neumann boundary cond
nNew(1) = nNew(2);
nNew(Ny) = nNew(Ny-1);
%%%% rs
rs = etas*max(0,g(sOld));
sNew = sOld - dt*(dy*sum(rs.*nOld)) - dt*lambdas*sOld + dt*S0;

%% Int - Mean - Variance
rho(kk+1) = dy*sum(nNew);
mean(kk+1) = dy*sum(ysp.*nNew)./rho(kk+1);
var(kk+1) = dy*sum(ysp.^2.*nNew)./rho(kk+1) - mean(kk+1).^2;

%% Updating
nOld = nNew;
sOld = sNew;
sMemory(kk+1) = sNew;
end

%% I-M-V with ODE Plot
figure()
subplot(1,3,1)
plot(tsp,rho,'--')
% axis([t0 tf min(min(rho))-.1e7 max(max(rho))+.1e7])
axis square
xlabel('time')
title('$$\rho(t)$$')
subplot(1,3,2)
plot(tsp,mean,'--')
% axis([t0 tf min(min(media))-.1 max(max(media))+.1])
axis square
xlabel('time')
title('$$\mu(t)$$')
subplot(1,3,3)
plot(tsp,var,'--')
% axis([t0 tf min(min(var))-.1 max(max(var))+.1])
axis square
xlabel('time')
title('$$\sigma^{2}(t)$$')
s = sprintf('I-M-V of solution n(y,s) with tf = %d, dt = %f',tf,dt);
suptitle(s)

%% s,a,b,h,g
figure()
subplot(1,5,1)
plot(tsp,sMemory)
title('S(t)')
xlabel('time')
axis square
subplot(1,5,2)
plot(tsp,a(sMemory))
title('a(S(t))')
xlabel('time')
axis square
subplot(1,5,3)
plot(tsp,b(sMemory))
title('b(S(t))')
xlabel('time')
axis square
subplot(1,5,4)
plot(tsp,h(sMemory))
title('h(S(t))')
xlabel('time')
axis square
subplot(1,5,5)
G = (gammas*sMemory./(alfas + sMemory)).*(1 - ysp'.^2);
plot(tsp,G(:,1))
title('g(S(t))')
xlabel('time')
axis square

• Julia's differential equations capabilities may be worth looking into (link). The language is intentionally similar to Matlab. – Richard Nov 22 '19 at 21:32
• I have not double-checked your discretization. I would highly recommend using the method of lines for this problem. Essentially you discretize in every variable except for time, which gives you a big system of ODEs. You can then input this into one of MATLAB's ODE solvers, such as ode45 or ode15s. This is advantageous because the ODE solvers can automatically handle stability problems by adaptively decreasing the step size $dt$ depending on which solver you use – whpowell96 Nov 23 '19 at 4:55
• Can you post your MATLAB code? – Richard Nov 23 '19 at 17:08
• i try now, it's a bit long. I upload on the question the code, i hope it's clear! for any question i'm here – AleManara Nov 23 '19 at 18:40
• Thank's all for the answers, i'm trying to understand – AleManara Nov 23 '19 at 18:49