University of Sheffield
School of Mathematics and Statistics
MAS286 Mathematics and Statistics in Action
What is pharmaco-kinetics?
Pharmaco-kinetics is concerned with how pharmaceutical drugs move around the body. In
reality this is a highly complex process, as drugs are absorbed, distributed, metabolised and
eliminated through various parts of the body. We will make various assumptions that mean we
can simplify these processes. In particular we will assume the body is made up of a small
number of ‘compartments’ through which the drugs move. We will also simplify the processes
by which drugs are administered. Our focus will be on how the concentration of a particular
drug in the body changes over time. We will thus construct mathematical models using
ordinary differential equations that describe how this occurs. The models will be made up of
functions that describe the key mechanisms that we think are involved. While data will guide
our knowledge of these mechanisms as well as ‘sensible’ parameter values to use, here we will
not be looking to fit trends to data-sets, or what we might call statistical modelling. Instead,
we will construct mechanistic mathematical models and analyse them.
Some details
In this topic you will see how we can build models to represent pharmaco-kinetics using sets of
first-order ordinary differential equations. The mathematical techniques will thus build on
content from MAS110 and MAS111. There will also be some programming in Python. I do not
assume any prior knowledge of Python, since some of you have never met it before. Sample
code that you can edit is available on Blackboard, as well as some instructions for those of you
who feel unfamiliar with Python. There will be some medical/biological terminology you will
need to learn, but this is secondary to the mathematical techniques. The overarching
methodology I am trying to introduce to you is mathematical modelling - specifying a
problem, translating it in to a mathematical model, analysing the model, and translating back.
1
1 Single intravenous bolus dose
We will start with the simplest (yet realistic) model we can think of, before gradually building
more complexity into it. Let us assume that:
• the body can be considered as a single compartment (effectively the bloodstream);
• there is a rapid dose of the drug into the body, called an intravenous bolus;
• only one dose of the drug is given;
• the drug is eliminated from the body at a rate proportional to the drug concentration.
Let C(t) be the concentration of the drug in the bloodstream at time t. If the rate of
elimination of the drug is k, we can write
dC
dt = −kC (1.1)
for t > 0, with some initial concentration C(0) = C0 > 0. Note that we only have a negative
term on the right-hand side of equation (1.1), indicating that the concentration of the drug is
always decreasing. This fits with our assumptions above, since we assume no additional drug
is added after t = 0, so all that can happen to the drug is it is eliminated from the body.
Equation (1.1) is a first-order, linear ordinary differential equation so can be solved by
separation of variables:
−kt (taking exponentials).
Imposing the initial condition C = C0 at t = 0, we find that e
A = C0, hence
C(t) = C0e
−kt
. (1.2)
We thus predict exponential decay of the drug concentration.
1.1 Calculating k
Suppose we wish to find out what the elimination rate is for a particular drug. All we need are
two data-points for the concentration at known times. One of these could be the starting
concentration (C0 at t = 0), and suppose we take a measurement of the concentration, C1
after t1 hours (hr). If we take logs of both sides of our solution (1.2) we find
ln(C(t)) = ln(C0) − kt =⇒ k =
ln(C0) − ln(C(t))
t
=⇒ k =
ln(C0) − ln(C1)
t1
(imposing C(t1) = C1). (1.3)
2
In fact, if we had multiple data points, we would simply need to know the gradient of the line
of ln(C(t)) vs t. Since k is a rate, it has units of hr−1
. Elimination rates for real drugs tend to
be in the range k ∈ [0.02, 0.4] hr−1
.
1.2 A note on C0
We assumed above that we know the initial concentration, C0. This might seem obvious since
we know the dose, but dose and concentration are not the same thing. The dose is the actual
amount of the drug administered, while the concentration is the relative amount of drug in
the body, as a proportion of the effective volume of the body. We can convert between these
two quantities by noting that
C0 =
drug dose
effective volume.
For complex biological reasons, the effective volume is not necessarily equal to a patient’s
actual body (or blood) volume, and varies from drug to drug. For simplicity, however, we will
assume it is a fixed value for each drug in every patient.
1.3 Drug half-life
Given that we expect exponential decay of the drug concentration over time, we cannot give a
precise time when the concentration reaches zero. However, we can calculate a ‘half-life’ for
the drug based on its elimination rate. If we know that the concentration is Ca after ta hours,
then to find the half-life we need to know the time when the concentration is Ca/2. Using
equation (1.3), we have
ln(Ca/2) = ln(Ca) − kt1/2 =⇒ t1/2 =
ln(Ca) − ln(Ca/2)
k
=⇒ t1/2 =
ln(2)
k
. (1.4)
1.4 An example problem
Suppose the half-life of a drug is known to be 3 hours, and the effective volume for the drug in
a patient is 30 litres (L). What should the initial dose be such that after 4 hours the
concentration of the drug in the patient is 2.5 mg/L?
First, we find the elimination rate k from the half-life, using equation (1.4):
t1/2 =
ln(2)
k
=⇒ k =
ln(2)
t1/2
=⇒ k = 0.231 hr−1
.
3
Next, we work out what the initial concentration must be. We know that k = 0.231 hr−1
and
that C = 2.5 mg/L at t = 4 hr. Using equation (1.2), we have
C(t) = C0e
−kt =⇒ C0 = C(t)e
kt
=⇒ C0 = 2.5e
0.231×4 = 6.297 mg/L.
The final step is to multiply the initial concentration by the effective volume to get the actual
dose, 6.297 × 30 ≈ 190 mg. A plot of the drug dynamics for these parameter values is shown
in Figure 1.
0 2 4 6 8
Time (hr)
0
2
4
6
8
Concentration (mg/L)
Figure 1: Exponential decay according to equation (1.2). Parameter values are k = 0.231 hr−1
,
C0 = 6.297 mg/L, and V = 30 L.
4
2 Repeated intravenous bolus doses
In Section 1 we assumed that a single dose of a drug is given to a patient, and studied how
the concentration of the drug decays over time. In some medical circumstances this may well
be all that happens. However, for ongoing conditions, a patient may need repeated doses of a
drug. This leads to an important question: what is the optimal time interval between doses for
a particular drug? We need to ensure the patient always has enough drug in their system so it
is effective, but not so much drug in their system that they overdose.
In Section 1 we could assume there was no drug in the bloodstream prior to the single dose.
After repeated doses, however, the drug will accumulate in the bloodstream. We know that
eventually the drug concentration decays towards zero after each dose, so if we leave enough
time between doses then there will be (effectively) no accumulation. However, if doses are
given more frequently, then some of the drug will already be present.
Consider the process demonstrated in Figure 2. A drug with a half-life of 6 hours is given to a
patient in 100 mg doses every 6 hours, with an effective volume of 25 L. The initial
concentration is thus 4 mg/L. The concentration halves to 2 mg/L after 6 hours, at which
point the second dose is given and the concentration increases to 6 mg/L. The concentration
then halves to 3 mg/L after a further 6 hours, at which point the third dose is given and the
concentration increases to 7 mg/L. As such the drug is accumulating in the bloodstream.
Notice, however, that this rate of accumulation starts to slow, and eventually the drug
concentration simply fluctuates between 8 mg/L and 4 mg/L. There is thus a limit to
accumulation, because eventually the constant amount of drug that is added at each dose
becomes equal to the density-dependent amount that is removed over the next 6 hours. We
can therefore design our dosage regime to ensure that the patient always has an effective, but
safe, amount of drug in their system.
Concentration (mg/L)
Figure 2: Repeated doses according to equation (2.13). Parameter values are k = 0.116 hr−1
,
C0 = 4 mg/L, V = 25 L, and τ = 6 hr.
5
In fact, this process is given by a geometric series. Recall from Section 1 that the dynamics of
a single dose are given by
C(t) = C0e
−kt
. (2.1)
Suppose that a dose of the drug is given every τ hours. We will also label the concentrations
with a subscript to the dose number. Thus τ hours after the first dose, but immediately before
the second dose, the drug concentration is
C1(τ ) = C0e
−kτ
. (2.2)
The second dose is then given. Immediately after this second dose, re-starting time as t = 0,
the drug concentration is
C2(0) = C1(τ ) + C0 = C0(1 + e
−kτ ). (2.3)
This quantity is therefore the new ‘initial condition’, and the dynamics for the next τ hours are
given by
C2(t) = C0(1 + e
−kτ )e
−kt
, (2.4)
and exactly τ hours after this second dose, but just before the third, we have
C2(τ ) = C0(1 + e
−kτ )e
−kτ
. (2.5)
This pattern continues. For example, immediately and τ hours after the third dose we have,
respectively,
C3(0) = C0(1 + e
−kτ )e
−kτ + C0 = C0(1 + e
−kτ + e
−2kτ ), (2.6)
C3(τ ) = C0(1 + e
−kτ + e
−2kτ )e
−kτ
. (2.7)
Let R = e
−kτ , then we can re-write (2.6)–(2.7) as
C3(0) = C0(1 + R + R
2
), (2.8)
C3(τ ) = C0(R + R
2 + R
3
). (2.9)
Extrapolating, the concentrations immediately and τ hours after the n-th dose are given by
Cn(0) = C0(1 + R + · · · + R
n−1
) = C0
Xn
i=1
R
i−1
, (2.10)
Cn(τ ) = C0(R + R
2 + · · · + R
n
) = C0
Xn
i=1
R
i
. (2.11)
These are clearly geometric series. We can express (2.10) as a simple fraction as follows:
Cn(0) = C0(1 + R + · · · + R
n−1
) =⇒ RCn(0) = C0(R + R
2 + · · · + R
n
)
=⇒ Cn(0) − RCn(0) = C0(1 − R
n
)
=⇒ Cn(0) = C0
1 − Rn
1 − R
= C0
1 − e
−nkτ
1 − e
−kτ
.
6
By similar reasoning, we find that
Cn(τ ) = C0R
1 − Rn
1 − R
= C0e
−kτ
1 − e
−nkτ
1 − e
−kτ
. (2.12)
This also leads to an expression for the concentration at any time, t, after the n-th dose,
Cn(t) = C0e
−kt
1 − e
−nkτ
1 − e
−kτ
. (2.13)
We can use these expressions to find the minimum and maximum concentrations of the drug
after many doses, i.e. as n → ∞. In particular, notice that e
−nkτ → 0 as n → ∞, meaning
we have
C∞(0) = C+ =
C0
1 − e
−kτ
for the maximum concentration, and
C∞(τ ) = C− =
C0e
−kτ
1 − e
−kτ
for the minimum.
2.1 An example problem
Suppose that a drug has a half-life of 4 hours, and is given in 200 mg doses every 6 hours with
effective volume of 25 L. Clinical guidelines suggest that, in the long-term, the maximum safe
concentration of the drug is 20 mg/L, and a minimum concentration of 7.5 mg/L is needed
for the drug to be effective. Does this drug dosage regime fit within these guidelines?
First, recalling Section 1.2, we calculate the initial concentration C0 = 200/25 = 8 mg/L.
Next, we use equation (1.4) to calculate the elimination rate k = ln 2/4 = 0.173 hr−1
. This
means we can find e
−kτ = 0.354. Substituting these values and τ = 6 hr into the expressions
for C+ and C−, we obtain
C+ =
8
1 − 0.354
= 12.38 mg/L
for the maximum concentration, and
C− =
6 × 0.354
1 − 0.354
= 4.38 mg/L
for the minimum concentration. Therefore this dosage regime is always safe (since the
maximum is less than 20 mg/L) but not always effective (since considerable time will be spent
at concentrations lower than 7.5 mg/L).
Now let’s think about the problem in a different way. Suppose we wish to design our dosage
regime so that it fits the clinical guidelines exactly - that is, a maximum concentration of 20
7
mg/L and a minimum of 7.5 mg/L. What should the dosage amount and time interval be?
We already know that k = 0.173 hr−1
. First, notice that
C∞(τ )
C∞(0) = R = e
−kτ =⇒
7.5
20
= e
−0.173τ
.
By taking logs and re-arranging this we get τ = − ln(0.375)/0.173 = 5.67 hr. Let’s round this
to the more realistic time frame of 6 hours. Then, by considering the maximum, we have
C0 = C∞(0)(1 − e
−kτ ) =⇒ C0 = 20(1 − e
−0.173×6
)
=⇒ C0 = 12.92 mg/L.
Given that V = 25 L, this means a dose of 322.92 mg, which we might round down (because
we definitely wish to avoid overdosing) to the more realistic 300 mg. Thus, based on our
knowledge of the drug, we would recommend a dose of 300 mg every 6 hours for maximum,
yet safe, efficacy. (A quick check shows that our rounding means this actually gives a
maximum of 18.58 mg/L and a minimum of 6.58 mg/L.) A time course of this regime is
shown in Figure 3(a).
2.2 Loading dose and maintenance dose
A downside of the approach we have just described is that it might take a long time for the
drug concentration to reach these long-term maximum and minimum values. For example,
Figure 3(a) shows that after the first dose the concentration is ineffective for over 3 hours.
What we might consider doing, then, is to initiate treatment with a larger ‘loading dose’ that
takes the concentration to (or near) the maximum concentration, with following dosages given
as a ‘maintenance dose’. The loading dose should come as close as possible to immediately
reaching the maximum concentration. For our example above, then, we might take a loading
dose of 20 × 25 = 500 mg. The maintenance dose and timing would then be identical to
previously (300 mg every 6 hours) since this regime is already balanced to fluctuate between
the maximum and minimum concentrations. Plots comparing the time-course of the
concentrations with and without the loading dose are shown in Figure 3.
We can still write a general solution, similar to equation (2.13), by considering the dynamics
over subsequent doses. Calling the loading dose CL and the maintenance dose CM, some time
t after the first dose we have
C1(t) = CLe
−kt
. (2.14)
Immediately after the second dose at time τ , re-starting time as t = 0, we have
C2(0) = C1(τ ) + CM = CM + CLe
−kτ
, (2.15)
and the dynamics for the next τ hours will now be given by
C2(t) = (CM + CLe
−kτ )e
−kt
. (2.16)
8
0 4 8 12 16 20 24 28 32 36
Time (hr)
0
10
20
Concentration (mg/L)
(a)
0 4 8 12 16 20 24 28 32 36
Time (hr)
0
10
20
Concentration (mg/L)
(b)
Figure 3: Repeated doses according to (a) equation (2.13) with no loading dose, and (b)
equation (2.20) with a loading dose. In each case, parameter values are k = 0.173 hr−1
,
C0 = CM = 12 mg/L, CL = 20 mg/L, V = 25 L, and τ = 6 hr. The maximum and minimum
concentrations given by the clinical guidelines are shown as red dashed lines.
Then, immediately after the third dose we have
C3(0) = CM + (CMe
−kτ + CLe
−2kτ ). (2.17)
Note that we can then re-write this as
C3(0) = CM(1 + e
−kτ + e
−2kτ ) + (CL − CM)e
−2kτ
. (2.18)
and for any time t after this we would have
C3(t) = CM(1 + e
−kτ + e
−2kτ )e
−kt + (CL − CM)e
−2kτ e
−kt
. (2.19)
We can see the first term is exactly the same as we had previously and will again be given by a
geometric sum. There is then the addition of a second term to do with the extra loading dose.
The final expression for the dynamics in this case is thus given by
Cn(t) = CMe
−kt
1 − e
−nkτ
1 − e
−kτ
+ (CL − CM)e
−k((n−1)τ+t)
. (2.20)
9
3 Single oral dose
So far we have assumed that as soon as the drug is administered it is immediately present in
the bloodstream. This is true for drugs that are administered intravenously. However many
common drugs, such as aspirin and paracetamol, are administered orally. In this case the drug
first reaches the gastrointestinal (GI) tract where it dissolves and is gradually absorbed in to
the bloodstream. To model the drug’s dynamics we need to include two ‘compartments’ - one
for the amount of the drug in the GI tract (which we shall denote by XG), and one for the
amount of the drug in the bloodstream (XB). Notice that we use the word amount here, as
the concentration would have little meaning in the GI tract.
Assume, as before, that there is a single dose of the drug. Then the amount of the drug in the
GI tract simply reduces as it is absorbed in to the bloodstream. Hence, we can describe the
dynamics of this first compartment by
dXG
dt = −aXG. (3.1)
For the second compartment we can assume all the drug that leaves the GI tract arrives in the
bloodstream, and this then decays as before as the body uses up the drug. Hence, the second
equation will be
dXB
dt = V
dC
dt = aXG − kXB = aXG − kV C. (3.2)
Even before deriving the solutions to these equations, we can get a good qualitative picture of
what might happen. The amount of drug in the GI tract will clearly decay exponentially down
to zero. Initially, the drug concentration in the bloodstream will be very low, suggesting little
being lost due to decay, but the amount being absorbed in to the bloodstream would be
relatively high, meaning the concentration will initially increase. As time goes on, the amount
of drug left in the GI tract will decrease until a point is reached that the absorption of new
drug is less than the decay of existing drug, and the bloodstream concentration will reduce.
Eventually we would expect the concentration to approach zero.
We can solve this pair of equations in two parts. The first quite simply yields
XG(t) = XG(0)e
−at
. (3.3)
We can then substitute this in to the second equation, to give
dXB
dt = V
dC
dt = aXG(0)e
−at − kV C. (3.4)
This can be re-arranged into the form
dC
dt + kC =
aXG(0)
V
e
−at
. (3.5)
Written in this form, we can see it is possible to solve this equation using the integrating
10
0 4 8 12 16 20 24
Time (hr)
0
2
4
6
8
Concentration (mg/L)
Figure 4: Concentration in the bloodstream after an oral dose, according to equation (3.7).
Parameter values are a = 1 hr−1
, k = 0.1 hr−1
, XG(0) = 250 mg, and V = 25 L.
factor e
kt, giving
e
kt dC
dt + e
ktkC =
aXG(0)
V
e
−ate
kt =⇒
d
dt
e
ktC
=
aXG(0)
V
e
kt−at
=⇒ e
ktC =
aXG(0)
V (k − a)
e
kt−at + A (for some constant A)
=⇒ C(t) = aXG(0)
V (k − a)
e
−at + Ae−kt
. (3.6)
We know that the initial concentration of the drug is zero. This allows a simple calculation to
find the constant A, meaning our solution is
C(t) = aXG(0)
V (k − a)
(e
−at − e
−kt). (3.7)
A plot of the concentration over time for some chosen parameter values is shown in Figure 4.
Notice that in this case the maximum possible concentration in the bloodstream would be
XG(0)/V = 10 mg/L but, due to the balance of absorption and decay, this value is never
reached. What, then, is the maximum concentration? This peak is where the curve of the
concentration becomes flat, which is by definition when dC/dt = 0. Hence we can find when
it occurs by finding
aXG(0)
V
e
−at − kC = 0 =⇒
aXG(0)
V
e
−at −
k
k − a
(e
−at − e
−kt)
= 0
=⇒ ke−kt − ae−at = 0
=⇒ ln
k
a
= (k − a)t
=⇒ t = ln
k
a
1
k − a
. (3.8)
For the parameter values used in Figure 4, expression (3.8) gives a value of t = 2.56 hr, so
just after 2.5 hours. By substituting this value of t back into our solution (3.7), we find that
this gives a maximum concentration of 7.74 mg/L.
1
0 4 8 12 16 20 24
Time (hr)
0
2
4
6
8
10
Concentration (mg/L)
a = 2 hr−1
a = 0.2 hr−1
a = 0.02 hr−1
Figure 5: Concentrations in the bloodstream after an oral dose for varying absorption rates.
Parameter values are a ∈ {0.02, 0.2, 2} hr−1
, k = 0.1 hr−1
, XG(0) = 250 mg, and V = 25 L.
In Figure 5 we compare the curves when the absorption rate of the drug is varied. Clearly the
faster the drug is absorbed in to the bloodstream, the faster the concentration increases and
the higher levels it reaches. However, this comes at the cost of faster drug loss. Plugging the
values of a into expression (3.8), we find the peak concentration moves from 1.6 hours at high
absorption rates, to 6.9 hours for medium rates and 20.1 hours for slow rates of absorption.
3.1 An example problem
A clinician is prescribing an orally delivered drug to a patient and wants to know how long it
will be before there is a higher amount of drug in the bloodstream than in the GI tract. The
drug is given in a 200 mg dose and the patient has an effective bloodstream volume, V = 20
L. The absorption rate is a = 0.5 hr−1
and the elimination rate is k = 0.1 hr−1
.
When plotted as functions of time, the amounts of drug in the bloodstream, XB(t), and GI
tract, XG(t), will be represented by two curves. The question we need to ask is when these
two curves cross. Before this point there is more drug in the GI tract, and after this point
there is more drug in the bloodstream. Notice that we are looking at amounts here, meaning
whereas we previously found a formula for the concentration C(t), we now need to consider
the actual amount XB(t) = C(t)V . The point when the two curves cross will then be given by
XG(t) = XB(t) =⇒ XG(0)e
−at =
aXG(0)
(k − a)
(e
−at − e
−kt). (3.9)
Cancelling the XG(0) terms and re-arranging this gives
k − 2a
k − a
e
−at =
−a
(k − a)
e
−kt
. (3.10)
12
Now we can cancel the (k − a) terms and again re-arrange to find
e
−at
e
−kt =
−a
k − 2a
. (3.11)
We can then take logs of both sides and re-arrange to find
t =
1
k − a
ln
−a
k − 2a
(3.12)
Substituting in the values we are given, we obtain t = 1.47 hours.
Note that our expression (3.12) gives only one solution. Therefore, once the amount in the
bloodstream is greater than the amount in the GI tract, the two curves never cross again.
Note also that we cannot take the log of a negative number. Thus, if k > 2a, there will be no
solution to (3.12), indicating that there is always more drug in the GI tract than in the
bloodstream.
13
4 Repeated oral doses
In Section 3 we assumed a single oral dose. However, much like the intravenous boluses, we
often give repeated oral doses of a drug. How will the drug concentration change over time in
this case? We might expect to see curves rather like the intravenous bolus case, with
accumulation of the drug until a balance of dose and decay is achieved. In fact, for the amount
of drug in the GI tract, XG(t), we will get exactly the same equations as for the bolus case.
For the concentration of the drug in the bloodstream we can think about things as follows.
Let us assume that a dose is given every τ hours and let t be the time since the last dose was
given (meaning we always have 0 < t < τ ). Sometime after the first dose, but before the
second dose, the concentration of the drug, from (3.7), is
C(t) = aXG(0)
V (k − a)
(e
−at − e
−kt). (4.1)
Now, sometime after the second dose the concentration will be made up of some contribution
from the first dose plus a contribution from the second dose, such that
C(t) = C1(t) + C2(t) = aXG(0)
V (k − a)
(e
−a(t+τ) − e
−k(t+τ)
) + aXG(0)
V (k − a)
(e
−at − e
−kt)
). (4.2)
Factoring out some terms and noting that e
−a(t+τ) = e
−ate
−aτ gives
C(t) = aXG(0)
V (k − a)
(e
−at(1 + e
−aτ ) − e
−kt(1 + e
−kτ )). (4.3)
Extrapolating, we can see that after n doses we have
C(t) = aXG(0)
V (k − a)
(e
−at(1 + e
−aτ + · · · + e
−(n−1)aτ ) − e
−kt(1 + e
−kτ + · · · + e
−(n−1)kτ )).
(4.4)
We can see that there are two geometric series, which means we can re-write this as
C(t) = aXG(0)
V (k − a)
e
−at
1 − e
−naτ
1 − e
−aτ
− e
−kt
1 − e
−nkτ
1 − e
−kτ . (4.5)
A plot of the resulting dynamics is shown in Figure 6. Again, we clearly see the gradual
approach to a steady fluctuation between the maximum and minimum drug concentrations
after repeated doses. As before, we can find these maximum and minimum values by
considering what happens as n → ∞ and noting that the minimum concentration will always
be as the new dose is taken, i.e. t = 0, finding
C∞(0) = aXG(0)
V (k − a)
1
1 − e
−aτ
−
1
1 − e
−kτ . (4.6)
14
0 4 8 12 16 20 24 28 32 36
Time (hr)
0
10
20
30
40
Concentration (mg/L)
Figure 6: Repeated oral doses, according to equation (4.5). Parameter values are a = 1 hr−1
,
k = 0.1 hr−1
, XG(0) = 500 mg, V = 25 L, and τ = 6 hr.
Further, if we assume that by this point there is negligible new drug being absorbed in to the
body (as the previous dose has virtually run out), then we can take e
−aτ ≈ 0 and therefore
C∞(0) = aXG(0)
V (k − a)
−
e
−kτ
1 − e
−kτ
. (4.7)
Calculating the maximum concentration is rather more tricky in this case. What we can do is
to find the average concentration of the drug during one dose period (τ ), defined as the area
under the curve for one dose period at steady state,
C¯ =
1
τ
Z τ
0
Css(t) dt. (4.8)
In fact we can make things easier by noting that this is equal to the whole area under the
curve of the first dose on its own,
C¯ =
1
τ
Z τ
0
Css(t) dt =
1
τ
Z ∞
0
C1(t) dt. (4.9)
Computing this integral, we find
C¯ =
1
τ
aXG(0)
V (k − a)
Z ∞
0
(e
−at − e
−kt) dt
=
1
τ
aXG(0)
V (k − a)
−e
−at
a
+
e
−kt
k
∞
0
=
1
τ
aXG(0)
V (k − a)
1
a
−
1
k
=
XG(0)
V kτ . (4.10)
This average concentration can be useful in its own right. For instance, notice this says that a
low dose at short intervals will give an equivalent average concentration as a high dose at long
intervals. It can also be used to construct a crude approximation of the maximum
concentration using C+ = C¯ + (C¯ − C−).
15
4.1 Example problems
Suppose an oral drug is prescribed for a patient with the parameter values used in Figure 6.
Suppose that guidelines state the maximum safe concentration of the drug is 30 mg/L and the
minimum effective concentration is 10 mg/L. Let us start by finding the minimum and
maximum concentrations reached after repeated doses. From equation (4.6) we have
C− =
aXG(0)
V (k − a)
−
e
−kτ
1 − e
−kτ
(4.11)
Substituting in the values used in Figure 6, this gives C− = 27.0 mg/L. We can also find the
average concentration at steady state using
C¯ =
XG(0)
V kτ , (4.12)
which gives C¯ = 33.3 mg/L and thus C+ = 39.6 mg/L. Comparing these values to the
guidelines, it is clear that this regime would reach doses that are much too high. Let us
therefore devise a new regime.
First, note that C¯ = (C+ − C−)/2 = 20 mg/L, so from (4.10) we find that XG(0) = 50τ .
Substituting back into the equation for the minimum concentration, we find
C− = 10 =
50τ
−22.5
−
e
−0.1τ
1 − e
−0.1τ
.
This equation is difficult to solve for τ analytically, but we can solve it numerically by plotting
the right-hand side for varying τ and look for where it equals 10. Doing this, we find that the
optimum timing is τ ≈ 14 hr, which means XG(0) ≈ 700 mg, C− = 10.2 mg/L and
C+ = 29.8 mg/L. We might consider 14 hours to be unrealistic and instead take τ = 12 hr,
giving XG(0) = 600 mg, C− = 11.5 mg/L and C+ = 28.5 mg/L. These dynamics are shown
in Figure 7.
0 4 8 12 16 20 24 28 32 36 40 44 48
Time (hr)
0
10
20
30
40
Concentration (mg/L)
Figure 7: Repeated oral doses according to equation (4.5). Parameter values are a = 1 hr−1
,
k = 0.1 hr−1
, XG(0) = 600 mg, V = 25 L, and τ = 12 hr.
16
5 A two-compartment bolus model
In our models of intravenous bolus administered drugs (Sections 1–2), we assumed the body
acted as a single ‘container’ for the drug, effectively considering the drug concentration in the
bloodstream. When we looked at the orally administered drug models (Sections 3–4) we
assumed there was a second ‘compartment’ where the drug had to enter the body (the GI
tract) before it could move into the bloodstream. In that case things simplified because the
dynamics of the first compartment were very simple - exponentially decreasing to zero - and
we only needed to focus on what happened in the bloodstream.
Now we will consider a more complicated case where there are two compartments, but the
drug can now move in and out of these compartments. A biological example might be the
bloodstream as our main compartment, and tissue as the second compartment. We will
assume that the drug is administered by an intravenous bolus and so appears immediately in
the bloodstream. The concentration of the drug in the bloodstream at time t will again be
given by C(t), and the concentration in the tissues will be given by D(t). While in the
bloodstream the drug is used up and is eliminated at rate k. The drug can also pass into
tissues at rate a and can pass back from tissues to the bloodstream at rate b. We will assume
the drug is not eliminated while in the tissue. Our mathematical model thus takes the form
dC
dt = −kC − aC + bD, (5.1)
dD
dt = aC − bD. (5.2)
Since both equations explicitly depend on both variables, in this case we cannot make a
simplification like we made in Sections 3–4. However, solving systems of linear ordinary
differential equations is still a relatively straightforward task.
First note that we can write equations (5.1)–(5.2) in matrix form, that is,
d
dt
C
D
!
=
−k − a b
a −b
! C
D
!
.
Let us seek solutions of the form
C(t)
D(t)
!
=
u1
u2
!
e
λ1t +
v1
v2
!
e
λ2t
(5.3)
for some λ1, λ2. We now see that this has become an eigenvalue/eigenvector problem since,
substituting (5.3) into our matrix equation we have
λ1
u1
u2
!
e
λ1t =
−k − a b
a −b
! u1
u2
!
e
λ1t
, (5.4)
and similarly for λ2.
17
To solve such a problem we first find the eigenvalues. Calling the matrix of parameters M, we
do this by setting
det(M − λI) = 0. (5.5)
This leads to the characteristic equation
λ
2 + (k + a + b)λ + kb = 0. (5.6)
This does not factorise particularly nicely, so let’s simplify things a little by choosing constants
α and β such that
α + β = k + a + b, (5.7)
αβ = kb. (5.8)
This helps as it means the characteristic equation is now
λ
2 + (α + β)λ + αβ = (λ + α)(λ + β) = 0, (5.9)
for which the solutions are clearly just λ = −α and λ = −β. If and when the time comes to
find out what values α and β take, we would need to substitute them into the quadratic
formula to find
α =
(α + β) + p
(α + β)
2 − 4αβ
2
,
β =
(α + β) −
p
(α + β)
2 − 4αβ
2
. (5.10)
We have therefore narrowed down our general solution to be