- Antoine Lavoisier’s discovery of the Law of Conservation of Mass (late 1700’s)
- Acetaminophen information at wiki, information, news, news.
- Clements, 1978 is a paper using linear ODEs to model clinical observations of acetaminophen absorption in 8 men.
- The History of Pharmacokinetics by Wagner, 1977.
- Use only as directed story by
*This american life*. - FDA page on Acetaminophen

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, and around 1.5 billion doses were sold in the United states in 2018. US CDC estimated there were 1,600 deaths from acetaminophen poisoning between 2001 and 2010. That’s a relatively small number given its heavy usage. For comparison, about 33,000 people died in motor vehicle accidents in 2010 in the US. But 160 per year is still a significant number of preventable deaths. 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 help prevent these deaths?

This is a standard question in the field of pharmacokinetics with which we can introduce compartmental modelling using first-order differential equation. But be aware that pharmacokinetics and the related area 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.

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{Input} - \text{Removal}.\]

The input rate is something that we can determine based on observation, be it a constant rate because it is being given to a hospital patient intravenously, or pulses associated with oral doses. Generically, we can represent all of these with a function \(i(t)\). Drugs are removed from the body by various processes like excretion and metabolization in ways that are not immediately observable. 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\) that 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 dimensions 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.

Using these hypotheses regarding input and removal, our single-compartment model can be written as \[dx/dt = i(t) - r x \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 inputs can be specified in terms of the initial condition: \[\frac{dx}{dt} = - r x, \quad x(0)\;\;\text{given}.\]

To turn these ideas into something practical, we need estimates the initial value \(x(0)\) and the elimination rate \(r\). A little research finds the following. An extra-strength dose of acetaminophen is 1,000 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]

