Help implementing 4th-order Runge–Kutta for ODEs with Dirac impulses (discontinuous system)

158 views (last 30 days)
Hello
  12 Comments
Torsten
Torsten on 3 Dec 2025 at 10:58
Edited: Torsten on 3 Dec 2025 at 10:58
I copied the article and code. @Sam Chak might be interested, too. From my side, you can remove it from Google Docs.
Sam Chak
Sam Chak on 3 Dec 2025 at 14:13
Thank you for sharing the paper. However, I am unable to replicate the results. It is possible that my approximation of the Dirac Delta function is inaccurate. I only copied the data for the injected doses for Mouse 1 to test it in the code referenced in my previous comment. Please rest assured that I did not download your code.
The authors of the paper did not provide the numerical data, only displaying a representation of the data in an image (Figure 1). Moreover, I suspect that the authors used the discrete-time unit impulse sample instead of the approximation of the continuous-time Dirac Delta function .
% Mouse 1
tDose1 = [ ...
3.92040 36.48975 43.42583 61.51992 64.53560 71.77322 ...
99.51751 106.75515 120.92884 127.56333 142.03861 156.21231 ...
162.84682 170.08446 177.02051 184.25815 191.19423 198.43185 ...
205.36793 212.30398 219.54164 223.46204 226.17613 ...
];
uDose1 = [ ...
4.00602 3.86145 3.30723 0.40361 0.19879 0.27711 ...
1.78313 0.13855 1.24096 1.01205 3.25904 0.92169 ...
1.98795 0.51807 0.28916 2.21687 1.46988 1.35542 ...
1.04217 1.27108 1.72289 2.71687 0.81325 ...
];
tMass1 = [ ...
4.27873 36.67482 43.39853 61.43032 64.79217 71.51590 ...
99.63325 106.96822 120.72128 128.05625 142.11492 156.47922 ...
163.20295 170.23227 177.26162 183.98533 191.62594 198.34965 ...
205.379 212.40835 219.43765 223.41077 226.467 ...
];
mass1 = [ ...
0.0226 0.02242 0.02334 0.02242 0.02214 0.02458 ...
0.02541 0.0243 0.02605 0.02490 0.02430 0.02421 ...
0.02462 0.02394 0.02315 0.02407 0.02389 0.02306 ...
0.02343 0.02297 0.02334 0.02403 0.02274 ...
];
tAdmin = 0:0.001:230;
% approximation of a time-shifted Dirac-type impulse δ(t - tau)
sigma = 9e-3;
amp = 1/(sigma*sqrt(2*pi));
tInject = (tDose1 + tMass1)/2;
delta = zeros(numel(tInject), numel(tAdmin));
for i = 1:length(tInject)
% standard Gaussian function
% delta(i,:) = amp*exp(- 1/2*((tAdmin - tInject(i))/sigma).^2)*uDose1(i)*mass1(i);
% standard Gaussian function
delta(i,:) = exp(- 1/2*((tAdmin - tInject(i))/sigma).^2)*uDose1(i)*mass1(i);
end
uDose = sum(delta);
figure
plot(tAdmin, uDose)
xlabel('t, (Days)')
ylabel('Dose, (mg)')
title('Injected doses for Mouse 1')

Sign in to comment.

Answers (1)

Torsten
Torsten about 17 hours ago
Edited: Torsten about 14 hours ago
Look at figure 3 (a) and (b) of the article. Each time the curve for U becomes discontinuous, a dosing is applied to the mice. If you compare this to figure 1, you will observe that the number of times dosing is applied is far greater than what figure (1) seems to suggest.
Thus what I would try if I were you is to take figure (3) and read out the time points when a jump appears and the height of the jump (which corresponds to Mi*ui in your code). This gives you two new arrays: "tdosei" and "jumpheighti". Then you can use these two arrays to repeat your computation.
If this procedure also gives results that far differ from the publication, you can additionally try to apply a simpler integration scheme with a fixed value of deltat = 0.2 days (although ode45 should of course be numerically more precise). I'm not sure whether the authors also used ode45. They report about a fixed value for the time step deltat = 0.2 days, but this could also mean that ode45 was forced to make outputs every 0.2 days (and computations were performed on a much finer grid). If the authors used a fixed time stepping scheme, this could also be a reason for results that differ from yours, and you could try such a simpler integration method.
Another thing that stroke me is that you use different arrays tDosei and tMassi. In my opinion, they must be equal because the weight of the mice was measured at the times when dosing was applied.
An indication whether you are on the right track is to check the tumor volume. It shouldn't increase far beyond 200 mm^3. In your present simulations, it reaches a maximum value of almost 1,2*10^6 mm^3 (for mouse 1) which is unrealistic.
I'd be interested if you could report about your future progress (hopefully).
  1 Comment
Torsten
Torsten on 4 Dec 2025 at 0:07
Edited: Torsten on 4 Dec 2025 at 0:08
when you suggest extracting the time points and jump heights directly from Figure 3, do you mean using a tool such as DataThief (as we did for Figure 1) to retrieve both the dosing times and the corresponding values?
I think if you don't want to contact the authors of the article, you don't have a better choice, do you ?

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!