Exploring Physiologically Based Pharmacokinetic Models

The author presents a Mathematica-based model for assessing the risks of styrene inhalation.

by Qi Zheng

Physiologically based pharmacokinetic (PBPK) models are rapidly winning   acceptance  among   risk analysts  and  regulatory agencies (see, e.g., [Travis 1989] and [Travis  et al. 1989]). PBPK models are in essence deterministic compartmental models employing  ordinary differential equations (ODE's) to explicitly describe the flow, disposition, absorption and metabolism within an animal or human body of the  chemical under consideration. A typical PBPK model consists of from a few to several dozens of  physiologically meaningful compartments representing the lung, fatty tissue, the liver,  etc.. Like any mathematical models mimicking complex biological systems, PBPK models require that many important physiological and biochemical parameters be known a priori or be directly measured in vivo. (Statistical estimation of some biochemical parameters are also commonly practiced, see, e.g., [Smith and Evans 1995]). Although prototypes of PBPK models can be traced many decades back, it was not until the 1980's that PBPK models became an important and practical component of  risk assessment when accurate measurement of the crucially important physiological and biochemical parameters ( e.g. alveolar ventilation rates, equilibrium partition coefficients and metabolic rates) became feasible. For more information about experimental measurement of such parameters, see [Gargas et al. 1986] and [Gargas et al. 1989]. 

Recent years have seen an explosive increase in theoretical research and applications of PBPK models (e.g., [Reitz et al. 1988],  [Corley et al. 1990], [Andersen et al. 1991], [Allen and Fisher 1993] and [Buckley 1995]), and consequently an array of slightly different versions of PBPK models burst into bloom. The earliest and simplest model that fits the currently accepted ``definition" of a typical PBPK model is perhaps that proposed by Ramsey and
Andersen [Ramsey andAndersen 1984] in their pioneering work on styrene inhalation pharmacokinetics.Their model (Fig. 1) consists of only 5 compartments and yet possesses all the salient features of today's PBPK models, and is therefore the focus of the present notebook for the sake of simplicity. For more details concerning the interpretation of the model depicted in Figure 1, the reader should consult the original paper by Ramsey and Andersen.

Figure 1:  THe PBPK model of Ramsey and Andersen.

Note that  the word "chemical" is employed throughout this notebook as a generic term to refer to any specific substance of interest (e.g., drugs, carcinogens, etc.).

The Mass Balance Equations

The starting point of PBPK modeling is the setting up of mass balance  equations.  For each organ or tissue represented by a  compartment, the  rate at which  the total amount of  chemical in the particular compartment changes is the difference between the rate of in-flow (through arterial blood) and the rate of out-flow (through venous blood). Take the fatty tissue compartment for example. The rate of in-flow of chemical is the velocity of blood running through this compartment  (QFat) times the  chemical concentration in the arterial blood (ConcArterical[t]).  The rate of out-flow  is calculated by analogy. Therefore, the mass balance equation for the fat compartment is of the form

AmtFat'[t]==QFat (ConcArterial[t]-ConcVenousFat[t]).

Since the lung compartment serves only as an interface to the outside exposure (styrene inhalation in this example), there are only four mass balance equations for the four compartments: the fatty tissue, the muscle, the richly perfused tissue and the liver.

MassBalanceEq={
 AmtFat'[t]==QFat (ConcArterial[t]-ConcVenousFat[t]),

AmtMuscle'[t]==QMuscle (ConcArterial[t]-
 ConcVenousMuscle[t]),

AmtRichly'[t]==QRichly (ConcArterial[t]-
 ConcVenousRichly[t]),

AmtLiver'[t]==QLiver (ConcArterial[t]-ConcVenousLiver[t])-(Vmax ConcVenousLiver[t])/(Km+ConcVenousLiver[t])}

Note that the last equation (for the liver) contains an extra term describing metabolism of styrene taking place in the liver (to be more precise, the word "liver" in this context includes the liver per se and other tissues with similar metabolic capacity). The metabolism is assumed to be saturable, satisfying the well-known Michaelis-Menten equation (Vmax is the maximunm reaction rate and Km is the Michaelis constant).
At first sight, the above four mass balance equations involve  nine unknowns. However, concentration of the chemical in a particular compartment  is related to the total amount of the chemical residing in that compartment (that is, Amt/Vol=concentration), and concentration in the venous blood running out of a particular compartment is linked  to the concentration in that particular compartment by the tissue:blood  partition coefficient.  

