Modelling a system of differential equations with dirac delta function using Matlab

Question

I have a system similar to the model in here. This is a simplified version with only two states, bacteria ($P$) and antibiotic ($C$). After 3 days have passed (the system is run from day 0) antibiotic treatment is started. And three doses are give ($d_1,d_2$ and $d_3$). $d_1$ is given at time $t_1=3$, $d_2$ is given at time $t_2=5$ and $d_3$ is given at time $t_3=10$. $\delta$ is a Dirac delta function. Each dose is a constant amount ($10$).

The system of ODEs are:

${dP\over dt}=rP(1-{P \over k})-CP-\theta P$
${dC\over dt}=d_1\delta(t-t_1)+d_2\delta(t-t_2)+d_3\delta(t-t_3)-gC$.
In this system $r,k,\theta,g$ are all fixed parameter values.

Before the start of treatment (before 3rd day) the system will be just ${dP\over dt}=rP(1-{P \over k})-\theta P$ .

I am trying to implement this in Matlab and this is the code that I have written,

p=1000;%initial size of the unattached bacteria
c=0;%initial antibiotic concentration 
time=1000;%time in hours

options = odeset('AbsTol',1e-13,'RelTol',1e-11);

[t,y]=ode45(@(t,y)DrugModel(t,y),0:time,[p,c],options);

figure(1)
plot(t/24,y(:,2))
title('antibiotic')

function s= DrugModel(t,y)
r=2.7726;
k=1000;
theta=0.02;
g=0.48;  

s=zeros(2,1);

if t>= 3*24 %treatment is initiated on 3rd day

    s(1)=(r*y(1))*(1-((y(1))/k))-0.25*y(2)*y(1)-theta*y(1) ;%rate of change of bacteria

    if (t-3*24)<0.01 %first dose
        s(2)=10;
    end
    if (t-5*24)<0.01%second dose
        s(2)=10;
    end
    if (t-10*24)<0.01%third dose
        s(2)=10;
    end
    s(2)=-g*y(2);%antibiotic concentration decays.

else %before the start of treatment

    s(1)=(r*y(1))*(1-((y(1))/k))-theta*y(1);%rate of change of bacteria
    s(2)=0;
end

end

But, with the way I have implemented it the antibiotic concentration remains at 0 the whole time period. However, there should be pulses when doses are included and in other times it should decay exponentially.

Can someone please let em know how I can implement this system in Matlab.
(The correct system that I have is complex than this one so it can't be solved analytically)

Answer 1

I would code this up as an "impulsive differential equation" so loop over your ODE solver as follows (might not be perfect, but hopefully helpful) *I made k=100:

M-file for ODE:

function DM = DM(t,y)

r=2.7726;
k=100;
theta=0.02;
g=0.48; 

DM = zeros(2,1);    
DM(1) = r*y(1)*(1-(y(1)/k))-y(2)*y(1)-theta*y(1);
DM(2) = -g*y(2);

Now for the file that runs the model:

close all

TIME=[3 2 5 20]; % length of time between doses , then run an extra 20 "days"
Dose=10;
y0=[100 0];
Y=[];

for i=1:length(TIME)
    t=0:1:TIME(i) % run until next dose 
    [t,ytemp] = ode45(@DM,t,y0);
    Y=[Y;ytemp];

    yneg=ytemp(end,:); % collect ode solution
    yplus=zeros(1,2);

    yplus(1)=yneg(1); % bacteria is unchanged
    yplus(2)=yneg(2)+Dose; % drug added 
    y0=yplus; % set new initial conditions 
end


%plot it 
plot(Y(:,1),'g','LineWidth',3)
hold on
plot(Y(:,2),'r','LineWidth',3)
legend('Bacteria','Drug')
set(gca,'FontSize',18)
set(gcf,'Color','w')