Trouble solving system of pdes with pdepe

I am trying to replicate Figure 2 from the Article "Mathematical Modeling of Alzheimer's Disease" by Hao & Friedman (2016). Below is my code to use pdepe to solve the system of 18 pdes from the paper. I keep receiving the error of "In pdepe (line 287)In math435replicate_pdepe (line 135) Warning: Time integration hasfailed. Solution is available at requested time points up to t=0.000000e+00."

main()
function main
% This is the main function, it encapsulates the whole code and allows
% parameters to be definied in the main body and still be accessible in
% other parts of the code
clear all
close all
% Set the parameter values
D_Ao = 4.32e-2; % diffusion coefficient of AbetaO
D_H = 8.11e-2; % Diffusion coefficient of HMGB-1
D_Talpha = 6.55e-2; % Diffusion coefficient of TNF-alpha
D_Tbeta = 6.55e-2; % Diffusion coefficient of TGF-beta
D_I10 = 6.04e-2; %Diffusion coefficient of IL-10
D_P = 1.2e-1; %Diffusion coefficient of MCP-1
lambda_beta_i = 9.51e-6; % production rate of Amyloid beta i
lambda_N = 8e-9; % production rate of amyloid beta O by neuron
lambda_A = 8e-10; % production rate of amyloid beta O by astrocytes
lambda_tau0 = 8.1e-11; % production rate of tau proteins in health
lambda_tau = 1.35e-11;% Production rate of tau proteins by ROS
lambda_F = 1.662e-3; % Production rate of NFT by tau
lambda_AT_alpha = 1.54; %Production/activation rate of astrocytes by TNF-α
lambda_A_AmyloidbetaO = 1.793; %Production/activation rate of astrocytes by amyloid beta O
lambda_AO = 5e-2; % Production rate of amyloid beta O
lambda_H = 3e-5; % Production rate of HMGB-1
lambda_MF = 2e-2; % Production/activation rate of microglias by NFT
lambda_MA = 2.3e-3; % Production/activation rate of microglias by astrocytes
lambda_M1_Tbeta = 6e-3; % Rate of M1 → M2
lambda_M1hat_beta = 6e-4; % Rate of M1 → M2
lambda_Tbeta_M = 1.5e-2; % Production rate of TGF-β by M
lambda_Tbeta_Mhat = 1.5e-2; % Production rate of TGF-β by Mhat
lambda_Talpha_M1 = 3e-2; % Production rate of TNF-α by M1
lambda_Talpha_M1hat = 3e-2; % Production rate of TNF-α by M1hat
lambda_I10_M2 = 6.67e-3;% Production rate of IL-10 by M2
lambda_I10_M2hat = 6.67e-3; % Production rate of IL-10 by M2 hat
lambda_PA = 6.6e-8; % Production rate of MCP-1 by astrocytes
lambda_PM2 = 1.32e-7; % Production rate of MCP-1 by M2
theta = 0.9; % M2 /M1 effectivity in clearance of Aoβ, estimated
alpha = 5; % Flux rate of macrophages, estimated
beta = 10; % Proinflammatory/anti-inflammatory ratio, estimated
gamma = 1; % I10 inhibition ratio, estimated
d_Abeta_i = 9.51; % DegradationrateofAiβ
d_AbetaO = 9.51; % Degradation rate of Aoβ
d_AbetaO_M = 2e-3; % Clearance rate of Aoβ by microglia
d_AbetaO_Mhat = 10e-2; % Clearance rate of Aoβ by macrophages
d_tau = 0.277; % Degradation rate of tau proteins
d_Fi = 2.77e-3; % Degradation rate of intracellular NFT
d_Fo = 2.77e-4; % Degradation rate of extracellular NFT
d_N = 1.9e-4; % Death rate of neurons
d_NF = 3.4e-4; % Death rate of neurons by NFTs
d_NT = 1.7e-4; % Death rate of neurons by TNF-α
d_Nd_M = 0.06; % Clearance rate of dead neurons by M
d_Nd_Mhat = 0.02; % Clearance rate of dead neurons by Mhat
d_A = 1.2e-3; % Death rate of astrocytes
d_M1 = 0.015; % Death rate of M1 microglias
d_M2 = 0.015; % Death rate of M2 microglias
d_M1hat = 0.015; % Death rate of M1 macrophages
d_M2hat =0.015; % Death rate of M2 macrophages
d_Ao = 0.951; % Degradation rate of Aβ O
d_H = 58.71; % Degradation rate of HMGB-1
d_Talpha = 55.45; % Degradation rate of TNF-α
d_tbeta = 3.33e2; % Degradation rate of TGF-β
d_I10 = 16.64; % Degradation rate of IL-10
d_p = 1.73;  % Degradation rate of MCP-1
R0 = 6; % Initial inflammation by ROS
M0 = 5e-2; % Monocytes concentration in blood
N0 = 0.14; % Reference density of neuron
Mg0 = 0.047; % Source of microglia
A0 = 0.14; %reference density of astrocytes
lK_AmyloidbetaO = 7e-3; % Michaelis-Mention coefficient for Ao
lk_Nd = 10e-3; % Michaelis-Mention coefficient for Nd
K_I10 = 2.5e-6; % Half-saturation of IL-10
K_Tbeta = 2.5e-7;% Half-saturation of TGF-β
K_M = 0.047; % Half-saturation of microglias
K_Mhat = 0.047; % Half-saturation of macrophages
K_M1 = 0.03; % Half-saturation of M1 microglias
K_M2 = 0.017; % Half-saturation of M2 microglias
K_M1hat = 0.04; % Half-saturation of M1hat macrophages
K_M2hat = 0.007; % Half-saturation of M2hat macrophages
K_Fi = 3.36e-10; % Half-saturation of intracellular NFTs
K_Fo = 2.58e-11;% Average of extracellular NFTs
K_AO = 1e-7; % Average of of Aβ O
K_P = 6e-9;% Half-saturation of MCP-1
K_Talpha = 4e-5; % Half-saturation of TNF-α 
d_AbetaO_M = 2e-3; % Clearance rate of Aoβ by microglia
lambda_beta_i = 9.51e-6; % production rate of Amyloid beta i
R = 6;
% Store all variable names in one place so that we can use for cycle for plotting etc 
Variable_name = {'Abeta_i','Abeta_o','tau','Fi','Fo','N','A','Nd','Ao','M1','M2','M1hat','M2hat','H','Tbeta','Talpha','I10','P'}; 
% Set the initial conditions
Init(1) = 10e-6; % initial condition for Abeta_i (u1)
Init(2) = 10e-8; % initial condition for Abeta_o (u2)
Init(3) = 1.37e-10; % initial condition for tau (u3)
Init(4) = 3.36e-10; % initial condition for Fi (u4)
Init(5) = 3.36e-11; % initial condition for Fo (u5)
Init(6) = 0.14; % initial condition for N (u6)
Init(7) = 0.14; % initial condition for A (u7)
Init(8) = 0; % initial condition for Nd (u8)
Init(9) = 0.02; % initial condition for Ao (u9)
Init(10) = 0.02; % initial condition for M1 (u10)
Init(11) = 0; % initial condition for M2 (u11)
Init(12) = 0; % initial condition for M1hat (u12)
Init(13) = 0; % initial condition for M2hat (u13)
Init(14) = 1.3e-11; % initial condition for H (u14)
Init(15) = 10e-6; % initial condition for Tbeta (u15)
Init(16) = 2e-5; % initial condition for Talpha (u16)
Init(17) = 10e-5; % initial condition for I10 (u17)
Init(18) = 5e-9; % initial condition for P (u18)
% Transform Init into a column vector
Init = Init';
% Set the division of space into pieces
x = linspace(0,1,20); % try something smaller than 101 if it takes too long to run;
% Set the division of time into pieces
t = linspace(0,10,40); % try to use different size than for x so that opne can track if needed
% Set the technical auxiliary variable m t0 0 based on the matlab convetion
% for the code of the equation
m = 0;
% Solve the pde (the functions  pdefun, pdeic and pdebc are defined lower
% in the code)
opts = odeset('RelTol',1e-2,'AbsTol',1e-2);
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t,opts);
% Plot the figure 2
for aux = 1:length(Init) % For all variables
    % Unpack the solution corresponding to that variable
    u = sol(:,:,aux);
    % Take the average of the solution with respect to location
    Avg_u = mean(u,2); 
    % Plot the average of u against time
    figure
    plot(t, Avg_u)
    xlabel('Time t')
    ylabel(['Average of ' Variable_name{aux}])
