Heaviside step function output as a parameter to another function matlab

Question

I am trying to implement the mathematical model described in this reference:

Chad Liu, Chuan-Yuan Li, Fan Yuan. Mathematical Modeling of the phoenix rising pathway. Plos One, Computational Biology.v.10, 2014.

I would like the output variable dPLAdt to be calculated numerically obeying the Heaviside function, which calculates the variable C3 and C7. As you can see from the code below, the Heaviside function will be activated after time 1440. The program is running but not producing results matching the attached picture from the paper![https://i.stack.imgur.com/37Vjo.jpg][1]. I would like to match the paper's result.

Here is my code:

1) The function

function dPLAdt = PLAprod(t,PLA)


dPLAdt = zeros(2,1);

k1 = 2.8*10^(-5);
k3 = 3.24*10^(-6);

kminus5 = 0.09/60;
k5 = 0.06/60;
k2short = 144/60;
k2 = 11;
k4short = 26/60;
k4 = 12;
k7 = 0.06;
p = [k1*heaviside(t-1440); k3*heaviside(t-1440)]; % CASP 3 e CASP7
C3 = p(1); 
C7 = p(2);



dPLAdt(1) = kminus5 - k5*(PLA(1)) - (k2short*(PLA(1))*(C3))/(k2 + PLA(1)) - (k4short*(PLA(1))*(C7))/(k4 + PLA(1));


dPLAdt(2) = (k2short*(PLA(1))*(C3))/(k2 + (PLA(1))) + (k4short*(PLA(1))*(C7))/(k4 + PLA(1)) - k7*(PLA(2));

end

2) Script to call the solver

[T,Y] = ode45(@PLAprod,[0,2880],[1500 500]);
plot(T,Y(:,1),'-r', T,Y(:,2),'-b')
legend('iPLA','aPLA')

Answer 1

The Heaviside function effectively makes this a (very) stiff system that ode45 isn't designed to handle. You could attempt to adjust integration tolerances or use ode15s or another stiff solver, but what you really appear to have is a system that is piecewise in time. You can accurately and efficiently integrate it using a general adaptive solver like ode45 if you break up time (tspan) and create two slightly different integration functions. Just replace the Heaviside function with a constant low or high value. Then simulate the first from 0 to 1440 and the second from 1440 to 2880 using as initial conditions the end state of the first. See my answer to this question for an example of this.

You can think of the Heaviside function as an if statement, and it's almost always a bad idea to put an if, or any other form of branching, inside of an integration function. ODE solvers assume a certain degree of smoothness in the functions they're integrating.