代写MAS286编程、代写Python程序

日期：2023-03-14 09:23

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

C(t)

D(t)

!

=

u1

u2

!

e

−αt +

v1

v2

!

e

−βt

. (5.11)

To find the eigenvectors we substitute each eigenvalue back in and solve:

−k − a − λ b

a −b − λ

! u1

u2

!

=

0

0

!

. (5.12)

Let us take λ = −α. Using the bottom line this gives the equation

au1 + (α − b)u2 = 0 =⇒ au1 = (b − α)u2. (5.13)

Thus for some constant, κu, taking u1 = κu(b − α) and u2 = κua will always satisfy this. The

solution for λ = −β will look almost identical. We can thus write our solution for the

concentration of drug in the bloodstream and tissue as

C(t)

D(t)

!

= κu

b − α

a

!

e

−αt + κv

b − β

a

!

e

−βt

.

18

We can determine the values of the two remaining constants, κu and κv, by substituting in

initial conditions. We know that at t = 0 the patient has just received the intravenous bolus.

Therefore the concentration of the drug in the bloodstream should be at its maximum (i.e.

the dose adjusted by the effective volume, C0) and the concentration in the tissue should be

zero. These two conditions tell us

C(0) = C0 = κu(b − α) + κv(b − β) (5.14)

D(0) = 0 = κua + κva (5.15)

From (5.15) we see that we have κu = −κv. Substituting back into (5.14), this gives

C0 = κu(b − α) − κu(b − β) =⇒ κu =

C0

β − α

.

The concentration of the drug in the bloodstream is therefore given by

C(t) = C0

 α − b

α − β



e

−αt −



β − b

α − β



e

−βt

. (5.16)

This is the sum of two negative exponential functions, suggesting that the concentration

eventually tends to zero, as we should expect. An expression for the drug concentration in the

tissue can be found in a similar way.

5.1 An example problem

Suppose a 200 mg dose of a drug is given to a patient with effective volume V = 10 L. The

drug is eliminated at rate k = 0.15 hr−1

and passes between the bloodstream and tissues at

rates a = 0.2 hr−1

and b = 0.15 hr−1

. What is the drug concentration in the bloodstream 4

hours after it has been administered?

To substitute these values into our general solution, there are a few values we need to

calculate first. These are

C0 = 200/10 = 20 (recalling Section 1.2), (5.17)

α + β = 0.15 + 0.2 + 0.15 = 0.5 (using equation (5.7)), (5.18)

αβ = 0.15 × 0.15 = 0.0225 (using equation (5.8)), (5.19)

hence

α =

1

2



0.5 + √

0.5

2 − 4 × 0.0225

= 0.45, (5.20)

β =

1

2



0.5 −

√

0.5

2 − 4 × 0.0225

= 0.05, (5.21)

using equation (5.10). We can therefore write down our solution as

C(t) = 20

0.75e

−0.45t + 0.25e

−0.05t

.

All that remains is to substitute t = 4 into the above equation to find that after 4 hours the

concentration of the drug will be 6.57 mg/L. The full timecourse of this example, showing the

concentrations in both the bloodstream and the tissue, is shown in Figure 8.

1

0.0 2.5 5.0 7.5 10.0

Time (hr)

0

5

10

15

20

Concentration (mg/L)

Bloodstream

Tissue

Figure 8: Concentrations in the bloodstream and tissue according to equation (5.16). Parameter

values are a = 0.2 hr−1

, b = 0.15 hr−1

, k = 0.15 hr−1

, and C0 = 20 mg/L.

20