Calculus: Project 35

Project 35: Drug Dynamics and Pharmacokinetics

Differential equations are useful in studying the dynamics of a drug in the body. The study of such dynamics is called "pharmacokinetics." Here is a basic example. Suppose a drug is introduced into the blood stream, say by an intravenous injection. The injection rapidly mixes with the whole blood supply and produces a high concentration of the drug everywhere in the blood. Several tissues will readily absorb the drug when its concentration is higher in the blood than in the tissue, so the drug moves into this second tissue "compartment." If this were the whole story, the concentration would eventually balance out so that the concentration in blood and tissue were both equal to the total amount of drug divided by the total volume. However, this usually is not the end of the story.

The kidneys remove the drug from the blood at a rate proportional to the blood concentration. This causes the blood concentration to drop, and eventually it drops below the tissue concentration. At that point, the drug flows from tissue back into the blood and is continually eliminated from the blood by the kidneys. In the long term, the drug concentration tends to zero in both blood and tissue. The speeds with which these various things happen is the subject of pharmacokinetics.

Figure 35.1: two-compartments

Why should we care about such dynamics? Some drugs have undesirable, or even dangerous, side effects if their concentration is too high. At the same time those drugs must be above a certain concentration to be effective for their intended use. As the drug is eliminated from the body, doses need to be given periodically in order to maintain the threshold level for effectiveness, yet doses cannot be too frequent or too large or the concentration will exceed a dangerous level.

This project studies the basic dynamics of the two-compartment drug model. Projects 15 and 16 study some important consequences of these dynamics. The "two-compartment model" corresponds to the previous story. More complicated dynamics might include an intramuscular injection that first diffuses into the blood stream, then into perfuse tissue, and finally is eliminated by the kidneys. Another dynamic might be a drug that is metabolized while it is in the tissue compartment.

The actual "compartments" can vary according to the drug. For example, "blood volume" often includes highly perfused lean tissue such as the heart, lungs, liver, kidney tissue, endocrine glands and occasionally the brain and spinal system (since these often present different barriers to drugs). In these cases, "tissue volume" consists of muscle and skin. Another possible "compartment" is the fat group, adipose and marrow. The anticoagulant warfarin (used as a blood thinner and in rat killer) becomes bound to a protein in the blood. It is not absorbed in tissue, but the bound state forms a second "compartment."

For now we just consider a drug that is introduced into the blood, diffuses into, tissue and is eliminated by the kidneys. Some of the drugs that fit the two-compartment model are aspirin (acetylsalicylic acid); creatinine, a metabolite of creatine produced by muscle contraction or degeneration; aldosterone; griseofulvin (an antifungus drug); and lecithin.

PRIMARY VARIABLES OF THE MODEL

The primary quantity we wish to measure is the drug concentration in the blood. That concentration is a function of time and is quickly and directly affected by the concentration of the drug in the tissue.

% Show that the units of the drug concentrations in mg/l are equal to units of micrograms per milliliter, so that we might measure the number of micrograms in a milliliter of a patient's blood, rather than the number of milligrams in a liter of the blood....

Each drug and each patient have certain important constants associated with them. In the complete story these will have to be measured.

PARAMETERS OF THE MODEL

Concentration is an amount per unit volume. Suppose that a patient has 2.10 liters of blood and that 0.250 grams of a substance is introduced into his blood. When mixed, the concentration becomes .

Secondary Variables
Give a general formula to convert the blood and tissue concentrations into the actual total amounts of the drug in the blood or tissue at a given time,

Suppose a patient has v_B=2.10 and v_T=1.30. If c_B[t] and c_T[t] were known functions, what would the amounts a_B[t] and a_T[t] be?

35.1 Derivation of the Equations of Change

The primary dynamic mechanism for moving the drug into tissue is a "concentration gradient." If the blood concentration is higher than the tissue concentration, there is a flow of the drug from the blood into the tissue. The simplest relationship would be a linear one,

We use the parameter k₁ for the constant of proportionality. Suppose that this patient has no kidneys, so there is no elimination (or metabolism) of the drug. Then we have the rate of decrease of c_B equal to the negative derivative,

% What are the units of ? What does this force the units of k₁ to be?

The drug that flows from blood to tissue is neither created nor destroyed in our model. We need to account for this amount in the equation for the rate of change of the tissue concentration, . Consider the case of our hypothetical patient with a blood volume of 2.10 liters and tissue volume of 1.30 liters. If c_B is reduced by 10.0mg/l in an hour, then the amount of drug that has left in the hour is (Note units: mg/l l = mg). This raises the amount of drug in the tissue by 21mg, but the concentration in the tissue goes up by 21/1.30=16.mg/l. Notice the units in our calculations and use them in the next exercise.

% What is the rate of change of the amount of drug in the blood, , in terms of patient parameters and ?
What is the rate of change of the amount of drug in the tissue, , in terms of patient parameters and ?
Show that