end
% Define auxiliary functions needed for solving the equations
    function [c,f,s] = pdefun(x,t,u,dudx)
        % c will be all 1s as all equations have (partial u)/(partial t) and no other crossterms
        c = ones(length(Init),1);  % vector of all 1s of the same length as u; 
        % f corresponds to the flux
        D_Ao = 4.32e-2; % diffusion coefficient of AbetaO
        D_H = 8.11e-2; % Diffusion coefficient of HMGB-1
        D_Talpha = 6.55e-2; % Diffusion coefficient of TNF-alpha
        D_Tbeta = 6.55e-2; % Diffusion coefficient of TGF-beta
        D_I10 = 6.04e-2; %Diffusion coefficient of IL-10
        D_P = 1.2e-1; %Diffusion coefficient of MCP-1
        f = [0; 0; 0; 0; 0; 0; 0; 0; D_Ao; 0; 0; 0; 0; D_H; D_Tbeta; D_Talpha; D_I10; D_P] .* dudx;
        s = zeros(1,18);
        
        s(1) = lambda_beta_i.*(1+R) - d_Abeta_i.*u(1);
        s(2) = u(1).*abs(1.8e-5) + (lambda_N .* (u(6)./N0)) + (lambda_A .* (u(7)./A0)) - ...
            (d_AbetaO_Mhat.*(u(12)+theta.*u(13)) + d_AbetaO_M.*(u(10)+theta.*u(11))).*(u(2)./u(2)+lK_AmyloidbetaO); % abs(1.8e-5) = abs (pn/pt)
        s(3) = (lambda_tau0 + lambda_tau.*R - d_tau .* u(3)) .* (u(6)./N0);
        s(4) = (lambda_F.*u(3) - d_Fi.*u(4)) .* (u(6)./N0);
        s(5) = u(4).*abs(1.8e-5) - d_Fo.*u(5);
        s(6) = -d_NF .* (u(4))./(u(4) + K_Fi).* u(6) - d_NT.*(u(16)./(u(16) + K_Talpha)).*((1 + ((u(17)./ K_I10))^-1) * u(6));
        s(7) = lambda_A_AmyloidbetaO.*u(2) + lambda_AT_alpha.*u(16) - d_A .* u(7);
        s(8) = d_NF .* (u(4))./(u(4) + K_Fi).* u(6) - d_NT.*(u(16)./(u(16) + K_Talpha)).*((1 + ((u(17)./ K_I10))^-1) * u(6)) - d_Nd_M.*(u(10)+u(11)).*(u(8)./u(8)+lk_Nd) - d_Nd_Mhat.*(u(12)+u(13)).*(u(8)./u(8)+lk_Nd);
        s(9) = lambda_AO.*u(2)-d_Ao.*u(9);
        s(10) = K_M .* ((lambda_MF .* (u(5)./(u(5)+K_Fo))) + lambda_MA.*(u(9)./(u(9)+K_AO))).* ((beta.*(u(16)./u(16)+K_Talpha)).*((beta.*(u(16)./u(16)+K_Talpha))+ (u(17)./u(17)+K_I10)))^-1 - (lambda_M1_Tbeta).* (u(15)./u(15)+K_Tbeta).*u(10) - d_M1.*u(10); % Number of Peripheral Anti-Inflammatory Microglias over Time
        s(11) = K_M .* ((lambda_MF .* (u(5)./(u(5)+K_Fo))) + lambda_MA.*(u(9)./(u(9)+K_AO))).*(u(17)./u(17)+K_I10).*((beta.*(u(16)./u(16)+K_Talpha))+(u(17)./u(17)+K_I10)).^-1 + (lambda_M1_Tbeta).* (u(15)./u(15)+K_Tbeta).*u(10) - d_M2.*u(11); % Number of Peripheral Anti-Inflammatory Microglias over Time
        s(12) = alpha.*(u(18)./u(18)+K_P).* (M0 - u(12)+u(13)).* ((beta.*(u(16)./u(16)+K_Talpha)).*(((beta.*(u(16)./u(16)+K_Talpha))+ (u(17)./u(17)+K_I10)))^-1) - lambda_M1hat_beta .* (u(15)./u(15)+K_Tbeta) .* u(12) - d_M1hat.*u(12);
        s(13) = alpha.*(u(18)./u(18)+K_P).* (M0 - u(12)+u(13)).* ((u(17)./u(17)+K_I10).*((beta.*(u(16)./u(16)+K_Talpha))+(u(17)./u(17)+K_I10)).^-1) + lambda_M1hat_beta.*(u(15)./u(15)+K_Tbeta).*u(12)-d_M2hat.*u(13);
        s(14) = lambda_H.*u(8) - d_H .* u(14);
        s(15) = lambda_Tbeta_M.*u(11) + lambda_Tbeta_Mhat .* u(13) - d_tbeta .* u(15);
        s(16) = lambda_Talpha_M1.*u(10) + lambda_Talpha_M1hat.*u(12) - d_Talpha.*u(16);
        s(17) = lambda_I10_M2 .*u(11) + lambda_I10_M2hat.*u(13) - d_I10.*u(17);
        s(18) = lambda_PA.*u(7) + lambda_PM2.*u(11) - d_p.*u(18); 
        s = s';
    end
    function u0 = pdeic(x)
        % Codes the initial conditions in the location x
        u0 = [10e-6 10e-8 1.37e-10 3.36e-10 3.36e-11 0.14 0.14 0 0.02 0.02 0 0 0 1.3e-11 10e-6 2e-5 10e-5 5e-9]';
    end
    function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
        % Codes the boundary conditions
        % We, at least initially, want the boundary conditions to be the same
        % as initial conditions, i.e., u1 constantly Init1 and the derivative w.r.t. x being 0
        % Code the conditions for x=0 (left)
        pl = ul - Init;  % This will make u on the left boundary the same as Init values
        ql = zeros(length(Init),1);  % vector of all 0s of the same length as u,  This will mean that derivatives are not considered on the left boundary
        % Code the conditions for x=1 (right)
        pr = zeros(length(Init),1);  % vector of all 0s of the same length as u, This will mean that the u values are not considered on the right boundary
        qr = ones(length(Init),1);  % vector of all 1s of the same length as u, This means that the derivatives of u will be 0 on the right boundary
    end
end

[Bill Greene]

If you use the debugger (or add a print statement) in your pdefun, you will see that s(8) is NaN the first time it is called. I.e. either your initial values or the s(8) equation is invalid

Beyond that, here is a more general comment.

I have a lot of experience with pdepe but have never seen it used for a system as complicated as 18 equations. It is of course valid for such a system but the difficulties in obtaining a good solution to such a complicated problem are significant!

I suggest starting with a much simpler model of this behavior, say one that can be represented with two or three equations, and then adding additional equations after validating the results of the simpler models.

Answer 1

[Torsten]

"pdepe" is a solver for parabolic-elliptic partial differential equations. The positions with zero elements in the f-vector indicate ordinary differential equations. Thus "pdepe" is not suited to solve your system.

Sometimes "pdepe" succeeds in solving additional ordinary differential equations if you specify du/dx = 0 at both ends of the interval of integration. You could give it a try.

If this doesn't work, you should again check the source terms s and model parameters for correctness.

Are the diffusive terms for equations (10) - (13) neglectable ?