A Spoonful of Medicine

Makes the Mathematics Go Down

Sheldon P. Gordon and Florence S. Gordon

Dr. Sheldon Gordon is Professor of Mathematics at SUNY Farmingdale. He is a member of a number of national committees involved in undergraduate mathematics education. He is the principal author of Functioning in the Real World and a co-author of the texts developed under the Harvard Calculus Consortium.

gordonsp@farmingdale.edu

________________________________________
Dr. Florence S. Gordon is Professor of Mathematics at New York Institute of Technology. She is a co-
author of Functioning in the Real World, co-author of Contemporary Statistics: A Computer Approach and co-editor of the MAA volume, Statistics for the Twenty First Century.

fgordon@nyit.edu

Introduction

For most of the last century, mathematics and the various scientific disciplines have become increasingly isolated from one another. This isolation has been reflected in our curricula, where the applications taught in mathematics courses tended to be highly artificial, and so those courses became increasingly less relevant to the other disciplines. In turn, the other disciplines began to teach more and more of the mathematics they think they need.

In the last few years, things have started to change for a variety of reasons. Traditional discipline lines are becoming increasingly meaningless as current research and developments occur in fields that overlap two or more disciplines. This is reflected strongly in the scientific and engineering workplace, which is increasingly interdisciplinary in nature. The reform movement in mathematics has, among other things, placed a renewed emphasis on realistic applications to motivate the mathematical developments. Finally, modern technology now allows us, trivially, to do things in the mathematics classroom that previously would have required a roomful of what not-so-long ago were state of the art computers.

The key to introducing such approaches into the mathematics curriculum is to focus on mathematical modeling as a unifying theme. Thus, instead of handing the students equations to solve or functions to differentiate or integrate, we would begin by presenting interesting situations that require the students to:

•	identify the mathematical component,

•	create an appropriate mathematical model (often a function or equation),

•	raise pertinent questions in the context of the situation,

•	use the mathematical model to answer those questions,

•	interpret the solutions to see if they make sense.

The mathematical techniques needed may arise and be developed in the context of the problem. Alternatively, these methods are often ones that the students likely have seen previously, though often have not mastered (probably because they saw little value in practicing skills without context or meaningful applications).

In the present article, we illustrate how one such topic, modeling the level of a drug or medication in the bloodstream, can be introduced and modeled at a variety of different levels throughout the curriculum from developmental mathematics through college algebra and precalculus, and then again at the calculus level. At all levels of the curriculum, this application demonstrates dramatically the usefulness of the mathematics the students are learning. In addition, it provides an excellent vehicle for introducing some exceptionally nice mathematical ideas and methods. To paraphrase Mary Poppins, a spoonful of medicine makes the mathematics go down, in a most delightful way.

Modeling Drug Responses

L-dopa is administered to people suffering from Parkinson's disease to relieve symptoms such as extreme tremors and rigidity. Consider the following set of data on the levels of L-dopa in the blood, in nanograms per milliliter, as a function of time, in minutes after the drug is administered.

t	0	20	40	60	80	100	120	140	160	180	200	220	240	300	360
L	0	300	2700	2950	2600	1550	1100	900	725	600	510	440	300	250	225
Table 1

A plot of these points is shown in Figure 1. (We note that comparable graphs and consequently the corresponding data values for many other drugs can be found at the websites for the major pharmaceutical companies.)

The pattern shown in this drug response curve is typical of what happens with any medication. The underlying assumption is that there is no medication in the bloodstream when the drug is first taken, that the level rises rapidly as the drug is absorbed into the blood, and the level then decreases as the drug is washed out of the blood by the kidneys.

We first develop models, in several different ways, for the decay in the drug level as a medication is washed out of the blood. We will come back to develop models for the overall pattern in the drug response curve later. We note that different drugs are "washed out" of the blood at different rates. In fact, every drug is rated by how long it takes for 50% to be eliminated from the body and this is known as its biological half-life or simply the half-life. For example, aspirin is removed quite rapidly; its level is reduced by about 50% every 29 minutes.

We first construct a mathematical model for the level of aspirin in the bloodstream once it has been fully absorbed into the blood. Based on the half-life of 29 minutes, in any one-hour period, approximately 76.2% of the aspirin in the blood is eliminated, leaving about 23.8%. If the initial dosage of aspirin is 500 mg, we can model the level of aspirin in the blood with the exponential decay function

A(t) = 500 (0.238)t ,