```
"""
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 = 1000.
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!
```

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
```

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 above, we predict concentrations as high as 200 micrograms per milliliter. Much less acetaminophen is getting into the blood than our model is predicting.

From 2 to 8 hours, the data shows an exponential decay in acetaminophen concentration (indicated by the roughly straight line on the right side of the right plot). That is in good agreement with our 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.

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 _ {10}\) 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.We can calculate the solution in closed-form using standard methods and show when there is no on-going input (\(i_0(t) = i_1(t) = 0\)), \[\begin{align} x_0(t) &= x_0(0) e^{-(r_0 + m_{10})t}, \\ x_1(t) &= \frac{x_0(0) m_{10}}{r_1 - r_0 - m_{10}} \left( e^{-(r_0 + m_{10})t} - e^{-r_1 t} \right). \end{align}\]

Considering the case of a single oral dose like the data above, we can take \(x _ 0(0) = 10^3\), \(x _ 1(0) = 0\), while keeping \(r _ 1 = 0.4\). We would like to estimate values for \(r _ 0\) and \(m _ {10}\) 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 _ {10} = 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 _ {10} + r_ 0) \approx 0.91\) of each dose being lost before it reaches the blood stream.

[Show code]

```
#!/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
m10 = 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+m10), \
X[0]*m10 - 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!
```

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 in 2011. We can use our model to make predictions about the maximum plasm concentrations under this dosing regime.

There is one catch we have to overcome in proceding – we can no longer use the initial condition to specify our full dosing regime.

It is situations like this for which we include the input terms. In our current scenario, there’s no direct input to the blood, so \(i_1(t)=0\), but for the gut we have to specify \(i_0(t)\) somehow. We know the patient will take one 500 mg dose every 2 hours. But the dimensions of \(i_0(t)\) are amount per time, *not* amount! What’s going on?

The problem is that the consumption of a pill imposes an “impulsive” forcing on the system – there is a very large jump in the gut’s acetaminophen level in an instant. A jump so large that the gut’s state variable \(x_0(t)\) isn’t continuous at the time the dose is ingested! Differential equations were initially intended for use in the study of systems where states only change a very little in each instant, so they can be treated as continuous, smooth functions which have a well-defined derivative. Most of the time, the amount of acetaminophen changes very little in each instant, *but* there are a few instants when it changes allot and these changes are important.

Fortunately, there is a way to handle this issue – an impulse or delta function \(\delta(x)\). An impulse function isn’t really a function – it’s sort-of the reciprocal concept to an infinitessimal. Infinitessimals are numbers so small that they are indistinguishable from zero on their own, but can ammount to something when you add just infinitely many enough of them. An impulse function is an idealized distribution that is zero everywhere except at the origin where it is infinitely large, but it is infinitely large in just the right way so that if you integrate the impulse funciton over any open interval containing the origin, you get exactly \(1\). Thus, it is like a normal distribution with no variance. Impulse functions can be shifted just like regular functions, so if we want to define an inflow function that inputs 6 doses of 500 mg of acetaminophen every 2 hours, we can use \[i_0(t) = 500 \sum _ {k=0}^{5} \delta(t - 2 k).\]

While impulse functions can be very useful conceptually and numerically, they are a bit of a challenge to deal with numerically – because they are so narrow and our numerical methods always use finite time steps, algorithms tend to skip right over them without noticing they were there. There are a number of ways we can fix this. One of the simplest is to replace the impulse function’s infinitessimal width with a finite width and force our differential equation integrator to use small enough step sizes to resolve these impulses. This approach tends hurt algorithm accurracy a little, but for our purposes of getting only a roughly accurate solution, this is acceptable. It all leads to the solution pictured below, and indicates that the blood’s acetaminophen concentration will max out around 13 \(\mu\) g/ml.

[Show code]

```
#!/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.
dt = 0.05
def delta(x,w=1/60):
# one-sided impulse after 0
return max(0,w-abs(x-w))/w**2
# two-sided impulse
#return max(0,w-abs(x))/w**2
def I(t):
ivec = array([0.,0.])
for a in range(0,11,2):
ivec[0] += dose*delta(t-a)
return ivec
M = array([[-m01,0],[m01,0]])
E = diag([r0,r1])
def f(X, t):
return I(t) + (M-E).dot(X)
tcrit = array(range(0,11,2))+0.5/60
t = linspace(0,25,4000)
x0 = zeros(2)
x = integrate.odeint(f, x0, t, tcrit=tcrit)
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()
```

Our two-compartment model’s predictions 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.

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.

If \(Y(x)\) is a continuous function solving the functional equation \(Y(a x) - a^2 Y(x) = 0\) for every positive value of \(a\) and \(x\), what is \(Y(x)\)?

Use our differential equation \(\dot{x} = r x\) to 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}\).

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.

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.

- Suppose the peak blood concentration of acetaminophen is observed to occur at \(t^*:=0.75\) hours and reaches a maximum of \(15\) \(\mu\)g/ml.
- Differentiate our formula for \(x_1(t)\) given above to show that \(r_0 + m_{10}\) solves the transcendental equation \[(r_0 + m_{10}) e^{-(r_0 + m_{10}) t^*} = r_1 e^{- r_1 t^*}.\]
- Find a numerical estimate for \(r_0 + m_{10}\).
- Assuming a person has \(V = 5,000\) milliters of blood, use the equation \(x_1(t^*) = 15 V\) to estimate \(m_{10}\).
- Estimate \(r_0\).

(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 squeeze a normal distribution’s variance to zero – a full unit of area concentrated at a single point.

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\]- A simple compartmental model is sometimes used to describe unemployment rates in economics. Let \(U\) be the fraction of the population unemployeed, and \(E\) be the fraction of the population employed, such that \(E+U = 1\). \[E \xrightarrow{f} U , \quad U \xrightarrow{s} E.\]
- What is the corresponding differential equation system?
- Find \(U(t)\), assuming \(U(0)\) is given.

One way we design better drugs is to attach polyetholene glycal to the molecule to slow it’s removal from the blood by the kidneys. This was a major scientific advance in the greatment of HIV and viral hepatitis. Modify the lecture example code for our two-compartment acetaminophen model to simulate an initial dose of 600 mg and 5 subsequent doses at 300 mg every four hours when the removal rate from the blood \(r_1 = 0.2\). Plot the amount in the gut and the concentration in blood plasma over 36 hours.

- Consider the following graphical model representation below, and answer these questions about it.
- What are the state variables?
- What are the rate parameters?
- Write this model out as a system of reactions.
- Write this model out as a system of ordinary differential equations.

- Little Blue Run is a coal ash pond in Pennsylvania that is leaking selenium into nearby water ways. Make a linear compartmental model that can describe the accumulation of selenium in a nearby lake’s, bottom sediments, and a fish living in the lake.
- Decide on a set of units to be used in your model.
- Clearly define a set of state variables for your model and their units.
- Identify a set of input, removal, and movement reactions for the selenium in the system, and define rate parameters for each reaction.
- Write down a system of ordinary differential equations for your model.
- Draw and label a graph of your model.