Of particular note in this regard is the assumption that the chemical entering a compartment via arterial blood distributes throughout the compartment  very rapidly, and hence an equilibrium ratio (called tissue:blood partition coefficient) of the concentration within the compartment  to  the concentration of the venous blood running out of this same compartment  is reached almost instantaneously. The following rules reflect this important assumption.  Application of these rules will accomplish the task of expressing the four venous blood concentration variables in terms of their corresponding total amount variables.

ConcToAmtRules={
ConcVenousFat[t]->AmtFat[t]/VolFat/PartCoefFat,
ConcVenousMuscle[t]->AmtMuscle[t]/VolMuscle/PartCoefMuscle,
ConcVenousRichly[t]->AmtRichly[t]/VolRichly/PartCoefRichly,
ConcVenousLiver[t]->AmtLiver[t]/VolLiver/PartCoefLiver};

Chemical enters the circulation system by  the arterial blood's absorbing the chemical from alveolar air in the lung. The air:blood partition coefficient bridges concentration in arterial blood with that in alveolar air, which in turn is directly linked with inhalation exposure. First, the simple relation between concentration in alveolar air and that in arterial blood.

AlveolarRule=ConcAlveolar[t]->ConcArterial[t]/PartCoefAirBlood;

The relation between inhalation exposure level and the total amount of chemical absorbed by the body is summarized in the the following VenousConcEqn equation.

VenousConcEqn=
QAlveolar (ConcInhale-ConcAlveolar[t]) == 
  Qt (ConcArterial[t]-ConcVenous[t])/.AlveolarRule;

This equation simply says that the total amount of chemical lost from the alveolar air at any moment is the same as the total amount of chemical absorbed by all tissues from  the  blood (measured as the difference between arterial blood and venous blood). Note that QAlveolar and Qt are the flow speeds of alveolar air and blood circulation through the lung.

SolnToVenousConcEqn=Solve[{VenousConcEqn},{ConcArterial[t]}][[1]];

In the above expression the variable ConcVenous[t] is the concentration in venous blood coming out of the lung, it is related to the concentrations in various compartments  in the following manner.

ConcVenousRule=ConcVenous[t]->(QFat ConcVenousFat[t]+
 QMuscle ConcVenousMuscle[t]
 +QRichly ConcVenousRichly[t]+QLiver ConcVenousLiver[t])/Qt;

Now one can express the ConcArterial variable in the four basic  variables using the following rule.

ArterialRule=SolnToVenousConcEqn/.ConcVenousRule;

Physiological parameter values must be incorporated into  the mass balance equations. Table 2 of [Ramsey and Andersen 1984] provides all these parameter values for mice, rats and humans. For illustration purposes, this section  deals with  rat models only. The interested reader can easily replace rat parameters with  mouse parameters. Human models will be dealt with in the subsequent section.

RatParameters={
QAlveolar->4.50,Qt->5.64,QLiver->2.11,QFat->0.5,
QMuscle->0.68,QRichly->2.35,VolLiver->0.012,VolFat->0.027,
VolMuscle->0.219,VolRichly->0.015,
PartCoefLiver->2.7,PartCoefFat->50,PartCoefMuscle->1.0,
PartCoefRichly->5.7,
Vmax->3.6,Km->0.36,
PartCoefAirBlood->40.0};

All the foregoing work now results in an ODE system of mass balance equations describing the pharmacokinetics of inhaled styrene within a laboratory rat.

RatODE=MassBalanceEq/.ArterialRule/.ConcToAmtRules/.RatParameters

Simulation Using  Rat Models

Now the four mass balance equations contain exactly four unknowns except for an indeterminate function ConcInhale, which is in fact the force function of the whole system to be determined by the inhalation exposure pattern, viz., how the external exposure changes with time. Once ConcInhale is specified, the ODE system of mass balance equations can be numerically solved, and various quantities of interest can be readily obtained. 

The objective of this section is to simulate the time course of styrene concentrations in the four tissue groups and the blood for some prescribed external exposure patterns. Following [Ramsey and Andersen 1984], interest is focused on eight-hour constant exposures with dose levels  ranging from 80ppm to 1200ppm. Since the concentration within  a compartment is expressed in mg/l rather than ppm as for inhalation concentration in the air, a function converting ppm to mg/l is handy.

ppmToMgl[ppm_]:= ppm 104.151/10^3/0.0821/298

