( previous, home, next )



There is a killer found in many household medicine cabinets -- acetaminophen, commonly known under the brand name Tylenol. Acetaminophen was discovered at Johns Hopkins in 1877, and marketed as a pain reliever beginning in 1950. US CDC estimated there were 1,600 deaths from acetaminophen poisoning between 2001 and 2010. Acetaminophen toxicity replaced viral hepatitis as the most common cause of acute hepatic failure and was the second most common cause of liver failure requiring transplantation. The United States government's Food and Drug Administration (FDA) provides recommendations for dosages and frequency of use. Is there something the FDA could do to reduce this number of deaths?

This is a classic problem in pharmacokinetics, and is a good problem with which to introduce compartmental modelling using first-order differential equation. The field of pharmacokinetics and the related field of pharmacodynamics are mature and sophisticated fields that use mathematics to merge physiology, pharmacology, and medicine for planning safe drug usage. We will only have space to touch on the basics here.

A linear rate law

Drugs are taken out of a compartment in the body by various processes like excretion and metabolism.
Let \(x(t)\) represent the bulk amount of drug in some compartment at time \(t\), where we measure time \(t\) in units of hours and the bulk amount \(x\) in units of milligrams. The rate of change in the bulk amount of drug \[dx/dt = - \text{Removal}.\] So, what formulas do we fill in for the removal?

First of all, we can assume the removal rate does not depend on time directly. The human body puts a fair amount of energy into keeping all its parts running at steady rates and constant temperatures. This is called homeostasis. Because of homeostasis, it is reasonable to suppose that the processes removing the drug via excretion or metabolization will be approximately independent of time.

Second, drugs are usually administered in small amounts relative to the size of a human body. An individual drug molecule will probably never encounter another drug molecule. So it seems like the rate at which individual drug molecules are removed from the body should be independent of the other drug molecules. If each molecule's chance of removal is the same, then halving the number of molecules should also halve the bulk output rate. Resaid as an equation, the rate of removal \(R(x)\) should satisfy \[R\left(\frac{x}{2}\right) = \frac{1}{2}R(x).\] In fact, reducing the amount by any fraction \(a\) should reduce the rate of removal by \(a\), \[R(a x) = a R(x).\]

We now have an equation for \(R\). This equation expresses a certain symmetry between the amount of drug and its removal rate. This particular kind of equation is called a functional equation because the unknown we are solving for is a function. In general, functional equations can be hard to solve, but this one isn't too bad. If we differentiate our functional equation with respect to the parameter \(a\) and evaluate at \(a = 1\), we get the differential equation \[x \frac{dR}{dx} = R.\] Using the separation of variables technique, we can solve this equation and show \[\begin{gather} \int \frac{d R}{R} = \int \frac{dx}{x}, \\ \ln {R} = C + \ln x, \\ R(x) = e^{C} x = r x, \end{gather}\] where we call the renamed integration constant \(r\) the elimination rate. The units of \(r\) will be reciprocal time. The expected time to removal is \(1/r\), while the half-life \(t _ {1/2} = (\ln 2)/r\). The elimination rate depends on the chemistry of the drug and physiology of the compartment, and is often estimated based on observational data.

A one-compartment model

The simplest pharmacokinetic model of the body is the one-compartment model -- the drug enters the body, but is then eliminated over time either through excretion or metabolization. Let \(x(t)\) represent the bulk amount of drug in the body, where we measure time \(t\) in units of hours and the bulk amount \(x\) in units of milligrams. The rate of change in the bulk amount of drug \[dx/dt = - \text{Removal} + \text{Input}\] Assuming our above result on removal and that the inputs are controlled by an arbitrary externally dosing regiment \(i(t)\), the general form of our single-compartment model is \[dx/dt = - r x + i(t) \quad \text{with $x(0)$ given.}\]

This is summarized by the flow diagram. The box represents the compartment, while the arrows represent flows in and out of the compartment.

If there is a single drug dose at the initial time \(t=0\), then the dosage can be specified in terms of just the initial condition: \[\frac{dx}{dt} = - r x, \quad x(0)\;\;\text{given}.\]

To turn these ideas into something practical, we must estimate the initial value \(x(0)\) and the elimination rate \(r\). A little research finds the following. A standard dose of acetaminophen is 325 milligrams (mg), and the plasma half-life of acetaminophen is \(1.75 \pm 0.75\) hours. A half-life \(t _ {1/2}\) is the amount of time it takes a dose to be reduced to half its initial value, and is related to the elimination rate by the formula \[e^{-r t _ {1/2}} = 1/2.\] By solving this equation, we find \(r \approx 0.4\). In order to compare our results with observations, we'll also need to convert total amount into concentrations. As a crude estimate, we'll assume our dose is administered to an adult with \(V=5\) liters of blood. Then the concentration in a blood plasma sample will be \(x(t)/V\).