where t measures the number of hours from when the aspirin was fully absorbed into the blood. Alternatively, if t measures the number of minutes from when the aspirin was fully absorbed, we can model the aspirin level by the exponential function

A(t) = 500 (0.9764)t .

The base or decay factor, 0.9764, indicates that almost 2.5% of the aspirin is eliminated every minute.

Once such an exponential model has been created, there are two kinds of questions that arise naturally:

(a) How much aspirin is left after any given length of time?

(b) How long does it take until the level of aspirin is down to a given amount?

Both types of questions make sense to students and they see the value in developing the model as a means of answering such questions. And, because the questions make sense, they are much more likely to expend the time and effort solving the corresponding algebraic equations than they are to solve a long list of context-free equations just for the sake of practice.

Alternatively, suppose we start with a collection of data on actual drug levels instead of being given a half-life. Consider again the table of data values on the level of L-dopa in the blood. We now focus on the data points while the level of L-dopa in the blood decays. We then have the following table where we re-label the starting time to be t = 0 when the medication is apparently at its peak level.

t	0	20	40	60	80	100	120	140	160	180	240	300
L	2950	2600	1550	1100	900	725	600	510	440	300	250	225
Table 2

We plot the data in Figure 2, which suggests that a decaying exponential function is a reasonable model. (A decaying power function would not be reasonable because there clearly is no vertical asymptote at the origin.)

The corresponding exponential regression equation is

L(t) = 2147.9(0.9909)t,

where t measures time in minutes from absorption. This model tells us that the level of L-dopa in the blood decreases by almost 1% every minute. We show this function superimposed over the data in Figure 3 and see that it is a reasonably good, though not an outstanding, fit. In addition, we note that the correlation coefficient r = -0.9545, which also suggests that the fit is good, especially as t increases.

Modeling the Drug Response Curve

Now let's consider the problem of modeling the entire drug response curve rather than just the portion of it where the level of medication is decreasing. Conventional wisdom [3] suggests that such curves are modeled by surge functions of the form

S(t) = A t p e-kt or	S(t) = A tp bt,

where A, p, and k are constants. The graph of
such a function is shown in Figure 4; we see
that it clearly exhibits the behavior pattern
displayed in the data on L-dopa in Figure 1.
We thus seek the equation of a surge function
that models this data.

From the plot of all the L-dopa data in
Figure 1, we see that the surge function should
reach its maximum at about t = 60, where the maximum value is approximately 3000. We also see that the surge function has two points of inflection, one on either side of the peak. From the data in Table 1, notice that the greatest increase in S occurs between t = 20 and t = 40, so we estimate that one inflection point occurs at about t = 30. The greatest decrease in L occurs between t = 80 and t = 100, so we estimate that the other inflection point occurs at about t = 90.

The general equation of a surge function, S(t) = A tp e-kt, involves three parameters, A, p, and k, whose values we have to determine. We differentiate this function with respect to t to obtain

S'(t) = A[ p t p-1 e-kt - k t p e-kt]

= A[ p - k t] t p-1 e-kt .	(1)

We set this expression to 0 and solve. Since the exponential function is never 0 and t p-1 is zero only when t = 0, the maximum occurs when

p - k t = 0 or t = p/k .

From our previous examination of the data and the corresponding plot, we decided that this occurs at about

t = p/k = 60 or p = 60k
(2)

We next find the points of inflection. When we differentiate Equation (1) for S '(t) above, we get

S "(t) = A{[ p (p - 1) t p-2 - k p t p-1] e-kt - k [ p t p-1 - k t p] e-kt}

= A{[ p (p - 1) t p-2 - k p t p-1]- k [ p t p-1 - k t p]}e-kt .

After some algebraic simplification, this reduces to

S "(t) = A{ p (p - 1) - 2k p t + k2 t2 }t p-2 e-kt .

The inflection points occur when the term inside the brackets is zero:

p (p - 1) - 2k p t + k2 t2 = 0.

This is a quadratic equation in t whose roots are

This tells us that the two inflection points are located at equal distances of ( p)/k on either side of the maximum point at t = p/k . Since we decided that the points of inflection occur at approximately t = 30 and t = 90, we must have

p /k = 30 or p = 30k.	(3)

Using Equation (2), we find that

p = 900k2 =	60k

so that

k = 1/15 and p = 60k = 4.

Finally, to estimate the value for A, we use the expression for the surge function

S(t) = A tp e-kt = A t4 e-t/15,

and the fact that S(60) = 3000 to get

S(60) = A 604 e-60/15 = 3000,

so that

A = 3000 e60/15 / 604 0.0126.

