Evolutionary dynamics in vascularised tumors, PDE-ODE coupled system

Question

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.

Here i upload the code:

    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