administer drug dose at certain time step in octave

Question

I am doing a simple task on ocatve. I have to administer drug dose at day 1,7,21 and 28..

i wrote a function like that:

   function xdot = f (x,t)

   a=[1;7;21;28]
   drug=0; ### initially drug is zero

   if (t==a)
   drug=drug+57.947;

  else
  drug=drug+0;
  endif

  xdot=(-0.4077)*(x)+ drug; passing the value of drug to differential equation

  endfunction

In the main file i called this function in lsode:

  t=linspace(0,30,30);
  x0=0;
  y=lsode(@ex,x0,t); ### ex is the file name where function is written
  plot(t,y,'o')

This program doesn't work.. it displays all the time zero value for drug. Can anybody help me that how to administer dug dose with certain time step by manipulating linspace function.

Answer 1

It looks like you have a simple clearance model, and at each time in a, you want the dose to be delivered instantaneously. That is, at each time in a, the amount of drug in the subject increases by 57.947.

If that is the model that you have in mind, implementing it in the formula for xdot will not work well. You would actually need to implement it as a "delta function", and lsode won't work with that.

Instead, I recommend solving the problem in stages, corresponding to the time intervals [0, 1], [1, 7], [7, 21], and [21, 28]. During one stage, the differental equation is simply xdot = -0.4077*x. In the first stage, the initial condition is 0. In the next stage, the initial condition is the final value of the previous stage plus the dosage amount 57.947.

Here's a script:

dose = 57.947;
a = [1 7 21 28 30];

x0 = 0.0;
t0 = 0;
t = [];
y = [];
for t1 = a
    tinterval = linspace(t0, t1, 2*(t1 - t0) + 1);
    yinterval = lsode(@(x, t) -0.4077*x, x0, tinterval);
    t = [t tinterval];
    y = [y yinterval'];
    t0 = t1;
    x0 = yinterval(end) + dose;
end

plot(t, y, "linewidth", 2);

The script creates this plot:

---

Note that the differential equation xdot = -k*x has the solution x0*exp(-k*(t-t0)), so the call to lsode could be replaced with

    yinterval = x0*exp(-0.4077*(tinterval - t0));

If you do this, also remove the transpose from yinterval two lines below:

    y = [y yinterval];

---

If you want to keep the implementation of the administration of the drug dose within the formula for xdot, you'll need to distribute it over a small time interval. It could be implemented as a short rectangular pulse of width w and height 57.974/w. You'll also need to ensure that lsode takes sufficiently small internal time steps (less than w) so that it "sees" the drug dose.

Answer 2

You probably want replace t==a by ismember (t, a). Also, you should erase the else clause as it has no effect on the answer.

---

UPDATE

Consider rewriting your function as:

function xdot = f (x,t)
  xdot = -0.4077*x + 57.947*ismember (t, [1 7 21 28])
endfunction