I am programing a mathematical model from an article about angiogensis but when I try to plot the functions with the given parameteres and initial values, the program doesn't run for the entirety of the timespan I defined, showing a part of the graph but ending abrupty and showing the message:
Warning: Failure at t=5.404287e+01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed
(1.136868e-13) at time t.
> In ode45 (line 360)
The link for the full article with all the functions and the graphics I am trying to copy is https://www.sciencedirect.com/science/article/pii/S0022519316001429
I have tried other ode beside ode45, but this one was the only one that showed some promise...
The code I have currently is this:
Function for the first imagine of Fig. 3
function dy = derivatives1(t,i)
%% Default parameter set (Table C1)
V0 = 8E-3; % V0 = Reference VEGF concentration (uM)
A0 = 8E-3; % A0 = Reference Ang1/2 concentration (uM)
p0 = 8E-3; % p0 = Reference PDGF concentration (uM)
E0 = 8E-8; % E0 = Reference EC density (uM)
P0 = 4E-8; % P0 = Reference PC density (uM)
phi = 10^5; % phi = Average receptors per cell
v_kon = 1.88E4; % v_kon = VEGF/VEDGFR2 association (uM^-1 h^-1)
v_koff = 0.99; % v_koff = VEGF/VEDGFR2 dissociation (h^-1)
v_kint = 0.86; % v_kint = VEGF internalisation (h^-1)
alphav = 3E-3; % alphav = VEGF delivery (uM^-1 h^-1)
deltav = 1; % deltav = VEGF decay (h^-1)
a1_kon = 417; % a1_kon = Ang1/Tie2 association (uM^-1 h^-1)
a1_koff = 1.25; % a1_koff = Ang1/Tie2 dissociation (h^-1)
a1_kint = 1; % a1_kint = Ang1 internalisation (h^-1)
a2_kon = 417; % a2_kon = Ang2/Tie2 association (uM^-1 h^-1)
a2_koff = 1.25; % a2_koff = Ang2/Tie2 dissociation (h^-1)
deltaa1 = 4.69; % deltaa1 = Ang1 decay (h^-1)
deltaa2 = 0.0625; % deltaa2 = Ang2 decay (h^-1)
deltaaw = 0.031; % deltaaw = Intracellular Ang2 decay (h^-1)
alphaa1 = 1E6; % alphaa1 = Ang1 production (h^-1)
alphaaw = 0.8; % alphaaw = Ang2 production (h^-1)
alphaa2 = 3; % alphaa2 = Ang2 release (h^-1)
alphap = 6.7E3; % alphap = PDGF production (h^-1)
p_kon = 1.88E4; % p_kon = PDGF/PDGFR-Beta association (uM^-1 h^-1)
p_koff = 0.86; % p_koff = PDGF/PDGFR-Beta dissociation (h^-1)
deltap = 1.2; % deltap = PDGF decay (h^-1)
theta = 0.0417; % theta = EC proliferation (h^-1)
rho = 0.02; % rho = EC death (h^-1)
kappa = 0.08; % kappa = PC proliferation (h^-1)
epsilon = 0.04; % epsilon = PC death (h^-1)
zeta = 5E5; % zeta = PC attachment (uM^-1 h^-1)
nu = 5E-3; % nu = PC detachment (h^-1)
eta = 0.01; % eta = Maturity increase (h^-1)
mu = 0.01; % mu = Maturity decrease (h^-1)
cep = 0.3; % cep = Threshold bound VEGF for EC proliferation
cd = 0.1; % cd = Threshold bound VEGF for EC death
ca = 0.3; % ca = Threshold bound Ang2 for EC proliferation
cpp = 0.02; % cpp = Threshold bound PDGF for PC proliferation
lambda = 8; % lambda = Ang2/Ang1 ratio for maturity decrease
%% Variable groupings
% Et = Total endothelial cells (uM)
Et = i(13) + i(14) + i(16);
% Pt = Total pericyte cells (uM)
Pt = i(15) + i(16);
% vb = Bound fraction of PDGFR-Beta receptors (None)
vb = i(3)/(i(2)+i(3));
% ab = Fraction of Tie2 receptors bound by Ang2 (None)
ab = i(9)/(i(7)+i(8)+i(9));
% pb = Bound fraction of PDGFR-Beta receptors (None)
pb = i(12)/(i(11)+i(12));
%% Model
dy = zeros(16,1);
% Biochemical Sub-Model
% Free VEGF
dy(1) = alphav/(((i(14)+i(16))/E0)+1) - v_kon*i(1)*i(2) + v_koff*i(3) - deltav*i(1);
% VEGFR2
dy(2) = - v_koff*i(1)*i(2) + (v_koff+v_kint)*i(3) + phi*theta*(vb^2/(vb^2+cep^2))*(1+(ab^2/(ab^2+ca^2)))*(1-(Et/E0))*i(13) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(2);
% VEGF/VEGFR2
dy(3) = v_kon*i(1)*i(2) - (v_koff+v_kint)*i(3) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(3);
% Free Ang1
dy(4) = alphaa1*Pt - a1_kon*i(4)*i(7) + a1_koff*i(8) - deltaa1*i(4);
% Intracellular Ang2
dy(5) = alphaaw*(1+vb)*((A0/E0)*i(13)-i(5)) - deltaaw*i(5) - alphaa2*i(5)*((i(5))^2/((i(5))^2+((i(13)*A0)/(2*E0))^2)) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(5);
% Free Ang2
dy(6) = alphaa2*i(5)*((i(5))^2/((i(5))^2+((i(13)*A0)/(E0*2))^2)) - a2_kon*i(6)*i(7) + a2_koff*i(9) - deltaa2*i(6);
% Tie2
dy(7) = - a1_kon*i(4)*i(7) + (a1_koff+a1_kint)*i(8) - a2_kon*i(6)*i(7) + a2_koff*i(9) + phi*theta*(vb^2/(vb^2+cep^2))*(1+(ab^2/(ab^2+ca^2)))*(1-(Et/E0))*i(13) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(7);
% Ang1/Tie2
dy(8) = a1_kon*i(4)*i(7) - (a1_koff+a1_kint)*i(8) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(8);
% Ang2/Tie2
dy(9) = a2_kon*i(6)*i(7) - a2_koff*i(9) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(9);
% Free PDGF
dy(10) = alphap*i(13) - p_kon*i(10)*i(11) + p_koff*i(12) - deltap*i(10);
% PDGFR-Beta
dy(11) = - p_kon*i(10)*i(11) + p_koff*i(12) + phi*kappa*(pb^2/(pb^2+cpp^2))*(1-(Pt/P0))*i(15) - epsilon*((i(6)/A0)+(cpp^2/(cpp^2+pb^2))*i(11));
% PDGF/PDGFR-Beta
dy(12) = p_kon*i(10)*i(11) - p_koff*i(12) - epsilon*((i(6)/A0)+(cpp^2/(cpp^2+pb^2)))*i(12);
% Cell Dynamics Sub-model
% Immature ECs
dy(13) = theta*(vb^2/(vb^2+cep^2))*(1+(ab^2/ab^2+ca^2))*(1-(Et/E0))*i(13) - rho*(cd^2/cd^2+vb^2)*i(13) - zeta*i(15)*i(13) + nu*i(16) + epsilon*((i(6)/A0)+(cpp^2/(cpp^2+pb^2)))*i(16) - eta*(lambda^2/(lambda^2+(i(9)/i(8))^2))*i(13) + mu*((i(9)/i(8))^2/((i(9)/i(8))^2+lambda^2))*i(14);
% Mature ECs
dy(14) = - zeta*i(15)*i(14) + eta*(lambda^2/(lambda^2+(i(9)/i(8))^2))*i(13) + mu*((i(9)/i(8))^2/((i(9)/i(8))^2+lambda^2))*i(14);
% Free Pericytes
dy(15) = kappa*(pb^2/(pb^2+cpp^2))*(1-(Pt/P0))*i(15) - epsilon*((i(6)/A0)+(cpp^2/(cpp^2+pb^2)))*i(15) - zeta*i(15)*(i(13)+i(14)) + nu*i(16);
% EC/PC Complex
dy(16) = zeta*i(15)*(i(13)+i(14)) - nu*i(16) - epsilon*((i(6)/A0)+(cpp^2/(cpp^2+pb^2)))*i(15);
end
Function for the second imagine of Fig. 3
function dy = derivatives2(t,i)
%% Default parameter set (Table C1)
V0 = 8E-3; % V0 = Reference VEGF concentration (uM)
A0 = 8E-3; % A0 = Reference Ang1/2 concentration (uM)
p0 = 8E-3; % p0 = Reference PDGF concentration (uM)
E0 = 8E-8; % V0 = Reference EC density (uM)
P0 = 4E-8; % P0 = Reference PC density (uM)
phi = 10^5; % phi = Average receptors per cell
v_kon = 1.88E4; % v_kon = VEGF/VEDGFR2 association (uM^-1 h^-1)
v_koff = 0.99; % v_koff = VEGF/VEDGFR2 dissociation (h^-1)
v_kint = 0.86; % v_kint = VEGF internalisation (h^-1)
alphav = 3E-3; % alphav = VEGF delivery (uM^-1 h^-1)
deltav = 1; % deltav = VEGF decay (h^-1)
a1_kon = 417; % a1_kon = Ang1/Tie2 association (uM^-1 h^-1)
a1_koff = 1.25; % a1_koff = Ang1/Tie2 dissociation (h^-1)
a1_kint = 1; % a1_kint = Ang1 internalisation (h^-1)
a2_kon = 417; % a2_kon = Ang2/Tie2 association (uM^-1 h^-1)
a2_koff = 1.25; % a2_koff = Ang2/Tie2 dissociation (h^-1)
deltaa1 = 4.69; % deltaa1 = Ang1 decay (h^-1)
deltaa2 = 0.0625; % deltaa2 = Ang2 decay (h^-1)
deltaaw = 0.031; % deltaaw = Intracellular Ang2 decay (h^-1)
alphaa1 = 1E6; % alphaa1 = Ang1 production (h^-1)
alphaaw = 0.8; % alphaaw = Ang2 production (h^-1)
alphaa2 = 3; % alphaa2 = Ang2 release (h^-1)
alphap = 6.7E3; % alphap = PDGF production (h^-1)
p_kon = 1.88E4; % p_kon = PDGF/PDGFR-Beta association (uM^-1 h^-1)
p_koff = 0.86; % p_koff = PDGF/PDGFR-Beta dissociation (h^-1)
deltap = 1.2; % deltap = PDGF decay (h^-1)
theta = 0.0417; % theta = EC proliferation (h^-1)
rho = 0.02; % rho = EC death (h^-1)
kappa = 0.08; % kappa = PC proliferation (h^-1)
epsilon = 0.04; % epsilon = PC death (h^-1)
zeta = 0.24; % zeta = PC attachment (uM^-1 h^-1)
nu = 5E-3; % nu = PC detachment (h^-1)
eta = 0.01; % eta = Maturity increase (h^-1)
mu = 0.01; % mu = Maturity decrease (h^-1)
cep = 0.3; % cep = Threshold bound VEGF for EC proliferation
cd = 0.1; % cd = Threshold bound VEGF for EC death
ca = 0.3; % ca = Threshold bound Ang2 for EC proliferation
cpp = 0.02; % cpp = Threshold bound PDGF for PC proliferation
lambda = 8; % lambda = Ang2/Ang1 ratio for maturity decrease
%% Variable groupings
% Et = Total endothelial cells (uM)
Et = i(13) + i(14) + i(16);
% Pt = Total pericyte cells (uM)
Pt = i(15) + i(16);
% vb = Bound fraction of PDGFR-Beta receptors (None)
vb = i(3)/(i(2)+i(3));
% ab = Fraction of Tie2 receptors bound by Ang2 (None)
ab = i(9)/(i(7)+i(8)+i(9));
% pb = Bound fraction of PDGFR-Beta receptors (None)
pb = i(12)/(i(11)+i(12));
%% Model
dy = zeros(16,1);
% Biochemical Sub-Model
% Free VEGF
dy(1) = alphav/(((i(14)+i(16))/E0)+1) - v_kon*i(1)*i(2) + v_koff*i(3) - deltav*i(1);
% VEGFR2
dy(2) = - v_koff*i(1)*i(2) + (v_koff+v_kint)*i(3) + phi*theta*(vb^2/(vb^2+cep^2))*(1+(ab^2/(ab^2+ca^2)))*(1-(Et/E0))*i(13) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(2);
% VEGF/VEGFR2
dy(3) = v_kon*i(1)*i(2) - (v_koff+v_kint)*i(3) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(3);
% Free Ang1
dy(4) = alphaa1*Pt - a1_kon*i(4)*i(7) + a1_koff*i(8) - deltaa1*i(4);
% Intracellular Ang2
dy(5) = alphaaw*(1+vb)*((A0/E0)*i(13)-i(5)) - deltaaw*i(5) - alphaa2*i(5)*((i(5))^2/((i(5))^2+((i(13)*A0)/(2*E0))^2)) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(5);
% Free Ang2
dy(6) = alphaa2*i(5)*((i(5))^2/((i(5))^2+((i(13)*A0)/(E0*2))^2)) - a2_kon*i(6)*i(7) + a2_koff*i(9) - deltaa2*i(6);
% Tie2
dy(7) = - a1_kon*i(4)*i(7) + (a1_koff+a1_kint)*i(8) - a2_kon*i(6)*i(7) + a2_koff*i(9) + phi*theta*(vb^2/(vb^2+cep^2))*(1+(ab^2/(ab^2+ca^2)))*(1-(Et/E0))*i(13) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(7);
% Ang1/Tie2
dy(8) = a1_kon*i(4)*i(7) - (a1_koff+a1_kint)*i(8) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(8);
% Ang2/Tie2
dy(9) = a2_kon*i(6)*i(7) - a2_koff*i(9) - rho*(cd^2/(cd^2+vb^2))*(i(13)/Et)*i(9);
% Free PDGF
dy(10) = alphap*i(13) - p_kon*i(10)*i(11) + p_koff*i(12) - deltap*i(10);
% PDGFR-Beta
dy(11) = - p_kon*i(10)*i(11) + p_koff*i(12) + phi*kappa*(pb^2/(pb^2+cpp^2))*(1-(Pt/P0))*i(15) - epsilon*((i(6)/A0)+(cpp^2/(cpp^2+pb^2))*i(11));
% PDGF/PDGFR-Beta
dy(12) = p_kon*i(10)*i(11) - p_koff*i(12) - epsilon*((i(6)/A0)+(cpp^2/(cpp^2+pb^2)))*i(12);
% Cell Dynamics Sub-model
% Immature ECs
dy(13) = theta*(vb^2/(vb^2+cep^2))*(1+(ab^2/ab^2+ca^2))*(1-(Et/E0))*i(13) - rho*(cd^2/cd^2+vb^2)*i(13) - zeta*i(15)*i(13) + nu*i(16) + epsilon*((i(6)/A0)+(cpp^2/(cpp^2+pb^2)))*i(16) - eta*(lambda^2/(lambda^2+(i(9)/i(8))^2))*i(13) + mu*((i(9)/i(8))^2/((i(9)/i(8))^2+lambda^2))*i(14);
% Mature ECs
dy(14) = - zeta*i(15)*i(14) + eta*(lambda^2/(lambda^2+(i(9)/i(8))^2))*i(13) + mu*((i(9)/i(8))^2/((i(9)/i(8))^2+lambda^2))*i(14);
% Free Pericytes
dy(15) = kappa*(pb^2/(pb^2+cpp^2))*(1-(Pt/P0))*i(15) - epsilon*((i(6)/A0)+(cpp^2/(cpp^2+pb^2)))*i(15) - zeta*i(15)*(i(13)+i(14)) + nu*i(16);
% EC/PC Complex
dy(16) = zeta*i(15)*(i(13)+i(14)) - nu*i(16) - epsilon*((i(6)/A0)+(cpp^2/(cpp^2+pb^2)))*i(15);
end
The actual plot of the graph and final equations needed:
%% Dependent variables and variable groupings (Table 1)
% Definition and default initial conditions used for numerical simulations
% of dependent variables and variable groupings
% Angiogenic factors (VEGF) bind to VEGFR-2 receptors
v = 8E-5; % v = Free VEGF concentration (uM)
rv = 6E-4; % rv = Total unbound VEGFR-2 receptors (uM)
bv = 6E-4; % bv = Total bound VEGFR-2 receptors (uM)
% Angiopoietin ligands (Ang1, Ang2) bind to Tie2 receptors
a1 = 8E-5; % a1 = Free Ang1 (uM)
aw = 8E-5; % aw = Total Ang2 stored in all WPBs (uM)
a2 = 8E-5; % a2 = Free Ang2 (uM)
ra = 6E-4; % ra = Total unbound Tie2 receptors (uM)
ba1 = 3E-4; % ba1 = Total bound Ang1/Tie2 (uM)
ba2 = 3E-4; % ba2 = Total bound Ang2/Tie2 (uM)
% PDGF-B binds to PDGFR-Beta receptors
p = 8E-5; % p = PDGF-B (uM)
rp = 4E-4; % rp = Total unbound PDGFR-Beta (uM)
bp = 4E-4; % rp = Total bound PDGFR-Beta (uM)
% Endothelial cells (EC) and pericyte cells (PC)
Ei = 4E-9; % Ei = Number density of immature ECs (uM)
Em = 8E-9; % Em = Number density of mature ECs (uM)
P = 8E-9; % P = Number density of free pericytes (uM)
Ep = 0; % Ep = Number density of ECs with pericyte coverage (uM)
%Initial condition vector
% 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
i = [v rv bv a1 aw a2 ra ba1 ba2 p rp bp Ei Em P Ep];
%% Integration
timespan = [0 200]; %200 days
[t1,y1] = ode45(@derivatives1,timespan,i);
[t2,y2] = ode45(@derivatives2,timespan,i);
%% Figure 3 plots
figure('Name','Numerical simulations of the full model')
subplot(1,2,1)
plot(t1, y1(:,13),t1,y1(:,14),t1,y1(:,15),t1,y1(:,16))
xlabel('Time (days)')
ylabel('Normalised Density')
subplot(1,2,2)
plot(t2,y2(:,13),t2,y2(:,14),t2,y2(:,15),t2,y2(:,16))
xlabel('Time (days)')
ylabel('Normalised Density')
