Hello everyone,
I am trying to reproduce the numerical results of a published tumour–immune ODE model that includes Dirac-type impulsive inputs. The system has four coupled variables [U,Y,Z,X][U, Y, Z, X][U,Y,Z,X]. Between injections, the dynamics are continuous, and at each injection time an instantaneous jump is applied to UUU:
Continuous part:
dU/dt = –kYU·U
dY/dt = kYU·U – kEY·Y
dZ/dt = kZY·Y – kEZ·Z
dX/dt = kX·X – kEX·X – η·Y·X·exp(–λ·Z)
At each dose time t=ti:
U(ti+)=U(ti−)+dosei×massi​
I have implemented this in MATLAB using two scripts:
A main script, where I loop over the dose times, call ode45 between impulses, and apply the jump in UUU. For example:
% Mouse 1 parameters and initial conditions
y_curr = [M1.U0; M1.Y0; M1.Z0; M1.X0];
t_curr = M1.t0;

for i = 1:length(tDose1)
ti = tDose1(i);
if ti > M1.tend, break; end

% integrate continuously until next dose
[t_seg, Y_seg] = ode45(@(t,y) fM1_fun(t,y,M1), [t_curr ti], y_curr, options);

% store results (code omitted here)
if isempty(t_all1)
t_all1 = t_seg;
Y_all1 = Y_seg;
else
t_all1 = [t_all1; t_seg(2:end,:)];
Y_all1 = [Y_all1; Y_seg(2:end,:)];
end

% state just before the impulse
y_curr = Y_seg(end,:)’;

% apply Dirac-type impulse to U
Mi = mass_fun1(ti); % interpolated mass at ti
ui = uDose1(i); % dose at ti (mg/kg)
y_curr(1) = y_curr(1) + ui * Mi; % jump in U

t_curr = ti;
end
2. ODE function files, one for each parameter set (Mouse 1 and Mouse 2). For example, for Mouse 1:
function dydt = fM1_fun(t,y,M1)

dU = -M1.kYU * y(1);
dY = M1.kYU * y(1) – M1.kEY * y(2);
dZ = M1.kZY * y(2) – M1.kEZ * y(3);
dX = M1.kX * y(4) – M1.kEX * y(4) …
– M1.etaEY * y(2) .* y(4) .* exp(-M1.lambda * y(3));

dydt = [dU; dY; dZ; dX];
end
The code runs without errors and the qualitative behaviour is reasonable, but the resulting trajectories do not match the published figures, even though the ODEs, parameters, dose times, and mass values were reproduced as closely as possible from the article.
Since all variables depend on U(t), any mismatch in handling the discontinuities propagates through the entire system.
I am unsure whether the mismatch is due to solver handling of the discontinuities, or to something subtle in the way I am implementing the impulses.
I would really appreciate if someone could look at this structure and tell me if I am handling the Dirac-type impulses and the splitting of the integration in a correct way for ode45, or if there is anything obviously wrong or missing in this implementation that could explain why my simulated curves do not match the published plots.
Thank you in advance for any help or suggestions.Hello everyone,
I am trying to reproduce the numerical results of a published tumour–immune ODE model that includes Dirac-type impulsive inputs. The system has four coupled variables [U,Y,Z,X][U, Y, Z, X][U,Y,Z,X]. Between injections, the dynamics are continuous, and at each injection time an instantaneous jump is applied to UUU:
Continuous part:
dU/dt = –kYU·U
dY/dt = kYU·U – kEY·Y
dZ/dt = kZY·Y – kEZ·Z
dX/dt = kX·X – kEX·X – η·Y·X·exp(–λ·Z)
At each dose time t=ti:
U(ti+)=U(ti−)+dosei×massi​
I have implemented this in MATLAB using two scripts:
A main script, where I loop over the dose times, call ode45 between impulses, and apply the jump in UUU. For example:
% Mouse 1 parameters and initial conditions
y_curr = [M1.U0; M1.Y0; M1.Z0; M1.X0];
t_curr = M1.t0;

for i = 1:length(tDose1)
ti = tDose1(i);
if ti > M1.tend, break; end

% integrate continuously until next dose
[t_seg, Y_seg] = ode45(@(t,y) fM1_fun(t,y,M1), [t_curr ti], y_curr, options);

% store results (code omitted here)
if isempty(t_all1)
t_all1 = t_seg;
Y_all1 = Y_seg;
else
t_all1 = [t_all1; t_seg(2:end,:)];
Y_all1 = [Y_all1; Y_seg(2:end,:)];
end

% state just before the impulse
y_curr = Y_seg(end,:)’;

