How to write a step function, x(t-tau), within a differential equation

I am trying to solve and plot the following system of equations, however I am struggling with plotting the T(t-tau) and E(t-tau) functions in Matlab

I have been trying to use the heaviside function but the resulting graphs are not changing with variance in tau. Any help would be greatly appreciated, thank you!

function final_project_4()

r1 = 0.1; k1 = 0.0001; r2 = 0.03; alpha = 24; k2 = 0.00005; k3 = 0.000024;
u2 = 0.3; gamma = 0.3; c = 3.8; u1 = 0.03; epislon = 0.1; N = 275;

tau = 15;

TSpan = [0 250];
Y0 = [50 780];
options = odeset('RelTol', 1e-10,'AbsTol',[1e-10,1e-10]);

[t,y] = ode45(@model,TSpan,Y0,options);

plot(t,y)
title('tau = 15')
xlabel('time (days)')
ylabel('concentration')
legend('T(t)-cancer cells','E(t)-effector cells')

    function res = model(t,y)
        T = y(1);
        E = y(2);

        Ttau = T*heaviside(t-tau) ;
        Etau = E*heaviside(t-tau) ;

        dTdt = r1*T - k1*T*E;
        dEdt = r2*T + alpha - u1*E - k1*T*E + (1-epislon)*k1*Ttau*Etau;
        res = [dTdt;dEdt];
    end
end

Answer 1

[Nivedita]

Hello Jack,

The issue you're encountering with the heaviside function is that it doesn't actually delay the functions "T(t)" and "E(t)" by "tau"; it merely multiplies the current values of "T" and "E" by 1 for (t > tau) (and 0 otherwise), which doesn't achieve the intended delay effect. To properly incorporate the delay, you need to store or calculate the values of "T(t-tau)" and "E(t-tau)" at the delayed times. Unfortunately, MATLAB's standard "ode45" solver does not directly support delay differential equations (DDEs).

Here is how you can try to approach solving your system using "dde23":

function final_project_4_dde()

% Parameters
r1 = 0.1; k1 = 0.0001; r2 = 0.03; alpha = 24; 
u1 = 0.03; epsilon = 0.1; 
tau = 15;

% Define the delays
lags = [tau];

% Initial conditions
Y0 = [50; 780];

% History function (for t <= 0)
history = @(t) Y0;

% Time span
tspan = [0 250];

% Solve the DDE
sol = dde23(@(t,y,Z) model(t,y,Z), lags, history, tspan);

% Plot the results
figure;
plot(sol.x, sol.y);
title(['System Dynamics with \tau = ', num2str(tau)]);
xlabel('Time (days)');
ylabel('Concentration');
legend('T(t) - cancer cells', 'E(t) - effector cells');

    function dydt = model(t,y,Z)
        T = y(1);
        E = y(2);
        
        % Delayed values
        Ttau = Z(1,1);
        Etau = Z(2,1);
        
        % System of equations
        dTdt = r1*T - k1*T*E;
        dEdt = r2*T + alpha - u1*E - k1*T*E + (1-epsilon)*k1*Ttau*Etau;
        
        dydt = [dTdt; dEdt];
    end
end

The "model" function now takes an additional argument "Z", which contains the delayed states of "T" and "E". The "lags" variable specifies the delay "tau" for both variables. The "history" function specifies the initial state of the system for (t<=0), which is necessary for "dde23" to compute the values at (t - tau).

You will have to adjust the "history" function as per your specific needs, especially if the initial conditions for "T" and "E" should be different for (t < 0).

For more information on "dde23", refer the following documentation link:

https://www.mathworks.com/help/matlab/ref/dde23.html