In the above ppmToMgl function, room  temperature (25C=298K) is assumed. The molacular weight of styrene is taken to be 104.151, and the gas constant is taken to be 0.0821 liter*atm/mole/K (c.f. [Masterton and Slowinski 1978, p.99]). For example, the following translates a concentration of 1200ppm of styrene into the corresponding concentration  expressed in mg/l.

ppmToMgl[1200]

Assuming that the rat has no previous exposure to styrene, that is, there is no styrene residue in the body of the rat, then  the following initial condition is self-evident.

initCond={AmtFat[0]==0,AmtMuscle[0]==0,
AmtRichly[0]==0,AmtLiver[0]==0};

Now a variety of exposure patterns can be investigated by simply specifying appropriate inhalation dose patterns. The dose patterns investigated by Ramsey and Andersen (1984) are two-piece step functions, that is, the rat is exposed to a constant styrene concentration (say, 1200ppm) in the time interval [0,6), and exposure is discontinued after 6 hours of exposure. This exposure pattern is now incorporated into the mass balance equations.

ExposurePattern={ConcInhale->If[t<6.0,ppmToMgl[1200.],0.]};
FinalRatODE=Join[RatODE,initCond]/.ExposurePattern//Simplify

Numerical solutions to the above ODE system can now be found to yield simulated styrene concentrations in all four tissue groups described by the mass balance equations. A 24-hour period is chosen for the simulation study, as was done by [Ramsey and Andersen 1984].

ODENames={AmtFat,AmtMuscle,AmtLiver,AmtRichly};
RatSoln=NDSolve[FinalRatODE,ODENames,{t,0,24.},MaxSteps->5000]

The simulation results contained in the above ODE solutions are best presented by graphics.

Needs["Graphics`Graphics`"];

The time course of styrene concentrations in the four tissue groups are now presented graphically.

pFat=Plot[Evaluate[AmtFat[t]/VolFat/.RatParameters/.RatSoln],{t,0,24},
FrameLabel->{"Time (Hours)","Concentration (mg/l)"},  
Frame->True,DisplayFunction->Identity,
PlotLabel->FontForm["Fatty Tissue",{"Helvetica",12}],
FrameTicks->{Range[0,24,4],Automatic}]; 

pLiver=Plot[Evaluate[AmtLiver[t]/VolLiver/.RatParameters/.RatSoln],{t,0,24},
FrameLabel->{"Time (Hours)","Concentration (mg/l)"},  
Frame->True,DisplayFunction->Identity,
PlotLabel->FontForm["Liver",{"Helvetica",12}],
FrameTicks->{Range[0,24,4],Automatic}];

pRichly=Plot[Evaluate[AmtRichly[t]/VolRichly/.RatParameters/.RatSoln],{t,0,24},
FrameLabel->{"Time (Hours)","Concentration (mg/l)"},  
Frame->True,DisplayFunction->Identity,
PlotLabel->FontForm["Richly Perfused Tissue",{"Helvetica",12}],
FrameTicks->{Range[0,24,4],Automatic}];

pMuscle=Plot[Evaluate[AmtMuscle[t]/VolMuscle/.RatParameters/.RatSoln],{t,0,24},
FrameLabel->{"Time (Hours)","Concentration (mg/l)"},  
Frame->True,DisplayFunction->Identity,
PlotLabel->FontForm["Muscle",{"Helvetica",12}],
FrameTicks->{Range[0,24,4],Automatic}];

Show[GraphicsArray[{{pFat,pLiver},{pRichly,pMuscle}}]];

The above graph offers a "faithful" presentation of the 24-hour period of styrene concentration in the four tissue groups. However, in pharmacokinetics literature,  graphics using semilog coordinates are  more commonly adopted. In the following, the built-in function LogPlot will be employed to comply with this practice.

pLogFat=LogPlot[Evaluate[AmtFat[t]/VolFat/.RatParameters/.RatSoln],{t,0,24},
Frame->True,PlotRange->{.1,3000},PlotStyle->Dashing[{0.01}],
FrameTicks->{Range[0,24,4],Automatic},
DisplayFunction->Identity]

Note that, since the concentration is zero at time t = 0, the error message simply informs that the logarithm of zero is a complex number, and hence the message can be safely ignored. Another variable of great interest is styrene concentration in (arterial) blood. Although this quantity is not directly involved in the mass balance equations, it is already defined in the previous section.

ConcArterialBlood=ConcArterial[t]/.ArterialRule/.
  ConcToAmtRules/.ExposurePattern/.RatParameters

