A model for an intramuscular injection of a medication.

Needs["Graphics`Colors`"];

Intravenous injection of a medication:

Injected drug in blood decays exponentially due to metabolic degradation and elimination

Differential equation model:
c'[t] = -β c[t]
c[0] = 50 mg

Unique solution:  c[t] = 50 exp[-βt]

β  =  -Log[ .9]

Clear[soln0,  c]
soln0  =  
            DSolve[{
c'[
                    t]  ==  -β c[t],  
c[
                    0]  ==  50},  
c[t],  t];

c[t_]  =  (c[t] /. soln0)[[1]]

Plot[{ c[t]},  {t,  0,  24},  
        PlotStyle -> {{Thickness[ .01], Blue}}]

Suppose we measured blood drug levels during a 24 hour period, and we saw a curve c[t]  which looked like our model curve.    We would like to "fit" the experimentally observed drug kinetics by our model .  How can we find the right parameter β from  the experimental data?

Notice the logarithm of c[t] has a straight plot.

Plot[{ Log[c[t]]},  {t,  0,  24},  
        PlotStyle -> {{Thickness[ .01], Blue}}]

We can recover the parameter β from the slope of  Log[c[t]]

h[
          t_]  =  (∂_t)(Log[c[t]]);

    h[10]

    h[20]

Intramuscular injection of a medication:

Model description:
Drug injected into muscle.
Two stage process occurs:

1)  exponential release of drug from muscle into bloodstream, with rate constant α

2)  exponential decay of drug in blood as before, with different rate constant β

Differential equation model:
A[t] = drug in muscle
B[t] = drug in bloodstream
A[0] = 50
B[0] = 0
A'[t] =  -α A[t]
B'[t] =  α A[t]  - β B[t]

This can be solved explicitly.  We'll let Mathematica do the work.

There are two parameters to specify,  α and β.  We give them some arbitrary values.

We specify and solve the differential equation model.

Clear[α,  β]
Clear[A,  B]

soln  =  
          DSolve[{
A'[
                  t]  ==  -α A[t],  
B'[
                  t]  ==  α A[t]  -  β B[
                      t], 
A[0]  ==  
                50,   
B[0]  ==  
                0},  
{A[t],  B[t]},  
            t];

A[t_]  =  
        A[t] /. soln[[1]];
B[
          t_]  =  Expand[
            B[t] /. soln][[1]];
α  =  -Log[ .2]
β   =  -Log[ .9]

A[t]

B[t]

Note that the solution B[t] to the model has the form
K ((e^(-β t)) -  
              (e^(-α t))), for some constant K which must be somehow related to the data of the problem, namely the 50 mg dose and the parameters α and β.

In fact the constant is  K = ((50 α)/(α - β)).

This observation is not at all obvious from the solution to the specific model with specified values of α and β,  but it follows from solving the model with unspecified parameters.  We are going to use this fact later.

((50 α)/(α - β))

We plot the solution:  
The green curve representents A[t],  drug in muscle.
The red curve represents B[t],  drug in blood.
The blue curve represents c[t],  the solution to the primitive model.

Plot[{A[t],  B[t],  c[t]},  {t,  0,  24},  
          PlotStyle -> {{Thickness[ .01], Green}, {Thickness[ .01], 
                Red},   {Thickness[ .01], Blue}}]

Let's pose two questions concerning our model:  
1.  Can we find the maximum level of the drug in the blood stream?
2.  Could we recover the parameters α and β from experimental data?

Maximum drug level in the bloodstream

What is the maximum value of drug in the bloodstream?      It has to occur where the derivative of the function B[t] is equal to  0.   So we find the maximum by finding the solution to B'[t] = 0.

B[t_]  =  Expand[(B[t] /. soln)[[1]]]

DB[t_]  :=  B'[t]

The following cell solves the equation 
B[t] =0
by a numerical method.  We will discuss such methods (using calculus) later.

maxplace  =  FindRoot[DB[t]  ==  0,  {t,  1}]

We compute the value of B[t] where it is maximum:

B[t] /.  maxplace

Our conclusion:  the maximum amount of medication in the blood occurs at t = 1.8126 hours, and the maximum amount is 41.308 mg.

Parameter estimation

Suppose we measured blood drug levels during a 24 hour period, and we saw a curve B[t]  which looked like the red curve.    We would like to try to model the drug kinetics by our model with parameters α and β.  How can we find the parameters?

In our model, the red curve looks very similar to the blue curve for large time, 
that is, the solution to our model is K ((e^(-β t)) -  
              (e^(-α t))), and for large t, the second exponential term is tiny in comparison with the first.    So the data looks like K (e^(-β t))  when we only look at large t.   We could find β by taking the derivative of the logarithm  of the data curve B[t] for large t.

Note that the graph of Log[B[t]] looks quite straight for large t, as we expect from the discussion above.

Plot[Log[B[t]],  {t,  0,  
        24}, 
PlotStyle -> {{Thickness[ .01], 
            Red},  {Thickness[ .01], Blue}}]