[Show code]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
"""
A one-compartment model for a single dose of acetaminophen
"""

from numpy import *
from scipy import integrate
from matplotlib.pyplot import *

r1 = log(2.)/1.75
dose = 325.
V_blood = 5000. # 5 liters approximately
mmg_per_mg = 1000.

def model1(X, t):
    return array([ -r1*X[0] ])

x0 = array([dose])
t = linspace(0, 15, 256)
x = integrate.odeint(model1, x0, t)

figure(1, figsize=(6,7))

subplot(2,1,1)
plot(t, x[:,0], 'b-', linewidth=3)
ylabel('Plasma content (mg)', fontsize=18)

subplot(2,1,2)
plot(t, x[:,0]*mmg_per_mg/V_blood, 'g-', linewidth=3, label='Model')
ylabel('Plasma concentration ($\mu$g/ml)', fontsize=18)
xlabel('Time (hours from initial dose)', fontsize=18)
tight_layout()
savefig('one_compartment_timeseries.png',transparent=True)
savefig('one_compartment_timeseries.pdf') # vector graphics formats are better!
Predictions of our one-compartment model

Predictions of our one-compartment model

Acetaminophen kinetics

To evaluate our model, we need something to which we can compare our predictions. Below is a data set and some plots showing average observed acetaminophen kinetics in blood plasma.

[ Data : hide , shown as table , shown as CSV shown as QR ]

# observations of Acetaminophen concentration in the blood
# over 8 hours, following a 1000 mg dose (2 extra strength 500mg caplets)
# http://images.rxlist.com/images/rxlist/apap2.gif
#
# time (hours), plasma concentration (mmg/ml)
#
0.001, 0.001
0.170, 2.965
0.292, 11.65
0.463, 14.75
0.682, 15.10
0.829, 13.03
1.000, 11.51
1.146, 10.82
1.341, 10.27
1.487, 9.655
2.000, 8.00
4.000, 4.55
6.000, 2.34
8.000, 1.31
Plot of average plasma concentration after 1,000 mg dose of 2 caplets over 8 hours, on an arithmetic scale (left), and a log scale (right)

Plot of average plasma concentration after 1,000 mg dose of 2 caplets over 8 hours, on an arithmetic scale (left), and a log scale (right)

Before we start reading in too much, we should remind ourselves to be cautious. This data is a set of "average" observations. Every patient will be different in some way, and we have no way of telling how different yet. From observations of many individuals, we can always calculate averages, but we can not infer distributions from averages alone. When only averages are reported, information has been lost.

Still, this data set indicates some important similarities and differences between our simple model and the observed pharmacokinetics. First thing we should check is the scales. In the data, the largest concentration we observe is around 15 micrograms per milliliter. In our model, even starting with 1/3rd of the dosage, we predict concentrations as high as 65 micrograms per milliliter. Much less acetaminophen is getting into the blood than this model is predicting. From 2 to 8 hours, we see an exponential decay (indicated by the roughly straight line on the right side of the right plot). That is in good agreement with our simple model. But over the first 2 hours, our model predictions do not match the observations. In the patient, the plasma concentration increases steadily to a peak at around 40 minutes, and then a small rapid drop before relaxing to the exponential decay we conjectured.

An improved model

What are we missing that can explain the 40 minute lag between the time of ingestion and the time of peek plasma concentration? One hypothesis is that the patients ingested their doses, those did not immediate enter the blood. Instead, they entered the stomach, and it took time for the drug to be absorbed from the stomach into the blood stream. There may also have been some break-down of the acetaminophen that could explain the lower observed concentrations than predicted.

The natural way to improve our model, then, is to use two compartments instead of one compartment. The compartmental model then takes the form \[dx _ k/dt = - \text{Removal} + \text{Input} \pm \text{Movement}\] for the amount \(x _ k\) in compartment \(k\). Let's let \(x _ 0(t)\) represent the amount of acetaminophen in the gut, and let \(x _ 1(t)\) represent the amount in the blood stream. In both compartments, there will be elimination either through decay, metabolization, or excretion. Let's say that the elimination rate from the gut is \(r_0\) and the elimination rate from the blood is \(r_1\). Depending on how drug is administered, there could be inputs to the gut or the blood stream, so let's call those \(i_0(t)\) and \(i_1(t)\) respectively.

