Using identical delay and length of integration domain for delay model (dde23)

Hi,

I've been trying to contsruct a biological process that is influenced by a drug every week, and causes few components of the immune system to react only 7-14 days later. Thus, the model consists a for loop with index of the number of weeks, and lags array for the delay. The results do differ from the model without delay, but if I change delay of 7 days to 14 days I get the same solution.

Is there any problem with my code? I just can't think of another option to model treatment without using for loop, but to me it seems that in this code there is something wrong with "retrieving" the history 7 or 14 days back. Thank you.

clc
a=6;
b=10;
c=9;
lags = [7 ];

y0 = [ 200 9.5*5*10^3 0.5*5*10^3 10^3 100 100];
tspan =[0,6/24];
n_weeks =6;


opt = odeset('AbsTol',1e-9,'RelTol',1e-6);


for i=1:3
  
sol =dde23(@(t,y,Z)ddefun(t,y,Z,a,b,c),lags,y0,tspan,opt);
semilogy(sol.x,sol.y(1,:),'r',sol.x,sol.y(2,:),'k',sol.x,sol.y(3,:),'b',sol.x,sol.y(4,:),'m',sol.x,sol.y(5,:),'g',sol.x,sol.y(6,:),'y')
legend('1','2','3','4','5','6')

xlabel('Time t');
ylabel('Solution y');
hold on

y0=sol.y(:,end);
tspan = [7*(i-1)+6/24,7*(i-1)+22/24];
m1=0;
sol =dde23(@(t,y,Z)ddefun(t,y,Z,a,b,c),lags,y0,tspan,opt);
semilogy(sol.x,sol.y(1,:),'r',sol.x,sol.y(2,:),'k',sol.x,sol.y(3,:),'b',sol.x,sol.y(4,:),'m',sol.x,sol.y(5,:),'g',sol.x,sol.y(6,:),'y')
legend('1','2','3','4','5','6')
xlabel('Time t');
ylabel('Solution y');
hold on
y0=sol.y(:,end);
mu1=0;
m1=0;
tspan = [7*(i-1)+22/24,7*i];
sol =dde23(@(t,y,Z)ddefun(t,y,Z,a,b,c),lags,y0,tspan,opt);
semilogy(sol.x,sol.y(1,:),'r',sol.x,sol.y(2,:),'k',sol.x,sol.y(3,:),'b',sol.x,sol.y(4,:),'m',sol.x,sol.y(5,:),'g',sol.x,sol.y(6,:),'y')
legend('1','2','3','4','5','6')
xlabel('Time t');
ylabel('Solution y');
hold on
y0=sol.y(:,end);
tspan = [7*i,7*i+6/24];

end


hold off
%-------------------------------------------
function dydt = ddefun(t,y,Z,a,b,c) %equation being solved
 ylag1 = Z(:,1);

 dydt = [ a*y(1);
b*ylag1(2)*y(1);
c*y(1)+a*ylag1(3);%
y(4);
ylag1(5)*ylag1(6);
y(6)];%

end

Answer 1

[Karan Singh]

Hi Marom,

The main issue with your code is that it does not properly account for the delay when updating "y0" and "tspan" within the loop. The delay in the biological response is not being modeled as intended, and this could be why changing the delay from 7 to 14 days does not affect the results.

When using "dde23" for delay differential equations, it's important to ensure that the history function (or the initial state y0) accounts for the delay period. However, in your loop, you're updating y0 with the end state of the previous iteration without considering the delay.

Here's a revised version of your loop that attempts to correct this by using a history function that can handle the delay:

clc
a = 6;
b = 10;
c = 9;
lags = [7]; % This should be in units consistent with tspan (e.g., days)
y0 = [200, 9.5*5*10^3, 0.5*5*10^3, 10^3, 100, 100];
tspan = [0, 6/24];
n_weeks = 6;
opt = odeset('AbsTol', 1e-9, 'RelTol', 1e-6);

% History function for the initial condition
history = @(t) y0;

for i = 1:n_weeks
    % Define the time span for this iteration
    tstart = 7 * (i - 1);
    tend = 7 * i;
    tspan = [tstart, tend];
    
    % Solve the DDE
    sol = dde23(@(t,y,Z) ddefun(t,y,Z,a,b,c), lags, history, tspan, opt);
    
    % Plot the results
    semilogy(sol.x, sol.y(1,:), 'r', sol.x, sol.y(2,:), 'k', sol.x, sol.y(3,:), 'b', sol.x, sol.y(4,:), 'm', sol.x, sol.y(5,:), 'g', sol.x, sol.y(6,:), 'y')
    legend('1', '2', '3', '4', '5', '6')
    xlabel('Time t');
    ylabel('Solution y');
    hold on
    
    % Update the history function to include the latest solution
    % This creates an anonymous function that interpolates the solution
    % for times within the last week (or whatever the delay period is)
    history = @(t) deval(sol, t);
end

hold off

%-------------------------------------------
function dydt = ddefun(t, y, Z, a, b, c) % equation being solved
    ylag1 = Z(:,1);
    dydt = [a*y(1);
            b*ylag1(2)*y(1);
            c*y(1)+a*ylag1(3);
            y(4);
            ylag1(5)*ylag1(6);
            y(6)];
end

Hope this resolves your doubt!