We can recover the parameter β from the slope of the curve Log[B[t]] for large t.

h[t_]   =  (∂_t) Log[B[t]];
Plot[
      h[t],  {t,  0,  
        24}, 
PlotStyle -> {{Thickness[ .01], 
            Red},  {Thickness[ .01], Blue}}]

    
h[10]
h[20]

One we know β, we can also find  K since the curve B[t]  looks like  
y = K (e^(-β t)) 
for large t.   Thus the logarithm of B[t] looks like
y =  Log[K] - β t
 for large t.    So if we  compute the tangent line to Log[B[t]] at some large t, this tangent line will also look like
 y =  Log[K] - β t.
 We can find   Log[K] as the y-intercept of this tangent line, and then find K as Exp[] of the y-intercept.
 
 Here is the tangent line:

ylin[t_]  :=   Log[B[20]]  +  
        h[20] (t - 20)

For illustration, we plot the tangent line together with the curve Log[B[t]].

Plot[{Log[B[t]],  ylin[t]},  {t,  0,  
        24}, 
PlotStyle -> {{Thickness[ .01], 
            Red},  {Thickness[ .01], 
            Blue}}]

We compute the y-intercept of the tangent line, and Exp[]  of this y-intercept.
We see that we recover the constant K = 53.5025.

ylin[0]

    Exp[ylin[0]]

Let's also recall that  K = 53.5025 =  ((50 α)/(α - β)).  Since we now know β and this constant K, we can solve for α.  (This we can do by hand.)

Effect of changing the parameters:

There are two parameters to specify,  α and β.  We give them some different arbitrary values.

We specify and solve the differential equation model.

Clear[α,  β]
Clear[A,  B]
soln  =  
        DSolve[{
A'[
                t]  ==  -α A[t],  
B'[
                t]  ==  α A[t]  -  β B[
                    t], 
A[0]  ==  
              50,   
B[0]  ==  
              0},  
{A[t],  B[t]},  
          t];

A[t_]  =  A[t] /. soln[[1]];
B[t_]  =  Expand[B[t] /. soln][[1]];

α  =  -Log[ .5]
β   =  -Log[ .9]

A[t]

B[t]

We plot the solution:  
The green curve representents A[t],  drug in muscle.
The red curve represents B[t],  drug in blood.
The blue curve represents c[t],  the solution to the primitive model.

Plot[{A[t],  B[t],  c[t]},  {t,  0,  24},  
          PlotStyle -> {{Thickness[ .01], Green}, {Thickness[ .01], 
                Red},   {Thickness[ .01], Blue}}]

Notice the qualitative change:  the red curve (quantity of medication in the blood)  has a broader and lower peak (which accords with our intuition, as the drug is "washed" out of the muscle tissue more slowly).


The first of our two questions is still interesting:  What is the the maximum drug level in the bloodstream?

Maximum drug level in the bloodstream

What is the maximum value of drug in the bloodstream?      It has to occur where the derivative of the function B[t] is equal to  0.   So we find the maximum by finding the solution to B'[t] = 0.

B[t_]  =  Expand[(B[t] /. soln)[[1]]]

DB[t_]  :=  B'[t]

maxplace  =  FindRoot[DB[t]  ==  0,  {t,  1}]

B[t] /.  maxplace

The  "intramuscular" model with variable parameters:

Let us consider the model with variable parameters α and β.

We specify and solve the differential equation model.

Clear[α, β]

Clear[A,  B]

soln  =  
          DSolve[{
A'[
                  t]  ==  -α A[t],  
B'[
                  t]  ==  α A[t]  -  β B[
                      t], 
A[0]  ==  
                50,   
B[0]  ==  
                0},  
{A[t],  B[t]},  
            t];

A[t_]  =  A[t] /. soln[[1]]

B[t_]  =  Expand[
          B[t] /. soln][[1]]
Simplify[%]

This justifies the assertion made earlier that 
B[t] = K ((e^(-β t)) -  
              (e^(-α t))),  for   K = ((50 α)/(α - β)).

The  "intramuscular" model with equal parameters:

Let us consider what happens to the model when α = β.

We specify and solve the differential equation model.

Clear[α, β]

α  =  β

Clear[A,  B]

soln  =  
          DSolve[{
A'[
                  t]  ==  -α A[t],  
B'[
                  t]  ==  α A[t]  -  β B[
                      t], 
A[0]  ==  
                50,   
B[0]  ==  
                0},  
{A[t],  B[t]},  
            t];

A[t_]  =  A[t] /. soln[[1]]

B[t_]  = B[t] /. soln[[1]]

Put in a particular value of β:

β  =  -Log[ .9]

B[t]

Plot[{ B[t],  c[t]},  {t,  0, 50}, 
        PlotStyle -> {{Thickness[ .01], Red},   {Thickness[ .01], 
              Blue}}, 
PlotRange  ->  {0, 50}]