Therefore, our model for the surge function is

S(t) = 0.0126 t4 e-t/15 .

The graph of this function superimposed over
the original data set for the level of L-dopa in
the blood is shown in Figure 5. We note that
it is a reasonably good fit for t between 0 and
about 120 minutes, although thereafter the
function decreases much more rapidly than the
level of the drug does.

We can improve on the fit somewhat
by applying the least squares criterion directly to minimize the sum of the squares of the deviations between the data values and the surge function rather than depending on inspection to estimate the location of the maximum and the inflection points. In particular, the sum of the squares corresponding to the function shown in Figure 5 is 1,716,480; the sum of the squares is reduced by about 10% to 1,555,393 using the similar function

S(t) = 0.012764 t3.962 e-0.0637 t .

Nevertheless, this is still a rather poor fit to the data, especially beyond about t = 140 minutes.

Using a Rational Function

An alternative approach to creating a function to fit this type of data is discussed by Brown and Timchek [1], who suggest using a rational function. At first thought, this likely seems unreasonable, since there are no vertical asymptotes. However, further thought might suggest that a rational function having no linear factors in the denominator could have this property. The fact that the curve passes through the origin might then suggest a rational function of the form

which dies out as x ∞. However, such a function does not have the appropriate behavior near the origin; we need a curve whose slope at 0 is 0 and which is concave up for x > 0, but the second derivative of this rational function is negative for 0 < x < b. After a bit of experimentation, it turns out that a rational function having the desired behavior pattern is

Brown and Timchek indicate that, using the curve fitting routines in Maple, the least-
squares criterion applied to such a function leads to the specific function

We show the graph of this function super
imposed over the L-dopa data in Figure 6, where we observe that it is actually quite accurate. The corresponding value for the sum of the squares is 620517.5, which is a significant improvement over the values obtained using the surge function.

Repeated Drug Doses

In the above developments, we considered the situation where a single dose of a medication is taken. In many cases, people take a dose of a medication on a repeated, often daily, basis. It might be insulin for diabetes, any of a variety of medications to lower high blood pressure, or drugs to reduce cholesterol. In each instance, once the medication has been absorbed into the blood, it is filtered out by the kidneys and excreted from the body. However, before it is totally eliminated, the next dose comes along. The problem then is: How do we model the level of the medication in the blood as a function of time?

To study this type of situation, suppose we consider Prozac, which is administered to individuals to counter the effects of extreme depression. A typical dose of Prozac is 40 mg and approximately 25% of the Prozac in the bloodstream is eliminated every 24 hours. If just a single dosage of Prozac were taken, the level in the bloodstream could be modeled by the exponential decay function

D(t) = 40(0.75) t,

where t measures the number of 24 hour periods since the Prozac was first taken.

More realistically, suppose that a person takes 40 mg of Prozac each day. After the first 24 hour time period, 25% of the Prozac, or 10 mg, is eliminated, leaving 30 mg and the next day's dose adds 40 mg to that. Therefore, after the first 24 hours, the amount of Prozac in the blood is

D1 = 30 + 40 = 70 mg.

During the second 24 hour period, the kidneys eliminate 25% of the Prozac present and the patient takes the next day's dose of 40 mg. Thus, after two 24 hour periods, the amount of Prozac in the blood is

D2 = 0.75(70) + 40 = 92.5 mg.

After three days, the level of Prozac is

D3 = 0.75(92.5) + 40 = 109.375 mg,

and so on, indefinitely. The corresponding sequence of Prozac levels is therefore

{40, 70, 92.5, 109.375, 122.031, 131.523, ... }.

Notice how the levels of the drug keep rising, but in a concave down manner.

Let's now look at a slightly more sophisticated way to describe this process. The initial dosage is D0 = 40 mg. During the first 24 hour time period, 25% of this amount is removed from the blood and the patient then takes the next dose of 40 mg. Thus

D1 = 0.75 D0 + 40.

Similarly, during the second day, 25% of the Prozac is eliminated and the person takes another 40 mg dose, so that

D2 = 0.75 D1 + 40

and, again, after three days,

D3 = 0.75 D2 + 40.

In general, at the end of n+1 days, for any value of n,

Dn+1 = 0.75 Dn + 40.

This equation, which shows the relationship between the level of Prozac on any two successive days, is called a difference equation or recurrence relation.

Since the initial dose is D0 = 40 mg, if we set n = 0, 1, 2, ... in this difference equation, we obtain, respectively,