where k₂=k₁v_B/v_T.

Finally, we shall hypothesize a linear elimination law for a patient with kidneys, that is,

so that the two terms reducing blood concentration are

The Dynamics
Begin your project with an explanation of the pharmacological meaning of the equations:

Give the units of your variables and parameters and explain why the parameters v_B, v_T, k₁ and k₃, but not k₂, need to be measured.

Now we want to begin with some numerical experiments to get a feel for the model.

Numerical Experiments with No Kidneys
Modify the AccDEsol program to find some explicit solutions of the concentration equations. Your initial experiment might look as follows: f := -(k1 + k3) cB + k1 cT;
g := k2 cB - k2 cT;
k1 = 0.17;
k2 = k1vB/vT;
k3 = 0.0 (*No kidney function.*);
ti = 0.0 (*Initial time.*);
cBi = 47.7 (*Initial blood concentration.*);
cTi = 0.0 (*Initial tissue concentration.*);
tfinal = 10 ;
dt = .1 ;
AccDEsol[{1,f,g},{t,cB,cT},{ti,cBi,cTi},tfinal,dt]; (*Use the AccDEsol pgm.*)

Once you have your program working properly, compute what happens to a patient with these parameters after 48 hours. Your graphs should look like Figure 35.8.

: Forty-eight Hours without Kidneys

Equilibrium without Kidneys
A patient with no kidneys, v_B=2.10, v_T=1.30, and k₁=0.17 is administered 100mg of a drug that is not metabolized. What will his blood concentration be after 48 hours? (Hints: What is the total amount of drug in his body? What is the total volume? What is the connection between amount and concentration?)
A patient is administered an initial amount a_I of drug. Show that the initial blood concentration is c_I=a_I/v_B. Let be an unknown amount for the equilibrium concentration. Show that the equilibrium amount of drug in the blood is and the equilibrium amount of drug in tissue is , for a total . Set this equal to a_I and solve for ?
The graphs above show that c_B and c_T both tend to a limit as t tends to infinity,

Why do both variables tend to this limit? What is a formula for this limiting concentration, , in terms of patient parameters and initial dose?

Your next task is to view these solutions from the point of view of the flow solutions. This lets us see many initial conditions at once.

The Flow Picture without Kidneys
Modify the Flow2D program to check a number of equilibria, but with many initial doses. Your program might look as follows: f := -(k1 + k3) cB + k1 cT;
g := k2 cB - k2 cT;
k1 = 0.17;
k2 = k1vB/vT;
k3 = 0.0 (*No kidney function.*);
initDoses = Table[{c,0},{c,0,100,5}]
initConcs = initDoses/vB;
tfinal = 48 ;
dt = .1 ;
numgraphs = 10 ;
flow2D[ {f,g}, {cB,cT}, initConcs, tfinal, dt, numgraphs ];

Figure 35.2: No Kidney Equilibria for Various Doses

Conjecture
Describe the flow animation of the patient with no kidneys. What are the apparent paths of the solutions? What is the locus of points where the solutions end up? Is there more than one equilibrium point?

The solutions of the flow animation are surprisingly simple and yet not "typical" of the main examples studied in Chapters 23 and 24 of the core text. However, what you see is indeed correct. Things do flow along parallel lines toward another line of equilibria in the phase plane.

Symbolic Proofs for a Degenerate System

1. Show that the constant solutions , for any constant , are equilibrium points of the dynamical system corresponding to no kidneys:

2. Show that the following is an analytical solution to our model:

where c₁ and are constant and h=k₁+k₂. In other words, plug the solutions and into the equations and verify that they do solve them. (Hint: c_B-c_T=(k₁+k₂)c₁e^{-(k₁+k₂)t}.)
3. Show that these analytical solutions satisfy .
4. Review the derivation of the equation for a vector parametric line X=P+tD (Chapter 17 of the core text) and use that idea to show that the solution

gives the motion along the parallel lines to the other line that you conjectured in the previous exercise. Notice that the scalar coefficient e^-ht starts at 1 and decreases to 0 as .

5. If the initial dose is the amount a_I, so c_B[0]=a_I/v_B and c_T[0]=0, show that and the constant c₁ from part (2) must be .

In summary, show that the unique solution to the drug equations without kidneys with initial condition (c_I,0) may be written in terms of the initial concentration, c_I, and parameters as

where

At this point we have a complete numerical and analytical solution to the dynamics of a patient without kidneys. This is a mathematical start but not the main idea we want to explore. However, it is a useful start, because elimination is typically a slower process than tissue absorption. You will see some similarity in the solutions with a small k₃ compared with k₃=0.