Now, we need to deal with movement from the gut to the blood stream. The same thought experiment that we used to motivate our terms for elimination apply to movement out as well -- each molecule should move independently and so, the net rate of movement out should be proportional to the mass of drug present. Thus, the net movement rate can be written be \(m _ {10} x _ 0\), where \(m _ {01}\) is a rate parameter for movement from the gut (compartment 0) to the blood stream (compartment 1), and has units of reciprocal time.

The difference between elimination and movement for our model is that movement conserves mass. So whatever moves out of the gut should reappear in the blood stream. This leads to the system of equations \[\begin{align*} \frac{dx_0}{dt} &= - r_0 x_0 - m_{10} x_0 + i_0(t), \\ \frac{dx_1}{dt} &= - r_1 x_1 + m_{10} x_0 + i_1(t). \end{align*}\] This model can be summarized with the following flow diagram.

Considering the case of a single oral dose like the data above, we can take \(x _ 0(0) = 1,000\), \(x _ 1(0) = 0\), and \(i _ 0(t) = i _ 1(t) = 0\), while keeping \(r _ 1 = 0.4\). We would like to estimate values for \(r _ 0\) and \(m _ {01}\) to match the peak concentration if we can. If you play around with parameter values (there is a better way to do this, which we will get to soon), you may find that if \(r _ 0 = 3\) and \(m _ {01} = 0.3\), we obtain a qualitatively reasonable fit, matching the timing and height of the peak concentration. This corresponds to a fraction \(r _ 0 /( m _ {01} + r_ 0) \approx 0.91\) of each dose being lost before it reaches the blood stream.

[Show code]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#!/usr/bin/env python3

"""
A two-compartment model for a single dose of acetaminophen
"""

from numpy import *
from scipy import integrate
from matplotlib.pyplot import *

rcParams.update({'font.size':12,'legend.framealpha':0.})


datafile='tylenol.csv'
data = loadtxt(datafile,delimiter=',')

r0 = 3.0
m01 = 0.3
r1 = 0.4
dose = 1000. # milligrams, mg
V_blood = 5000. # milliliters, ml
mgtommg = 1000.


def f(X, t):
    return array([ \
        -X[0]*(r0+m01), \
        X[0]*m01 - X[1]*r1 \
        ])

x0 = array([dose, 0.])
t = arange(0, 8.25, 0.1)
x = integrate.odeint(f, x0, t)

subplot(2,1,1)
plot(t, x[:,0], 'b-', linewidth=3)
ylabel('Gut (mg)', fontsize=14)
ylim(0,dose*1.1)
xlim(0,t[-1])

subplot(2,1,2)
plot(t, x[:,1]*(mgtommg/V_blood), 'g-', linewidth=3, label='Prediction')
plot(data[:,0], data[:,1], 'ro', markersize=4, label='Data')
ylabel('Plasma ($\mu$g/ml)', fontsize=14)
xlabel('Time (hours from initial dose)', fontsize=14)
xlim(0,8)
legend()

tight_layout()
savefig('two_compartment_timeseries.png', transparent=True)
savefig('two_compartment_timeseries.pdf') # vector graphics formats are better!
Our 2-compartment model's numerical prediction, together with the data

Our 2-compartment model's numerical prediction, together with the data

To help prevent overdose deaths and avoid potential FDA regulation, the producers of acetaminophen dropped the daily maximum number of dosage of 500mg acetaminophen tablets from 8 per day to 6 per day, e.g. from 4,000mg per day to 3,000mg per day. We can use our model to make predictions about the maximum plasm concentrations.

[Show code]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#!/usr/bin/env python3

"""
A two-compartment model for multiple doses of acetaminophen
"""

from numpy import *
from scipy import integrate
from matplotlib.pyplot import *

r0 = 3.0
m01 = 0.3
r1 = 0.4
dose = 500. # milligrams, mg
V_blood = 5000. # milliliters, ml
mgtommg = 1000.

dosepts = list(range(0,12,2))+[25]

dt = 0.05

def f(X, t):
    return array([ \
        -X[0]*(r0+m01), \
        X[0]*m01 - X[1]*r1 \
        ])

x0 = zeros(2)
xs = []
ts = []
i = 0
while i < len(dosepts)-1:
    x0 = x0 + array([dose,0])
    t = linspace(dosepts[i], dosepts[i+1], \
            int((dosepts[i+1]-dosepts[i])/dt)+1)
    ts.append( t )
    xs.append( integrate.odeint(f, x0, t) )
    x0 = xs[-1][-1,:]
    i += 1
t = concatenate(ts)
x = vstack(xs)