If n = 0:	D1 = 0.75 D0 + 40	= 0.75(40) + 40 = 70 mg.

If n = 1:	D2 = 0.75 D1 + 40	= 0.75(70) + 40 = 92.5 mg.

If n = 2:	D3 = 0.75 D2 + 40	= 0.75(92.5) + 40 = 109.375 mg

and so forth, so that we generate the same sequence

{40 , 70, 92.5, 109.375, 122.031, 131.523, ... }

as above. This sequence is the solution to the difference equation.

Notice that each successive term has grown by somewhat less than the term before. The curve drawn through the points in Figure 7 is concave down and the successive values seem to be leveling off. If you continue the above process numerically, you will see that the limiting amount L appears to be very close to 160 mg.

This recursive pattern can be generated in a trivial manner on most graphing calculators (or with a spreadsheet). With a calculator, start by keying in the initial value, say 40, and pressing ENTER. Then enter the expression

.75 * 2nd ANS + 40

and, each time you press ENTER thereafter, the following term in the sequence is displayed.

In terms of the original situation, this limit L represents the maximum level of Prozac that will ever be reached in the blood. It is known as the maintenance level for the drug. Furthermore, once that level of Prozac has been reached, the amount in the blood will remain constant at that level every 24 hours, so long as the same dosage is taken repeatedly.

In practice, medical researchers determine that a specific level L of a medication is most effective, considering factors such as safety and effectiveness. An initial dose of 40 mg of the Prozac means that for some period of time, the amount in the bloodstream is below the optimal level. Because of this, doctors often prescribe an initial dose that is above the normal repeated dose so that the drug level approaches the maintenance level L more rapidly. For example, an initial dose of 120 mg of Prozac followed by daily doses of 40 mg will achieve the maintenance level much more rapidly. However, this strategy must be balanced with the safety issue involved in taking such a large dosage, especially as the first dose of the drug.

What happens if a patient takes an overdose, so that the level of drug in the system exceeds the maintenance level? For instance, suppose that a person takes 400 mg of Prozac initially and thereafter takes the usual 40 mg daily dose. We use the same difference equation

Dn+1 = 0.75 Dn + 40 ,

but now with an initial dose D0 = 400, so that

if n = 0: D1 = 0.75 D0 + 40 = 340 mg

if n = 1: D2 = 0.75 D1 + 40 = 295 mg

if n = 2: D3 = 0.75 D2 + 40 = 261.25 mg

and so on.

The levels of the drug in the bloodstream are now given by the solution sequence

{400, 340, 295, 261.25, 235,937, 216.953, 202.715, 192.036, 184.027, ... }.

We plot these points in Figure 8 and observe that they fall into a decreasing, concave up pattern that apparently converges to the same limiting value L = 160 mg, but they approach it from above rather than from below. Consequently, we see that if the drug level rises too high, there are some counteracting effects that reduce the level, assuming the overdose does not cause any other problem.

Notice that by changing the initial condition D0 in the difference equation, we obtain a different solution sequence. In fact, for every possible initial value D0, there will be a different solution.

Determining the Maintenance Level L

We can determine the limiting value L for Prozac precisely using the following argument. Suppose that for some value of n, Dn reaches the limit L, so that all successive levels of Prozac are the same. Thus for n large enough, we assume that Dn+1 = L and Dn = L. Substituting these values into the difference equation

Dn+1 = 0.75 Dn + 40 ,

for the Prozac drug model, we obtain

L = 0.75L + 40

so that

0.25L = 40

and hence the limiting value L is

L = 40 /0.25 = 160 mg.

More generally, this model can be used for any medication. Suppose that, for a certain drug, a fraction a is eliminated from the bloodstream during a given time period, so that the fraction remaining in the blood after that length of time is b= 1 - a. If the repeated dosage is C, then the corresponding difference equation is

Dn+1 = b Dn + C

and the associated maintenance level for the medication is

L = C/(1 - b) = C/a.

Finding a Formula for the Solution

We now attempt to find a closed-form expression for the solution sequence of the difference equation

Dn+1 = 0.75Dn + 40

for Prozac and the more general difference equation model

Dn+1 = b Dn + C.

Although there are some powerful techniques for finding formulas for the solution sequences to such difference equations (they mirror the methods used to find solutions of differential equations), we will not go into them here. The interested reader can find these methods in either [1] or [2]. Instead, we will construct such a solution using some fundamental ideas about functions.

Let's look again at the terms we calculated above for the solution sequence to the Prozac difference equation model based on an initial value of D0 = 40 mg:

{40 , 70, 92.5, 109.375, 122.03125, 131.523, ... } .

We can think of these values as a set of data values:

n	0	1	2	3	4	5
Dn	40	70	92.5	109.375	122.031	131.523
(We could use more values, but these will suffice.) Since the values for Dn approach a limiting value of L = 160 mg as a horizontal asymptote, we can shift each Dn value to obtain the corresponding values of 160 - Dn, as shown in the following table and in Figure 9.

n	0	1	2	3	4	5
Dn	40	70	92.5	109.375	122.031	131.523
160 -Dn	120	90	67.5	50.625	37.969	28.477

Since Dn approaches 160 as n increases, the values of 160 - Dn approach 0 as n gets larger. It therefore makes sense to fit a decaying exponential function to 160 - Dn as a function of n. (We would not use a decaying power function because the data starts with n = 0.) The exponential function that best fits this data using any of the standard calculators or Excel is

160 - Dn = 119.99961 (0.7500021) n,

with a correlation coefficient of r = -1, which suggests a virtually perfect fit. If we solve for the level of Prozac, Dn, we obtain

Dn = 160 - 119.99961 (0.7500021) n.

The numbers in this expression suggest that the "correct" formula for the solution might be

Dn = 160 - 120(0.75) n.	(4)

We will come back to this problem later to show that the formula truly holds for every possible value of n, not just for the few particular values of n we used in constructing the best-fit exponential function.

We now extend the above formula to solve the comparable difference equation for any medication with any fixed periodic dose. The corresponding difference equation is

Dn+1 = b Dn + C.

For any value of b between 0 and 1 and any positive value for C, the successive terms in the solution sequence for Dn will have the comparable behavior to that shown in Figure 7: the solution will be an increasing, concave down function that approaches a horizontal asymptote. If your calculator displays graphs of solutions of difference equations, select some typical values and check out the behavior of the solutions.

We found above that the maintenance level is given by

L = C/a = C/(1 - b).

If you examine the solution (4) we created for the level of Prozac in the blood, it is apparent that a more general formula for the solution would likely be

Dn = L - (L - D0) bn

for any value of n, with any parameters a, C, L = C/(1 - b), and D0. To verify that this expression for Dn is actually a formula for the solution, we must show that it satisfies the difference equation

Dn+1 = b Dn + C	(5)

for every value of n. (This is essentially an application of mathematical induction.) We substitute both the above assumed expression for Dn = L - (L - D0) bn and the corresponding expression for Dn+1 when n is replaced by n + 1 into the difference equation. The left-hand side of the difference equation (5) can be expressed as

Dn+1 = L - (L - D0) bn+1.

The right-hand side is

b Dn + C = bL - b(L - D0) bn + C

= bL - (L - D0) bn+1 + C.

However,

L = C/ (1 - b) so that C = L(1 - b).

Therefore, the right-hand side becomes

b Dn + C = bL - (L - D0) bn+1 + C

= bL - (L - D0) bn+1 + L(1 - b)

= L - (L - D0) bn+1,

which is identical to the assumed expression for Dn+1 on the left-hand side. Thus, the expression

Dn = L - (L - D0) bn

is truly a formula for the solution sequence Dn and it holds for all possible values of n.

This solution applies to any difference equation of the form

xn+1 = b xn + C,

not just the drug model discussed here. In particular, any situation that can be represented in a decay-replenish or growth-diminish model can be treated in the identical fashion. As several examples, consider removing a fixed amount of money from a retirement fund that grows at a given rate or harvesting a fixed number of animals from a herd that is growing at a given rate or contributing a fixed amount, say $2000, annually to an IRA account growing at a given rate.

References

1. Gordon Sheldon P., Florence S. Gordon, et al, Functioning in the Real World: A PreCalculus Experience, Addison-Wesley, 1997, Reading, MA.

2. Mickins, Ronald, Difference Equations: Theory and Applications, Van Nostrand Reinhold, 1990, New York.

3. Hughes-Hallett, Deborah, Andrew Gleason, et al, Applied Calculus for Business, Social Science, and Life Science Students, John Wiley & Sons, New York,1998.

4. Brown, Douglas and Thomas Timchek, Levodopa to the Brain: A Rate-Change Study (unpublished manuscript).

Acknowledgment

The work described in this article was supported by the Division of Undergraduate Education of the National Science Foundation under grant #DUE-9555401 for the Long Island Consortium for Interconnected Learning. However, the views expressed are not necessarily those of either the Foundation or the project.