Numerical Experiments with Kidneys
Modify the AccDEsol and Flow2D programs to solve the drug equations when . For example, f := -(k1 + k3) cB + k1 cT;
g := k2 cB - k2 cT;
k1 = 0.17;
k2 = k1vB/vT;
k3 = 0.03;
ti = 0.0 (*Initial time.*);
cBi = 45.0 (*Initial blood concentration.*);
cTi = 0.0 (*Initial tissue concentration.*);
tfinal = 72 ;
dt = 0.5 ;
AccDEsol[{1,f,g},{t,cB,cT},{ti,cBi,cTi},tfinal,dt]; (*Using the AccDEsol pgm.*)
initConcs = Table[{c,0},{c,0,100,5}];
numgraphs = 10 ;

flow2D[ {f,g}, {cB,cT}, initConcs, tfinal, dt, numgraphs ];

: Seventy-two Hours with Kidneys

Equilibrium with Kidneys
If k₁>0 and k₃>0, prove that (0,0) is the only equilibrium for the drug equations.

Explain intuitively why this equilibrium must be an attractor.

Notice that our concentration equations have the form of a linear system

where a_x=-(k₁+k₃), a_y=k₁, b_x=k₂, and b_y=-k₂. This implies that the concentration solutions may be written in terms of linear combinations of the exponential functions e^-h₁t and e^-h₂t where -h₁ and -h₂ satisfy the characteristic equation:

The Analytical Form of the Solutions

1. Use Theorem 24.1 from Chapter 24 of the core text to show that the solutions c_B[t] and c_T[t] may be written as a linear combination of e^-h₁t and e^-h₂t,

where

2. Use the Local Stability Theorem from Chapter 24 to prove that (0,0) is an attractor when k_j>0, for j=1,2,3. (Typically, k₃<<k₁, although this is not needed for your proof.)
3. Show that

4. What happens to h₁ and h₂ in the case of nonfunctioning kidneys (k₃=0)?
5. Let k₁=0.17, v_B=2.10, v_T=1.30. Compute h₁ and h₂ in two cases: k₃=0.03 and k₃=0.0.
What we know from the theory is that

for some constants b₁ and b₂, where h₁ and h₂ are computed from the parameters by the quadratic formula of the previous exercise. The exact values of b₁ and b₂ are given in Hint 35.1.14. The details are messy but amount to some formulas that the computer can easily compute for us.
Reasonable approximations for b₁ and b₂ are easy to guess in the case h₁>>h₂. If the fast absorption into the tissue is so fast that the kidneys do nothing during this period, we know that the solution looks like where

% Justify this last claim, , for the constants above. What are the necessary assumptions? (Hint: See a previous exercise.)
Show that k₁>>k₃ implies that h₁>>h₂ for some reasonable meaning of >>. (If you want to be formal, take k₁ (and k₂) positive but not infinitesimal and . Show that , while . Connect your formalization with the real problem.)
In the case where k₃=0, the full solution could be written

The actual exact symbolic solution when can be written in the vector form

% Show that if , the two solutions immediately before this exercise are close to one another for finite time.
If we let

then we get the initial conditions

Show that this choice of constants c₁ and c₂ immediately before this exercise may be written in the form

where

These formulas give the exact symbolic unique solution of the full drug dynamics equation with initial condition (c_I,0)

% Compare the approximate solution function c_A[t]=u₁e^-h₁t+u₂e^-h₂t to the AccDEsol solution of the drug equations with initial concentration c_I=1 and to the exact symbolic solution c_B[t]=w₁e^-h₁t+w₂e^-h₂t. (Use the values of the parameters k₁, k₃, etc., from the exercises above.) Calculate the associated constants u₁ versus w₁ and u₂ versus w₂, using the formulas from Exercise 35.1.8 and the previous exercise.
What happens when k₃=0? What happens when k₃ is nearly as large as k₁?
k1 = 0.17;
vB = 2.1;
vT = 1.3;
k3 = 0.03;
k2 = k1 vB/vT;
h1 = ((k1 + k2 + k3) + Sqrt[(k1 + k2 + k3)2 - 4 k2 k3])/2;
h2 = ((k1 + k2 + k3) - Sqrt[(k1 + k2 + k3)2 - 4 k2 k3])/2;
u1 = vT/(vT + vB)
w1 = (h1 - k2)/(h1 - h2)
u2 = vB/(vT + vB)
w2 = (k2 - h2)/(h1 - h2)
cB[t_] := cI (w1 Exp[ -h1 t] + w2 Exp[ -h2 t]);
cB[t]

: Computer Solution and Approximate Coefficients