% apply Dirac-type impulse to U
Mi = mass_fun1(ti); % interpolated mass at ti
ui = uDose1(i); % dose at ti (mg/kg)
y_curr(1) = y_curr(1) + ui * Mi; % jump in U

t_curr = ti;
end
2. ODE function files, one for each parameter set (Mouse 1 and Mouse 2). For example, for Mouse 1:
function dydt = fM1_fun(t,y,M1)

dU = -M1.kYU * y(1);
dY = M1.kYU * y(1) – M1.kEY * y(2);
dZ = M1.kZY * y(2) – M1.kEZ * y(3);
dX = M1.kX * y(4) – M1.kEX * y(4) …
– M1.etaEY * y(2) .* y(4) .* exp(-M1.lambda * y(3));

dydt = [dU; dY; dZ; dX];
end
The code runs without errors and the qualitative behaviour is reasonable, but the resulting trajectories do not match the published figures, even though the ODEs, parameters, dose times, and mass values were reproduced as closely as possible from the article.
Since all variables depend on U(t), any mismatch in handling the discontinuities propagates through the entire system.
I am unsure whether the mismatch is due to solver handling of the discontinuities, or to something subtle in the way I am implementing the impulses.
I would really appreciate if someone could look at this structure and tell me if I am handling the Dirac-type impulses and the splitting of the integration in a correct way for ode45, or if there is anything obviously wrong or missing in this implementation that could explain why my simulated curves do not match the published plots.
Thank you in advance for any help or suggestions. Hello everyone,
I am trying to reproduce the numerical results of a published tumour–immune ODE model that includes Dirac-type impulsive inputs. The system has four coupled variables [U,Y,Z,X][U, Y, Z, X][U,Y,Z,X]. Between injections, the dynamics are continuous, and at each injection time an instantaneous jump is applied to UUU:
Continuous part:
dU/dt = –kYU·U
dY/dt = kYU·U – kEY·Y
dZ/dt = kZY·Y – kEZ·Z
dX/dt = kX·X – kEX·X – η·Y·X·exp(–λ·Z)
At each dose time t=ti:
U(ti+)=U(ti−)+dosei×massi​
I have implemented this in MATLAB using two scripts:
A main script, where I loop over the dose times, call ode45 between impulses, and apply the jump in UUU. For example:
% Mouse 1 parameters and initial conditions
y_curr = [M1.U0; M1.Y0; M1.Z0; M1.X0];
t_curr = M1.t0;

for i = 1:length(tDose1)
ti = tDose1(i);
if ti > M1.tend, break; end

% integrate continuously until next dose
[t_seg, Y_seg] = ode45(@(t,y) fM1_fun(t,y,M1), [t_curr ti], y_curr, options);

% store results (code omitted here)
if isempty(t_all1)
t_all1 = t_seg;
Y_all1 = Y_seg;
else
t_all1 = [t_all1; t_seg(2:end,:)];
Y_all1 = [Y_all1; Y_seg(2:end,:)];
end

% state just before the impulse
y_curr = Y_seg(end,:)’;

% apply Dirac-type impulse to U
Mi = mass_fun1(ti); % interpolated mass at ti
ui = uDose1(i); % dose at ti (mg/kg)
y_curr(1) = y_curr(1) + ui * Mi; % jump in U

t_curr = ti;
end
2. ODE function files, one for each parameter set (Mouse 1 and Mouse 2). For example, for Mouse 1:
function dydt = fM1_fun(t,y,M1)

dU = -M1.kYU * y(1);
dY = M1.kYU * y(1) – M1.kEY * y(2);
dZ = M1.kZY * y(2) – M1.kEZ * y(3);
dX = M1.kX * y(4) – M1.kEX * y(4) …
– M1.etaEY * y(2) .* y(4) .* exp(-M1.lambda * y(3));

dydt = [dU; dY; dZ; dX];
end
The code runs without errors and the qualitative behaviour is reasonable, but the resulting trajectories do not match the published figures, even though the ODEs, parameters, dose times, and mass values were reproduced as closely as possible from the article.
Since all variables depend on U(t), any mismatch in handling the discontinuities propagates through the entire system.
I am unsure whether the mismatch is due to solver handling of the discontinuities, or to something subtle in the way I am implementing the impulses.
I would really appreciate if someone could look at this structure and tell me if I am handling the Dirac-type impulses and the splitting of the integration in a correct way for ode45, or if there is anything obviously wrong or missing in this implementation that could explain why my simulated curves do not match the published plots.
Thank you in advance for any help or suggestions. #bioinformatic #ode #rungekutta #dirac MATLAB Answers — New Questions

​