Problem involving drug delivery and ODE

I am trying to model drug delivery over time using ODE45 to simultaneously solve three differential equations, and it is giving me an exponential decay function, which is what I want. However, I am now trying to model extra doses being applied to the system. For example, I want to increase the concentration of the drug by a set amount at a certain time interval while having the exponential decay still apply. Basically, I want it to look like this:

What can I do to accomplish this?

Answer 1

[Star Strider]

See odeset and particularly the Event Location Property section. An events function will allow you to stop the integration according to specific criteria, change something, then restart the integration.

In your pharmacokinetics application, you can do this either with respect to time or drug concentration.

[Alberto]

Is there a way that you can do this by setting the criteria according to time? the documentation only describes the event happening for when y' = 0.

[Torsten]

If you know the times at which an event happens, you don't need the "events" option of the ode integrators. You can simply integrate up to time 1, change your initial conditions and/or equations, restart the solver and integrate up to time 2, change your initial conditions and/or equations,...

[Sam Chak]

Hi @Alberto

Could you please post a new question on the MATLAB Answers forum and provide the drug delivery ODE (Differential Equation) for academic purposes? You can compare the ode45 / ode15s results with the SimBiology.

When considering the dosage from an absolute time perspective, there will be a slight time delay (no matter how small) between the required drug concentration and the actual drug concentration being absorbed into the bloodstream.

If you can model such continuity in the dosage (where the curve is smooth and numerically differentiable by the ODE solver), then I believe the event() function may be unnecessary.

Update: Code is edited to illustrate the dosage model.

t   = 0:0.0001:10;
% d1  = heaviside(t - 5);   % previous code
d1  = heaviside(t - 3.5);
d2  = (tanh(2*(t - 5)) + 1)/2;
plot(t, d1), hold on
plot(t, d2), grid on, ylim([-0.5, 1.5])
legend('idealized dosage model', 'practical dosage model')

[Sam Chak]

I'm not an expert. But I was trying to illustrate that the OP's simulation (Figure 1) was too idealized, most probably due to the modeling of instantaneous drug administration (using step function or if / else conditional approach). Figure 2 depicts the effect of true pharmacodynamics where there is a slight time delay after drug administration (for every 24 hours) and before the drug plasma concentration reaches the peak.

Figure 1

Figure 2

Answer 2

[Sean de Wolski]

Have you tried using SimBiology?

Perhaps this webinar could get you started.

[Arthur Goldsipe]

Yes, SimBiology is designed just for this sort of problem. You can write your model in terms of reactions instead of ODEs and doses instead of events. If you have access to SimBiology and want to see how to translate this problem into a SimBiology model, please post another question on MATLAB Answers. It should just be a few lines of code.