35.2 Where Do We Go from Here?
The exact symbolic solution to a single dose [or initial conditions (c_I,0)] given in the previous exercise is the starting point for the Project 15. That project explores things like the peak tissue concentration, the time interval during which the concentration remains effective, and the total integral of effective concentration. You could go to that project now to round out your thinking on pharmacokinetics.
Project 16 uses this single-dose symbolic solution and some basic work on logs from Project 15 to show how drug concentration data can be used to measure the rate parameters of this model. You might go to that project now instead of working Project 15.
Other alternatives for completing this project are to study more complicated dosing regimens. Usually a single dose of a drug is not given, but rather doses are either given periodically or a constant flow of drug is maintained in an intravenous fluid. The analysis of these two regimens follows as your third and fourth alternatives. Finally, a fifth alternative completion of the project is to study dosing by intramuscular injection. In this case, the drug must diffuse from the muscle into the blood before the two-compartment dynamics takes effect. This amounts to a three-compartment model.
35.3 Periodic Intravenous Injections (Optional Project Conclusion)
Suppose a patient is given a first dose of a drug, a₁ (mg) at time t=0. Typically, this dose is larger than subsequent doses. After a time interval , say 6 hours, a second dose a₂ is administered. We know that the exact solution from time t=0 to is given by the exponential sum,

where c_I=a₁/v_B and

or we could simply use AccDEsol to find the solution.
At time , the blood concentration is increased by a₂/v_B, but the old drug is not absent from the blood (and we often want to be sure it does not drop below a threshold of effectiveness). Say the concentration just before the second dose is (using the formula above for c_B[t].) Our new initial concentration becomes .
Also, the concentration in the tissue is not zero at time . We can use the formula above to compute . As a result, we want to solve the new initial value problem

As an example, let us take a₁=200 (mg) and a₂=100 (mg) using the parameters k₁=0.17, v_B=2.10, v_T=1.30, k₃=0.03, as before. The values of concentration before the second dose are c_B[6]=52.6 and c_T[6]=49.7.

Figure 35.3: Concentrations up to Time 6
This makes the new initial value problem beginning with the second dose

Notice that the initial value of c_T is no longer zero.
We can put this into AccDEsol as a new initial value problem and solve for six more time units, then repeat the process. Or we can find the exact symbolic solution to this initial value problem and compute the concentrations six units later, add the third dose and continue with that procedure. In either case we get the long-term drug dosages and the question is: Will the concentrations build up? What is the long-term max? The min?
That's your project.
Here's some help if you want to use exact solutions which aren't really any better than numerical ones from AccDEsol. The general solution to the differential equations may be written in the vector form

If we reset our clock at the second dose, we want to have

Notice that substitution of t=0 into the general vector form yields

or in matrix form

For example, these equations can be solved with Mathematicaby typing e1 = 100.;
e2 = 49.7; {c1,c2} = LinearSolve[{{k1,k1},{h2 - k2 , h1 - k2}},{e1,e2}]
cB[t_] := c1 k1 Exp[ -h1 t] + c2 k1 Exp[-h2 t];
cT[t_] := c1 (h2 - k2) Exp[ -h1 t] + c2 (h1 - k2) Exp[-h2 t];
cB[6]
cT[6]
Our calculations with the parameters and doses above gave the following:

: Concentrations between Doses 2 -- 3, 3 -- 4, 4 -- 5, 5 -- 6, 6 -- 7, 7 -- 8

What are the concentrations at the end of 48 hours, just before the ninth dose? Do the concentrations continue to build up over time, or do they reach some maximum level and stop growing?
35.4 Steady Intravenous Flow
The previous problem started with an initial dose of 200 mg and then followed that with doses of 100 mg every 6 hours. What would happen to the patient's drug concentrations if we constantly fed the drug into the blood stream at the rate of 100/6 = 16.7 mg/hr? This is much easier for us to analyze , because the "forcing term" or dosage is continuous.
First, we need to modify the basic equations,

Now we have an additional term of a constant growth r added to the blood concentration. In the case described in the previous paragraph, r = 16.7/vB = 7.94 mg/l/hr. In general, this constant dose rate makes our initial value problem:

% Solve this initial value problem (using the computer with AccDEsol or Flow2D or an exact method by changing variables) and show that its equilibrium level is

: A Year of 16.7 mg/l/hr Dosing

It seems reasonable to guess that the average concentration of the drug levels in the patient dosed every 6 hours (in the previous example) would build up to the equilibrium level of the continuously dosed patient. How high is this level?
35.5 Intramuscular Injection - A Third Compartment
Our final example of dosing is an intramuscular injection. In this case the drug is injected in high concentration in a muscle and then the drug diffuses into the blood. Once in the blood, it diffuses into the tissue compartment and is eliminated by the kidneys.
PARAMETERS FOR IM INJECTION v_M - the volume of the injection site k₄ - the rate parameter for diffusion into the blood k₅=k₄v_M/v_B

% Explain the meaning of the system of differential equations:

in terms of the intramuscular injection. What is the amount of drug in the injection?
Use the computer to solve this system with v_M=0.250 and k₄=1.5 and the parameters from the previous examples.

Previous project Next project Close this window