( previous, home, next )

include a second application to evironmental contamination modelling


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 300 deaths from acetaminophen poisoning per year 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 and predict safe drug usage. and we will only have space to touch on the the basics.

Linear rate laws

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 = - R\] where \(R\) is the removal rate. So, what formulas do we fill in for the removal?

Drugs are taken out of a compartment in the body by various processes like excretion and metabolism. 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 eliminated from the body should be independent of the other drug molecules. If each molecule’s elimination risk 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\) should satisfy \[R\left(\frac{x}{2}\right) = \frac{1}{2}R(x).\] In fact, this should hold for any constant \(a\): \[R(a x) = a R(x).\]

We now have an equation for \(R\). 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 R' = 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 elimination rate depends on the chemistry of the drug and physiology of the compartment, and is often estimated based on observational data.

A first compartmental model

To better understand the dangers of acetaminophen and the FDA’s recommendations, it helps to build a model of the drugs kinetics in a human body.

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 = I - R\] where \(I\) is the input rate and \(R\) is the output rate. So, what formulas do we fill in for the inputs and outputs?

The inputs are controlled externally by the dosing regiment, and so, can be treated as a function of time. Thus, the general form of our single-compartment model is \[dx/dt = i(t) - r x \quad \text{with $x(0)$ given.}\] This is summarized by the following diagram

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 mass into concentrations. As a crude estimate, we’ll assume our dose is administered to an adult with 5 liters of blood.

[Show code]
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))

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

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)
savefig('one_compartment_timeseries.pdf') # vector graphics formats are better!
Predictions of our initial model
Predictions of our initial model

To evaluate our model, we need something to which we can compare our predictions. Fortunately, there is some data available. Below are a pair of plots showing some observed acetaminophen kinetics in the blood.

[ 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)
# 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 does 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 does of 2 caplets over 8 hours, on an arithmetic scale (left), and a log scale (right)

First, before we start reading in too much, we should to reminder ourselves to be cautions. This set 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 individiuals, we can always calculate averages, but not vice-versa, so when averages only 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. So much less acetaminophen is getting into the blood than we predict. From 2 to 8 hourse, we see an exponential decay (indicated by the roughly straight line on the right side of the right plot). So that in surprisingly 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 intuitive explanation 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 in the stomach 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 general form of a compartmental model takes the form \[dx _ i/dt = \text{Input} \pm \text{Movement} - \text{Removal}\] for the amount \(x _ i\) in compartment \(i\). 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 as well – each molecule should move independently and so, the net rate of movement should be proportional to the mass of drug present. If \(m _ {01}\) is the rate of movement from the gut to the blood stream, then the net movement rate will be \[ m_{10} x_0.\] The difference between elimination and movement for our model is that movement conserves mass (though not concentration!). 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*}\]

Considering the case of a single oral dose like the data above, we can take \(x_1(0) = 0\) and \(i_0(t) = i_1(t) = 0\).

[Show code]

A two-compartment model for a single dose of acetaminophen

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

# 1.75 hour halflife
r1 = log(2.)/1.75
r0 = log(2.)/4.75
m01 = 2.77 # if r0=log(2.)/1.75
#m01 = 3.17 # if r0=0
dose = 500.*2
V_blood = 5000. # 5 liters approximately

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

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

plot(t, x[:,0], 'b-', linewidth=3)
ylabel('Stomach content (mg)')

plot(t, x[:,1]*1000./V_blood, 'g-', linewidth=3)
ylabel('Plasma concentration ($\mu$g/ml)')
xlabel('Time (hours from initial dose)')

savefig('two_compartment_timeseries.png', transparent=True)
savefig('two_compartment_timeseries.pdf') # vector graphics formats are better!
Our 2-compartment model’s numerical solution
Our 2-compartment model’s numerical solution

More information

\[dx/dt = i(t) + M x - Ex\] where \(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 one so that mass is conserved by movement.

A full body flow network model

Include discussion of effects of compartmental volume and adjusting flow rates to take volume differences into account.


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 excretion by the kidneys, tissue absorbtion, and toxicity effects. 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.

To help reduce these deaths, the Food and Drug Administration (FDA) took steps in 2011. The recommended daily maximum number of doses of 500mg acetaminophen tablets were dropped from 8 per day to 6 per day, i.e. from 4,000mg per day to 3,000mg per day. The FDA also asked manufacturers of prescription combination drugs to limit the amount of acetaminophen to 325 mg per tablet or capsule by January 14, 2014. In March 2014, the FDA announced all manufacturers had stopped marketing products with more than 325 mg of acetaminophen. We’re still waiting to learn the outcomes of these changes.


  1. Imagine a compartmental model consisting of 3 compartments. Drug can between any of

  2. (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.

  3. 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\]

( previous, home, next )