
Problem solving a biological model by ode45
3 views (last 30 days)
Show older comments
Marom Yossef
on 6 Feb 2022
Commented: Davide Masiello
on 6 Feb 2022
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:
- The system is stiff. However, i'm getting similar results using ode15s solver.
- 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).
- 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
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
0 Comments
Accepted Answer
Davide Masiello
on 6 Feb 2022
Edited: Davide Masiello
on 6 Feb 2022
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.
4 Comments
Star Strider
on 6 Feb 2022
Definitely impressive!
+1
Out of curiosity, what (substance or drug) is being modeled here, and what are the variables (or compartments)?
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!