pLogBlood=LogPlot[Evaluate[ConcArterialBlood/.RatSoln],{t,0.001,24},
Frame->True,PlotRange->{0.1,100},
FrameTicks->{Range[0,24,4],Automatic},
DisplayFunction->Identity]

Combining  the two concentration curves reproduces one of the graphics in Figure 2 of Ramsey and Andersen (1984).

Show[pLogFat,pLogBlood,Graphics[Text["Blood",{15,1.3}]],
Graphics[Text["Fat",{19,2.9}]],
PlotLabel->FontForm["Styrene Concentration in Rats",
   {"Palatino-Bold",18}],
FrameLabel->{" Time (hours)","Log Concentration (mg/l)"},
DisplayFunction->$DisplayFunction]

Simulation Using Human Models

The most encouraging aspect of PBPK models is that they allow extrapolation of pharmacokinetic behavior mong different species by appropriately scaling the physiological and biochemical parameters. Ramsey and Andersen provide parameters for humans solely by rescaling the parameters for rats. For simplicity, the parameters for humans in Table 2 of Ramsey and Andersen are taken directly without repeating the allometric scaling process.

HumanParameters={
QAlveolar->230,Qt->290,QLiver->108,QFat->26,
QMuscle->34,QRichly->121,VolLiver->3.32,VolFat->7.47,
VolMuscle->60.6,VolRichly->4.15,
PartCoefLiver->2.7,PartCoefFat->50,PartCoefMuscle->1.0,
PartCoefRichly->5.7,
Vmax->184,Km->0.36,
PartCoefAirBlood->52.0};

The mass balance equations for humans are obtained by simply replacing the rat parameters with the human parameters in the  mass balance equations for rats obtained in the second section.

HumanODE=MassBalanceEq/.ArterialRule/.ConcToAmtRules/.HumanParameters;

To keep the simulation simple, all conditions are kept  the same except that the external dose is now 80ppm and the simulation is now performed on a 48-hour period.

HumanExposurePattern={ConcInhale->If[t<6.0,ppmToMgl[80.],0.]};
HumanFinalODE=Join[HumanODE,initCond]/.HumanExposurePattern//Simplify;
HumanSoln=NDSolve[HumanFinalODE,ODENames,{t,0,48.},MaxSteps->5000]

Styrene concentration in the venous blood is again of great interest. The expression remains the same as that for rats, only some physiological and biochemical parameters need changing.

HumanConcArterialBlood=ConcArterial[t]/.ArterialRule/.
  ConcToAmtRules/.HumanExposurePattern/.HumanParameters

pHuman1=LogPlot[Evaluate[HumanConcArterialBlood/.HumanSoln],{t,0,48.},
PlotRange->{0.00001,10},Frame->True,
FrameTicks->{Range[0,48,8],Automatic},
DisplayFunction->Identity];

Another variable of interest in studying the human models is styrene concentration of exhaled air. According to the assumption of Ramsey and Andersen, alveolar respiration accounts for 70% of total respiration, and hence the following definition.

ConcExhaledAir=(0.7 ConcAlveolar[t]+0.3 ConcInhale)/.AlveolarRule/.ArterialRule/.
ConcToAmtRules/.HumanParameters/.HumanExposurePattern

pHuman2=LogPlot[Evaluate[ConcExhaledAir/.HumanSoln,{t,0,48}],
PlotRange->{0.00001,1},
Ticks->{Range[0,48,8],Automatic},
DisplayFunction->Identity];

Combining the above two curves gives a figure similar to Figure 5 of Ramsey and Andersen (1984).

Show[pHuman1,pHuman2,Graphics[Text["Exhaled Air",Scaled[{0.5,0.3}] ]],
  Graphics[Text["Blood",Scaled[{0.7,0.55}]  ]],
PlotLabel->FontForm["Styrene Concentration in Humans",
           {"Palatino-Bold",18}],
FrameLabel->{" Time (hours)","Log Concentration (mg/l)"},
DisplayFunction->$DisplayFunction];

Finally, consider the simulation of an 8-hour workday styrene exposure at 80 ppm dose level for a week of five workdays. This simulation is of some practical interest because human exposure to styrene occurs mainly in the industrial workplace.

WorkdayExposure=ConcInhale->If[t<8||
(t>24)&&(t<32)||(t>48)&&(t<56)||
 (t>72)&&(t<80)||(t>96)&&(t<104),ppmToMgl[80],0];
 
 WorkdayODE=
 Join[HumanODE,initCond]/.WorkdayExposure//Simplify;
 
 WorkdaySoln=NDSolve[WorkdayODE,ODENames,{t,0,168.},
  MaxSteps->5000]

