Problem solving a biological model by ode45

I am taking part in a biomathematics project, and really new to solving ODEs in MATLAB.

trying to understand the way to solve and analyze ode, but I keep getting a solution that is derived by ignorig negative data (which is unrealitic in biology).

I guess my coding approach is wrong, so I thought that maybe it to deal with one of the following:

1. The system is stiff. However, i'm getting similar results using ode15s solver.
2. Breaking time-span into intervals with repeated inital data for one variable (the first equation is required to model drug injection each 7 days by Delta-Dirac function).
3. The model is supposed to consider delay in time for few variables but I couldn't find any coding approach to make these translation.

I would be grateful for any advice.

mu1=1;
mu2=0.41;
mu=mu2/mu1;
p1=1.25/mu1;
p2=0.285;
p3=1.1;
p4=0;
p5=0.003;
alpha=0.52;
beta= 0.011;
b=5;
r=0.032;
tspan =linspace(0,7);
tspan2 =linspace(7,14);
tspan3 =linspace(14,21);
tspan4 =linspace(21,28);

y0 = [5 mu1/p2 mu1/p2 mu1/p2 ];


[t,y] = ode45(@(t,y) odefcn (t,y,mu,r,p1,p2,p3,p4,p5,alpha,beta), tspan, y0);
semilogy(t,y(:,1),'b',t,y(:,2),'k',t,y(:,3),'r',t,y(:,4),'g')
legend('B','E','Ti','Tu')

hold on

l0=y(end,:);
l0(1)=5;

[t,y] = ode45(@(t,y) odefcn (t,y,r,mu,p1,p2,p3,p4,p5,alpha,beta),  tspan2, l0);
semilogy(t,y(:,1),'b',t,y(:,2),'k',t,y(:,3),'r',t,y(:,4),'g')
legend('B','E','Ti','Tu')

hold off

Warning: Negative data ignored

function dydt = odefcn(t,y,mu,r,p1,p2,p3,p4,p5,alpha,beta)
  dydt = zeros(4,1);
  B=y(1);
  E=y(2);
  Ti=y(3);
  Tu=y(4);
 
  dydt(1) = B*(-1-p1*E-p2*Tu);
  dydt(2) =E*(-mu+p4*B-p5*Ti)+alpha*Ti;
  dydt(3) =-p3*E*Ti+p2*B*Tu;
  dydt(4) =Tu*(-p2*B+r*(1-beta*Tu));

end

Answer 1

[Davide Masiello]

Note that the solution is not ignoring negative data, but you are plotting it using a logarithmic y-axis, reason for which it cannot plot negative values.

In order to account for the drug injection, you could do something like this

clear,clc,clf

mu1=1;
mu2=0.41;
mu=mu2/mu1;
p1=1.25/mu1;
p2=0.285;
p3=1.1;
p4=0;
p5=0.003;
alpha=0.52;
beta= 0.011;
b=5;
r=0.032;
tspan = [0,7];
n_weeks = 4;
y0 = [5 mu1/p2 mu1/p2 mu1/p2 ];
opt = odeset('AbsTol',1e-9,'RelTol',1e-6);

for i = 1:n_weeks
sol(i)  = ode15s(@(t,y) odefcn (t,y,mu,r,p1,p2,p3,p4,p5,alpha,beta), tspan, y0, opt);
tspan   = [tspan(2),tspan(2)+7];
y0      = sol(i).y(:,end); y0(1) = 5;
end
plot([sol(:).x],[sol(:).y])
legend('B','E','Ti','Tu')
xlabel('Time, s')

function dydt = odefcn(t,y,mu,r,p1,p2,p3,p4,p5,alpha,beta)
  dydt = zeros(4,1);
  B=y(1);
  E=y(2);
  Ti=y(3);
  Tu=y(4);
 
  dydt(1) = B*(-1-p1*E-p2*Tu);
  dydt(2) =E*(-mu+p4*B-p5*Ti)+alpha*Ti;
  dydt(3) =-p3*E*Ti+p2*B*Tu;
  dydt(4) =Tu*(-p2*B+r*(1-beta*Tu));
end

Regarding the delay, I wouldn't know what to say without more info, but I would recommend checking the following documentation.

https://it.mathworks.com/help/matlab/delay-differential-equations.html

[Marom Yossef]

It is very helpful, thank you!

I just don't understand how to code an impluse or something like delta-Dirac function each 7 days.

[Davide Masiello]

Apologies, I have edited the answer to consider the other issues.

[Star Strider]

Definitely impressive!

+1

Out of curiosity, what (substance or drug) is being modeled here, and what are the variables (or compartments)?

[Davide Masiello]

Thanks :)