subplot(2,1,1)
plot(t, x[:,0], 'b-', linewidth=3)
ylabel('Gut (mg)', fontsize=14)
ylim(0,dose*1.1)
xlim(0,t[-1])

subplot(2,1,2)
plot(t, x[:,1]*(mgtommg/V_blood), 'g-', linewidth=3)
ylabel('Plasma ($\mu$g/ml)', fontsize=14)
xlabel('Time (hours from initial dose)', fontsize=14)
xlim(0,t[-1])


tight_layout()
fname = 'mdose_2box'
savefig(fname+'_timeseries.png', transparent=True)
savefig(fname+'_timeseries.pdf') # vector graphics formats are better!
show()
Predictions of plasma concentration when 6 doses of 500 mg are taken over 12 hours.

Predictions of plasma concentration when 6 doses of 500 mg are taken over 12 hours.

Discussion

Our two-compartment model's predicts visibly differ from the observations. The peak concentration is a little larger than predicted, and the long-term clearance is a little slower than predicted. Part of this difference may be fixed by a better choice of parameter values. But the errors in our predictions may also be because we are only using two compartments. The actual physiology is more complicated. Absorption rates can be different in different parts of the gut. Acetaminophen in the gut must pass through the liver before it gets to the blood stream. Once it the blood stream, it can get absorbed and re-released by variety of tissues and organs. Our model could be expanded to account for many of these factors.

A flow diagram for the body capturing the major components of drug pharmacokinetics.

A flow diagram for the body capturing the major components of drug pharmacokinetics.

These possibilities can be represented by the general vector equation \[\frac{dx}{dt} = I(t) + M x - E x,\] where \(x(t)\) is the vector of compartment amounts, \(I(t)\) is a vector of input rates, \(E\) is a diagonal matrix of elimination rates, and \(M\) is a movement matrix with non-negative off-diagonal entries and columns that sum to zero (so that mass is conserved by movement). These models are a special case of the general class of linear inhomogeneous systems, so the theory associated with their solution (eigenvalues, matrix exponentials, Green functions) applies.

The pharmacokinetic models we've developed here so-far are quite preliminary. While they do a descent job of recovering the observed average dynamics of acetaminophen in the blood stream, we're still far from a validated model that might be used for designing policy. There are a number of physiological things we have not yet addressed, such as circadian rhythms, metabolization, and toxicity. Even more importantly, there can be large variations in response to a particular dose between patients, based on weight, gender, race, age, and other factors, and when designing a public-policy recommendation, it's important to design recommendations that will be safe for close to 100% of the population, not 50% or 90% of people. It is commonly the case that constructing a highly accurate and reliable model can involve 10 or 100 times as much labor as making a crude model of common circumstances because of the extensive data collection and validation needed.

Exercises

  1. If the function \(y(x)\) is a function such that for every positive value of \(a\), \(y(a x) - a^2 y(x) = 0\), what is \(y(x)\)?

  2. Using our simple explain why the formula half-life formula \(t _ {1/2} = \ln 2 / r\) correctly relates the elimination rate \(r\) and the halflife \(t _ {1/2}\).

  3. Imagine a compartmental model consisting of 3 compartments. Drug can between any of 3 compartments, but and only be input into the first. Draw a flow diagram for this model.

  4. Suppose a hospital patient is given a 325 mg dose of acetaminophen orally, and then two hours later is put on an intravenous drip of 10 mg per hour directly into the blood stream. Use our 2-compartment model to predict the acetaminophen plasma concentration over the 12 hours since the oral dose.

  5. (to complete) Alternatively, all the drug could be delivered in a single does at the 1-hour mark, in which case \[dx/dt = a \delta(t-1) - rx,\] where \(a\) is the size of the drug dose and \(\delta()\) represents Dirac's delta-function. Dirac's delta-function is a special function that integrates to one, yet is zero everywhere but the immediate neighborhood of \(0\). It sometimes helps to think of Dirac's delta-function as what you get when you squash a normal distribution's variance down to zero -- a full unit of area concentrated at the origin.

  6. Our single-dose two-compartment model given above can be solved in closed form using the standard theory of linear systems. Find a closed-form equation for \(x_1(t)\) in closed-form.

  7. Show by differentiation that the following equation solves the general linear pharmacokinetics system.
    \[x(t) = e^{(M-E)t} x(0) + \int _ {0} ^t e^{(M-E)(t-s)} I(s) ds\]

  8. include a second application illustrating environmental contamination modelling

    The general ideas of pharmacokinetic modeling can also be applied to the modelling of environmental contamination.

( previous, home, next )