ConcWorkdayBlood=ConcArterial[t]/.ArterialRule/.
  ConcToAmtRules/.WorkdayExposure/.HumanParameters

For comparison purposes, the semilog scale in not adopted here. The following graph is similar to Figure 4 of Ramsey et al. (1980), although the two were drawn using slightly different models.

Plot[Evaluate[ConcWorkdayBlood/.WorkdaySoln],
 {t,0,168.},Frame->True,
 PlotLabel->FontForm["Human Exposure in Workplace",
   {"Palatino-Bold",18}],
FrameLabel->{"Time (hours)","Concentration (mg/l)"},
FrameTicks->{Range[0,168,24],Automatic}  ];

Closing Remarks

The Ramsey-Anderson  model presented here is only one of the simple yet mature PBPK models. In principle, PBPK models are chemical dependent in that both the structure of the model and the biochemical parameters are determined by the chemical under investigation. Following this notebook, the reader can try many more sophisticated PBPK models alone the same line. The unique value of this Mathematica notebook lies in that it makes all the conceptually important details transparent, yet the reader need not worry too much about the purely mathematical technicalities.

Bibliography

Andersen, M.E., Clewell, H.J. III, Gargas, M.L., MacNaughton, M.G., Reitz, R.H., Nolan, R.J., and McKenna, M.J.  Physiologically based pharmacokinetic modeling with dichloromethane, its metabolite, carbon monoxide and blood carboxyhemoglobin in rats and humans, Toxicology and Applied Pharmacology 108:14-27 (1991).

Allen, B.C.  and Fisher, J.W.  Pharmacokinetic modeling of
trichloroethylene and trichloroatcetc acid in humans,  Risk Analysis 13:71-86 (1993).

 Buckley, L.A., Biologically-based models of dioxin pharmacokinetics, Toxicology 102:125-131 (1995).

Gargas, M.L., Andersen, M.E., and Clewell, H.J. III, A physiologically based simulation approach for determining metabolic constants from gas uptake data. Toxicology and Applied Pharmacology 86:341-352 (1986).

Gargas, M.L., Burgess, R.J., Voisard, D.E., Cason, G.H., and Andersen, M.E., Partition coefficients of low-molecular-weight volatile chemicals in various liquids and tissues. Toxicology and Applied Pharmacology 98:87-99 (1989).

Corley, R.A., Mendrala, A.L., Smith, F.A., Staats, D.A., Gargas, M.L., Conolly, R.B., Andersen, M.E.  and Reitz, R.H., Development of a physiologically based pharmacokinetic model for chloroform, Toxicology and Applied Pharmacology 103:512-527 (1990).

Masterton, W.L.  and Slowinski, E.J., Chemical Principles with Qualitative Analysis, W.B. Saunders Company, Philadelphia, 1978.

Ramsey,  J.C.  and Andersen, M.E., A physiologically based
description of the inhalation pharmacokinetics of styrene in rats and humans. Toxicology and Applied Pharmacology 73:159-175 (1984).

Ramsey, J.C., Young, J.D., Karbowski, R.J., Chenoweth, M.B., McCarty, L.P.  and Braun, W.H., Pharmacokinetics of inhaled styrene in human volunteers. Toxicology and Applied Pharmacology 53:54-63(1980).

Reitz, R.H.,  McDougal, J.N., Himmelstein, M.W.,  Nolan, R.J.  and Schumann, A.M., Physiologically based pharmacokinetic modeling with methylchloroform: Implications for interspecies, high dose/low dose, and dose route extrapolations. Toxicology and Applied Pharmacology 95:185-199 (1988).

Smith, A.E.  and Evans, J.S. , Uncertainty in fitted estimates of apparent in vivo metabolic constants for chloroform, Fundamental and Applied Toxicology 25:29-44 (1995).

 Travis, C. C. (ed), Biologically Based Methods for
Cancer Risk Assessment, Plenum Press, 1989.
 
 Travis, C.C., White, R.K., and Arms, A.D., A physiologically
based pharmacokinetic approach for assessing the cancer risk of 
tetrachloroethylene, In Paustenback, D. J. (ed),The Risk Assessment of Environmental and Human Health Hazards,  Wiley, 1989.

ABOUT THE AUTHOR

Qi Zheng is a researcher at the National Center for Toxicological Research.

(*address to come*)

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file Styrene.nb containing Mathematica code for the material described in this article.