5
$\begingroup$

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
$\endgroup$
  • 4
    $\begingroup$ Julia's differential equations capabilities may be worth looking into (link). The language is intentionally similar to Matlab. $\endgroup$ – Richard Nov 22 '19 at 21:32
  • $\begingroup$ 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 $\endgroup$ – whpowell96 Nov 23 '19 at 4:55
  • 1
    $\begingroup$ Can you post your MATLAB code? $\endgroup$ – Richard Nov 23 '19 at 17:08
  • $\begingroup$ 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 $\endgroup$ – AleManara Nov 23 '19 at 18:40
  • $\begingroup$ Thank's all for the answers, i'm trying to understand $\endgroup$ – AleManara Nov 23 '19 at 18:49

Your Answer

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

Browse other questions tagged or ask your own question.