Compartmental Modeling in The Analysis of Biological Systems

See discussions, stats, and author profiles for this publication at: https://www.researchgate.
net/publication/286537976
Compartmental modeling in the analysis of biological systems
Article in Methods in molecular biology (Clifton, N.J.) · January 2012

DOI: 10.1007/978-1-62703-50-2_17
CITATIONS READS
6 1,394
4 authors, including:
James Bassingthwaighte Bartholomew E Jardine

University of Washington Seattle University of Washington Seattle
213 PUBLICATIONS 5,912 CITATIONS 9 PUBLICATIONS 58 CITATIONS
SEE PROFILE SEE PROFILE
Gary Raymond
University of Washington Seattle
15 PUBLICATIONS 845 CITATIONS
SEE PROFILE
Some of the authors of this publication are also working on these related projects:
Hepatic Clearance of Indocyanine green View project
All content following this page was uploaded by James Bassingthwaighte on 03 May 2016.
The user has requested enhancement of the downloaded file.

Chapter 17
Compartmental Modeling in the Analysis

of Biological Systems
James B. Bassingthwaighte, Erik Butterworth, Bartholomew Jardine,
and Gary M. Raymond
Abstract
Compartmental models are composed of sets of interconnected mixing chambers or stirred tanks. Each
component of the system is considered to be homogeneous, instantly mixed, with uniform concentration.
The state variables are concentrations or molar amounts of chemical species. Chemical reactions, transmem-
brane transport, and binding processes, determined in reality by electrochemical driving forces and con-
strained by thermodynamic laws, are generally treated using first-order rate equations. This fundamental
simplicity makes them easy to compute since ordinary differential equations (ODEs) are readily solved
numerically and often analytically. While compartmental systems have a reputation for being merely descrip-
tive they can be developed to levels providing realistic mechanistic features through refining the kinetics.
Generally, one is considering multi-compartmental systems for realistic modeling. Compartments can be
used as “black” box operators without explicit internal structure, but in pharmacokinetics compartments are
considered as homogeneous pools of particular solutes, with inputs and outputs defined as flows or solute
fluxes, and transformations expressed as rate equations.
Descriptive models providing no explanation of mechanism are nevertheless useful in modeling of many
systems. In pharmacokinetics (PK), compartmental models are in widespread use for describing the
concentration–time curves of a drug concentration following administration. This gives a description of
how long it remains available in the body, and is a guide to defining dosage regimens, method of delivery,
and expectations for its effects. Pharmacodynamics (PD) requires more depth since it focuses on the
physiological response to the drug or toxin, and therefore stimulates a demand to understand how the
drug works on the biological system; having to understand drug response mechanisms then folds back on
the delivery mechanism (the PK part) since PK and PD are going on simultaneously (PKPD).
Many systems have been developed over the years to aid in modeling PKPD systems. Almost all have
solved only ODEs, while allowing considerable conceptual complexity in the descriptions of chemical
transformations, methods of solving the equations, displaying results, and analyzing systems behavior.
Systems for compartmental analysis include Simulation and Applied Mathematics, CoPasi (enzymatic
reactions), Berkeley Madonna (physiological systems), XPPaut (dynamical system behavioral analysis),
and a good many others. JSim, a system allowing the use of both ODEs and partial differential equations
(that describe spatial distributions), is used here. It is an open source system, meaning that it is available for
free and can be modified by users. It offers a set of features unique in breadth of capability that make model
verification surer and easier, and produces models that can be shared on all standard computer platforms.
Brad Reisfeld and Arthur N. Mayeno (eds.), Computational Toxicology: Volume I, Methods in Molecular Biology, vol. 929,
DOI 10.1007/978-1-62703-050-2_17, # Springer Science+Business Media, LLC 2012
391
392 J.B. Bassingthwaighte et al.
Key words: Physiological and pharmacologic modeling, PKPD, Pharmacokinetics–pharmacody-

namics, Compartmental systems, Systems biology, Physiome, JSim, CellML, SBML, Reproducible
research, Unit checking, Verification, Validation, Ordinary and partial differential equations,
Optimization, Confidence limits
1. Introduction
1.1. Overview Compartmental analysis implies the use of linear first-order differential
of the Topic operators as analogs for describing the kinetics of drug distribution
and elimination from the body. Concentrations are measured in
accessible fluids, usually the plasma, and the concentration–time
curve is used to provide a measure of how long the drug concentra-
tion remains at a therapeutic level. This is the basis of pharmacoki-
netics (PK). The influences on efficacy and utility are considered by
the term ADME, administration, distribution, metabolism, and
elimination. Drugs are given in chemically significant amounts;
they bind to enzymes, channels, receptors or transporters, changing
reaction rates and fluxes in concentration-dependent fashion. The
effects on the biological system is termed pharmacodynamics (PD).
Precise mathematical statements about the kinetics and the body’s
responses comprise the combination pharmacokinetics–pharmaco-
dynamics (PKPD).
Compartmental analysis had its historical start with the use of
tracers. Tracer-labeled compounds were used in order to determine
kinetics when the drug concentrations were too low to be measured
chemically. Radioactive tracers were given in such low concentra-
tions relative to those of native non-tracer mother substances that
the kinetics were in fact linear. Consider a reaction rate, k(C), that is
dependent upon the concentration of the mother substance of
concentration C(t):
Flux of mother substance ¼ kðCÞ CðtÞ: (1)
When tracer of concentration C 0 is added to the system, then
Total flux; mother, and tracer ¼ kðC þ C 0 Þ ½CðtÞ þ C 0 ðtÞ: (2)
When C 0 << C, then the rate constant is determined solely by
C, as k(C þ C 0 ) k(C), and the rate constant is independent of
the tracer concentration:
Tracer flux of C 0 ðtÞ ¼ kðCÞ C 0 ðtÞ; (3)
0
where the flux is first order in C when the background non-tracer
mother substance concentration is constant. When only the tracer is
changing concentration, the k(C) is constant and the system is first
order and linear. In general then, one can look upon tracers in
compartmental systems as being linear, first-order systems, though
nowadays they can go far beyond that. The originators and later
proponents of compartmental analysis (Berman (1); Jacquez (2, 3);
17 Compartmental Modeling in the Analysis of Biological Systems 393
Cobelli et al. (4)) used this simplification, but were always aware of
the greater possibilities of allowing nonlinear coefficients. Berman’s
classic 1963 article (1) provides much more than solutions to ordi-
nary differential equations (ODEs) for he outlines an important
philosophic approach to modeling in general. Jacquez’ books and
many articles, and the book by Cobelli et al. (4), give detailed
mathematical approaches and explicit applications. The desire to
use linear kinetics was not so much to avoid solving nonlinear
equations as it was to use linear algebra to solve the differential
equations. A system of linear differential equations can be solved
by matrix inversion and can provide the much desired analytical
solutions. As we shall see below, analytical solutions are still desired,
for they serve as verification that the numerical solutions produced
by modern simulation systems are correct in specific reduced cases,
and thereby imply that the nonlinear system solutions in that neigh-
borhood of parameter space are also correct. But, because most
biological phenomena are nonlinear, such that the rate coefficients
vary with the concentrations of one or usually more solutes, tem-
perature, and pH, we have to acknowledge right at the start that
using linear compartmental systems analysis is an approximation.
Compartmental analysis was mostly descriptive, not mechanis-
tic. It was “Black Box,” not attempting to define enzymatic reac-
tions mechanistically, but to describe the time course, “White Box”
modeling is where the innards of the operational analysis attempt to
describe mechanism, not just the kinetics of a relationship. Never-
theless, the descriptive level was a success; in FDA reviews quanti-
tative descriptions are valuable, for they distinguish groups of
responses and allow categorization even when they cannot provide
a physiological interpretation. The plasma concentration–time data
are very useful in choosing methods of administration and in defin-
ing dosage regiments.
Modern molecular biology and emerging integrative multi-scale
modeling analysis likewise are changing the game. Personalized
medicine is pointedly mechanistic, with cell and molecular physiol-
ogy dominating in the strategies of Administration, Distribution,
Metabolism, and Elimination (ADME). Fortunately the huge
increase in the rates of acquisition of data, causing a demand for
detailed, informative but complex simulation analysis has been more
than compensated for by the increases in computational speed and in
improved software facilitating modeling analysis. Most importantly,
software sharing is now relatively easy, and of much increased impor-
tance since comprehensive models may take years in development.
Now, highly nonlinear complex systems, spatially distributed or
lumped, are handled with faster computation, and can include kinet-
ics and detailed physiological pharmacodynamics (the PD of
PKPD), which is the systematic analysis of the body’s responses to
the drug or toxin. Given the relevance of physiological transport
processes (diffusion, flow, transmembrane exchange, binding) in
both administration and distribution, a relatively new term, PBPK,

physiologically based pharmacokinetics, has arisen to recognize the
importance of incorporating anatomy and physiology into PK.
Drugs are toxins. It is only a matter of dosage. The struggle to
distinguish acceptable toxicity from unacceptable toxicity is the
central conflict, not just for cancer therapy but also in defining
drug usage in general. Aspirin, ibuprofen, oxycodone, sugar, and
water, all create problems when in excess. Distinguishing “Thera-
peutic Dose” from “Toxic Dose” depends on the drug and on the
particulars of the patient (size, age, body fat level, other drugs,
physiological state and past history, genetic heritage). Many sub-
stances in our environment augment the difficulties, adding other
sources of specific or general toxicity.
Drugs and toxins have many common features; for example the
lipid solubility that allows easy permeation of cell membranes, so
desirable in drugs, is the source of the problems with inhaled
toxicants. While the body has evolved a rather general system for
dealing with foreign toxicants, the P450 system in the liver that
handles hundreds of different chemicals, it is also good at degrad-
ing and excreting drugs as well. As a corollary, hepatotoxicity can be
a problem with drugs. Likewise, renal damage can be a risk from
those drugs excreted in the urine, like ibuprofen.
Computer modeling includes all phases of ADME. The
method of Administration by swallowing a pill is the commonest
but is not as fast as intravenous (i.v.) injection. Intravenous is the
method with the best defined administrative kinetics, followed by
intramuscular (i.m.). There are a host of other local injection types,
slow release subdermally, suppositories, inhalation, sublingual
absorbance, etc., all of which have different rates of drug delivery
into the circulation and to the target. Since the exposure of the
target to the drug is most often measured in terms of the AUC
(Area Under the Curve of the drug’s plasma concentration versus
time) the differences amongst methods of delivery are important.
The AUC is influenced by dilution in the circulation (part of
Distribution), degradation by hydrolysis or metabolism (the M of
ADME), and elimination or excretion (the E of ADME), so a
pharmacokinetic model must include all of these. See earlier chap-
ters on Fundamentals of ADME and on Modeling of Absorption.
In the neighborhood of the target (enzyme, receptor, channel,
transporter, or transcription regulator) much depends on physico-
chemical attributes of the drug or toxin. Does it bind to plasma
proteins? What is the affinity of the target proteins compared to that
of other competing proteins, or DNA or membrane-bound pro-
teins? What are the on- and off-rates of any competing binding
sites? What is the drug’s tissue/blood partition coefficient or its
solubility in body fat? One of the standard ways of getting clues on
the fate of the drug is to do whole body distribution studies on rats
to see where it is deposited at a succession of points in time. The
distribution sites of positron-labeled drug may be shown with high

resolution in reconstructed 3D imaging using MicroPET (Positron
Emission Tomography) systems. Such observations give data on
Distribution and Excretion and sometimes Metabolism. Since the
retention times in specific tissue locations influence, and may dom-
inate, the AUC, this information is crucial for optimizing effect on
the target while minimizing toxic side effects. Modeling analysis
then ideally provides specific information on all of these aspects of
the kinetics, and further provides the information essential for a
critical understanding of the pharmacodynamics. A by-product of a
good understanding of the PK is that it allows one to consider using
combined drug therapy wherein a second drug is used in advance to
protect a sensitive binding site by binding to it harmlessly and
preventing the toxicant drug from binding. Variants on this
theme might be combining therapy with a drug preventing the
activation of a receptor whose binding site accepts the drug of
interest.
1.2. Model Types, 1. Closed system: No sinks or sources, literally, all fluxes between
Topologies, and inside and outside are zero, and external driving forces are all
Equation Types zero.
1.2.1. In Terms of Input 2. Open system: There are external sinks and/or sources for some
and Output Characteristics of the constituents in some of the compartments or cells. (Sinks
We May Classify are defined as operations via which a substance vanishes from
Compartmental Systems the system; sources are operations generating a substance.)
as Closed or Open
1.2.2. In Terms 1. Catenary system: Two or more compartments arranged in

of Interconnections Amongst series.
Compartments or Cells, the 2. Cyclic system: Three or more in series, with the last connected to
Topology of the Network Is the first allowing a circulating flux.
Useful in Defining
Mathematical or Analytical
3. Mammillary: Two or more peripheral compartments connected
Approaches
to a central compartment and having no cyclic components, e.g.,
a blood compartment connected to each organ in the body.
Equating the circulatory system, for example, to a mammillary
compartmental system raises questions. “How can this be a rational
description?” Total circulatory mixing in humans requires many
minutes. Is the rate of solute escape from blood so slow that mixing
throughout the whole circulation is fast in comparison to the rate of
exchange within the organs? Solute-binding to plasma proteins
might so retard escape from the blood that the approximation is
adequate. Alternatively, consider the mouse, where the circulatory
mixing time is a few seconds, but permeabilities are similar to those
in humans, so equilibration in tissues is slow compared to circula-
tion times; here the idea that the system is mammillary is more
reasonable. The basic compartmental premises, instantaneous
mixing throughout the compartment, and therefore uniform inter-

nal compartmental concentrations are almost never truly valid;
being alert to this implies the next step: evaluating the error due
to failure to fulfill the requirements.
In terms of types of equations, models may be expressed in
many forms, but physiological models set up for solving by numer-
ical methods are mainly in the form of ODEs, Differential Alge-
braic Equations (DAEs), and partial differential equations (PDEs).
First-order ODEs are central to compartmental modeling, and
may be linear and nonlinear, for example a Michaelis–Menten
equation. A set of first-order equations set up with a single variable
to the left of each equal sign is said to be in state variable form. See
the Background chapters on Linear Algebra and ODEs. ODEs and
DAEs may be mixed together in a model, where the DAEs define
variables used in the ODEs, either implicitly or explicitly. Seeking
solutions sometimes requires prior analysis, but solvers like JSim
and Matlab can handle many implicit forms. PDEs are becoming
more common in PKPD problems now that computation is fast.
Especially recently it is being recognized that there are usually
concentration gradients along the lengths of capillaries for con-
sumed substrates and for drugs during the uptake phase; the
existence of gradients violates the compartmental assumption of
uniform concentration (5, 6) and causes errors in the estimation of
permeabilities and conductance parameters of all sorts. One-
dimensional PDEs usually suffice for capillary–tissue exchange
since in well-perfused organs the radial distances between blood
and cells are so short that radial diffusional retardation is negligible
compared to membrane permeation across endothelial cells or
parenchymal cell membranes.
1.3. Distribution Distribution, the D of ADME, is by convection, diffusion, and

from the Site permeation, and must precede the drug’s therapeutic and degradative
of Administration reactions at target sites or metabolic sites. Distribution includes the
reversible transfer of a drug between one compartment and another.
Some factors affecting drug distribution include regional blood flow
rates, volumes of interstitial and cellular spaces, molecular size, polar-
ity, and binding to serum proteins or other nontarget sites, forming
complexes. In using the AUC of plasma concentration versus time,
one should keep in mind that plasma protein-bound drug is having
no effect at the targeted site. Administration and Distribution are
best treated together as it is the combination of processes that
governs the effective concentration reaching the target.
Whatever the administration method (delivery on nanoparticles,
laser-induced controlled release, ultrasound enhancement of intra-
dermal diffusion, pill decomposition rates, location in gut, etc.), the
next phases of Distribution are physiological processes that are fairly
well understood (convection, permeation, diffusion, transport across
membranes or intracellular microtubular transport), but not
necessarily well characterized for the particular drug. Heterogeneities

in regional flows, capillary densities, and tissue composition may have
to be accounted for. These precede the reaction, binding, inhibition,
or receptor-mediated responses that compose the desired pathophys-
iological responses. Metabolic reactions, sequestration, uptake, and
excretion by epithelial cells (liver, kidney, saliva, skin, etc.), as by the
liver’s P450 system, or other reactions which inactivate (glucuronida-
tion, glycosylation) are all a part of the ADME model and have to be
considered in any PK modeling. For compounds which are degraded,
the metabolic products have to be assessed for long-term effects.
Compounds excreted by the liver or kidney become concentrated,
even 100-fold, in the process of elimination in urine or bile, and if the
drug is not conjugated or inactivated in some way there may be
damage to the excretory organ, compromising the organ function
and changing the pharmacokinetics after prolonged usage, and rais-
ing the AUC following each administration. Both renal glomerular
filtration and tubular secretion and hepatic biliary excretion create
high concentrations of the drug or metabolites.
At the level of the target, Quantitative structure–activity rela-
tionship (QSAR) (sometimes quantitative structure–property rela-
tionship: QSPR) comes into play. This is the process by which
chemical structure is quantitatively correlated with biological or
chemical reactivity. Biological activity can be expressed quantita-
tively as in the concentration of a substance required to give a
certain biological response, but in the context of quantitative PK
analysis, one would prefer that the PD (pharmacodynamics or
biological response) also be expressed in mechanistic terms. Addi-
tionally, when the physicochemical properties are expressible in
numbers, one can formulate the quantitative structure-activity rela-
tionships in the form of a mathematical model. The mathematical
model may also predict the biological response to related chemical
compounds. But it is the concentration of the bound agent com-
plexed to the target site that determines the level and duration of
the response.
There are a good many useful measures in considering the PK
of ADME. These include, in addition to AUC (exposure), Cmax
(maximum concentration), Tmax (Time to Cmax), half-life, clearance
route and rate, volume of distribution, and bioavailability in
unbound form. In the steady state of long-term administration
one assesses the plasma concentrations, total accumulation, linear
or nonlinear PK, time-dependent changes in kinetics, and metabo-
lites, their identity, and their PK. It is the combination of PK and
PD that allows knowledgeable optimization of dosage regimens.
A well-developed PKPD model will account for most of these
considerations and thereby be predictive and advisory to the thera-
pist. A clear description of the kinetics allows the planning of
efficient dosage regimens.
2. Software
for Systems
Description,
Simulation, Pharmacokinetic models have been written in virtually every
and Data Analysis computer language that exists, and it is a field that has stimulated
the development of a large set of relatively specialized simulation
systems. A partial list of simulation software for compartmental
2.1. Common Software
analysis goes back half a century:
and Methods Used
in the Field SAAM, Simulation, and Analysis Modeling, was the first, developed by
Mones Berman at NIH for analyzing tracer kinetics (http://depts.
washington.edu/saam2/). It exists still as SAAMII.
SIMCON, a general simulation control system (7), now evolved
into JSim, was used to solve FORTRAN-based models of all
sorts.
XPPAUT, from Bard Ermentrout (http://www.math.pitt.edu/
~bard/xpp/xpp.html). XPPAUT is particularly good for bifur-
cation analysis of dynamical systems.
Gepasi, now CoPasi, from Pedro Mendez (http://www.softpedia.
com/progDownload/Gepasi-Download-167140.html;
http://www.copasi.org/tiki-view_articles.php). CoPasi is espe-
cially good for enzymatic reactions and biochemical systems,
allowing a menu of choices for reaction types.
Modelica, http://www.openmodelica.org/, is excellent for linking
operators and presenting the forms of model networks (http://
www.ida.liu.se/labs/pelab/modelica/OpenSourceModelica-
Consortium.html).
Jarnac is designed for symbolic or diagrammatic entry for biochem-
ical and gene regulatory reactions (from Herbert Sauro (8, 9),
http://sys-bio.org/jarnac/).
JSim: Developed from SIMCON and XSim (for X-windows linux
systems, http://www.physiome.org/software/xsim/) into
JSim (10) (http://www.physiome.org/jsim/). JSim was devel-
oped by Erik Butterworth (11); it provides automated unit
balance checking and unit conversion, thus avoiding errors
due to inconsistency in the units used in the code.
Non-MEM, nonlinear mixed effects modeling, is a commercial
software package providing the capability to use a wide variety
of pharmacokinetic models. It is particularly designed for the
analysis of sparse data sets using combinations of single patient
and population data, (http://www.iconplc.com/nonmem).
BioSPICE (derived from SPICE, for biology: http://sourceforge.
net/projects/biospice) is designed for molecular biology. See
also http://jigcell.cs.vt.edu/software.php for JigCell, based on
BioSpice.
Cellular Open Resource (COR) is for a Windows environment:

http://cor.physiol.ox.ac.uk/ for cellular level physiological
systems. PCEnv for physiological systems is being developed
from it.
Stella (http://www.iseesystems.com, a commercial system) for
networks of operators such as in compartmental systems.
StochSim (http://www.ebi.ac.uk/~lenov/stochsim.html) is writ-
ten explicitly for the treatment of molecular interaction when
there are few molecules and interactions occur stochastically.
These simulation systems have one purpose in common: to
make the programming of models simpler and to facilitate the
analysis of experimental data in terms of the parameterized descrip-
tions of the kinetics. They vary considerably with respect to their
representation of physical/chemical mechanisms. Such modeling
and analysis systems do not displace FORTRAN, C, C++, and Java
as mainstream languages, but rather they replace the front-end
entry to formulate the models and interface them to data sets.
None of these have the general capabilities of Matlab or Mathema-
tica, nor do they attempt algorithmic manipulation as in Maple, but
are more directly tuned to the user’s needs, as will be described.
2.2. A Preferred JSim is perhaps the most general of these simulation analysis systems,
Simulation designed for the analysis of experimental data. It is built around a
System, JSim “project file, .proj,” that may hold many data sets, several different
models, and results of multiple analyses. JSim handles not only the
ODEs around which traditional compartmental modeling is built,
but also DAEs, implicit functions, PDEs, and stochastic equations.
JSim, uniquely, and from its beginning in 1999, uses unit balance
checking and automated unit conversion. (Unit balance checking
assures that the units of the expressions on the left of the equal sign
are the same as those on the right. Automated unit conversion means
that when time is expressed in minutes a velocity expressed in cm/
s will be converted to cm/min by multiplying cm/s by 60 s/min.)
This pair of features is a great boon in programming since in the first
phase of compilation it automates the first stage of verification of the
model’s mathematical implementation by making sure that every
equation has unitary balance. The second phase of compilation parses
the details of the equations, and sequences them for efficient compu-
tation. The run-time code is compiled into Java, which now runs
almost as fast as FORTRAN and C. (On a cardiovascular–respiratory
system model JSim ran exactly 300 times faster than a Matlab–Simu-
link version of the identical model.) JSim’s advantages over the ODE-
based systems listed above are the following:
1. Runs on Linux, Macintosh, and Windows.
2. Is free and downloadable from www.physiome.org. On the
Macintosh it takes about 30 s to download and install, and
another 10 s to bring up a model.
3. Is the only one that solves PDEs and offers an assortment of solvers
for both PDEs (three available now) and ODEs (eight available).
4. Imports and exports both SBML and CellML archival forms.
5. Provides sensitivity analysis of two types, relative and absolute.
6. Graphical output is immediately available during the simulation
run and setup in seconds.
7. Has seven built-in optimizers for excellent power in parameter
adjustment to fit data.
8. Provides the covariance matrix giving the correlation among
free parameters and estimates of parameter confidence limits.
9. Use project files that allow the analysis of many experimental
data sets in one file.
10. Stores parameter sets so that individualized parameter sets for
each data set can be stored.
11. Allows the use of several models within one project file so that
competing hypotheses (models) can be compared and evaluated.
12. Is structured so that the front-end parameter control and
graphical user interface (GUI) can be framed explicitly for any
model.
13. Has linear and log line graphs, 2D contour plots representing
3D, and phase-plane plots.
14. Has “looping” capability, allowing discrete successive jumps of
the values of one or two parameters at a time in order to
explore system sensitivities visually and rapidly.
15. Uses a Mathematical Modeling Language, MML, in which one
writes the equations directly, for simultaneous solution, and in
which the order of the equations is not specified.
There are no special requirements for the JSim software or for
its methods of use for model building and exploring or use in
analysis with respect to hardware, computing platform, or operating
system. It has important limitations, not being a procedural lan-
guage but a declarative mathematical language. This means there is
no equivalent of a FORTRAN DO-loop (or GO TOs or jumps). It
cannot yet do matrix inversions (except through a special mecha-
nism), and is in a continuing state of development. JSim 2.0,
released in February 2011, is based on a new compiler providing
many new features described at nsr.bioeng.washington.edu/JSim/.
The features listed above, and others not listed, have been
implemented in JSim because the years of experience with a large
variety of models, with teaching graduate and undergraduate clas-
ses, and postdoctoral and faculty workshops have led to a detailed
understanding of how people use modeling in scientific research.
Experiment design, hypothesis testing, and system parameteriza-
tion are given priority in the conveniences provided.
Thus JSim, since it is designed around the analysis of experimen-

tal data, is our preferred software and will be used for the compart-
mental modeling shown next.
3. Compartmental
Modeling
3.1. The Modeling The overall process in the experiment/model hypothesis iteration
Process loop of Platt (12) is as follows: (1) express the hypothesis in
quantitative terms, as a mathematical model, with units on every-
thing; (2) use the model to determine the best experiment that
might contradict the predictions of the model, or, better yet,
develop an alternative model that is seemingly as good but makes
different predictions, and then design the experiment that clearly
distinguishes between the models; and (3) do the experiment and
analyze the data. One of the two competing models, maybe both,
must be proven wrong, and so science is advanced.
The normal data analysis using models begins by putting the
data to be analyzed in a “project file,” modelname.proj, and dis-
playing them on the JSim plot-pages. The second stage, coding and
model verification in accord with standards (http://www.imagwiki.
nibib.nih.gov/mediawiki/index.php?title¼Working_Group_10),
is building and testing the model, incorporating reference analytical
solutions if appropriate to verify the solutions as being mathemati-
cally accurate, and representing the equations. The “project file”
may contain two or more models so that alternative model forms
can be compared directly by examining the solutions (changing
parameters and rerunning, using “loops” to automatically change
parameters, using behavioral analysis, plotting in various forms
including phase-plane plots, contour plots). The verification stage
is to show that the model solutions are computed correctly, done by
testing different solvers, using different time-step sizes, and com-
paring with analytical solutions in special cases.
The validation stage is to test the fitting of the model solutions
against the experimental data. The word “validation” is truly opti-
mistic, because a good fit of the model solution to the data does not
really validate the model, but merely fails to invalidate it. It is the
failure of the model that leads to the scientific advances by forcing
new ideas to be incorporated. Nevertheless, fitting of the model to
the data provides characterization of the data, augments diagnostic
acuity, assesses progress of disease or evidence of successful therapy,
and is generally useful in reconciling the working hypothesis
(model) with observations.
The final phase is preparing the model so that it can be reproduced
by others. This is not only critical from a tutorial point of view but also
in fact is a requirement for any scientific publication. Anything that
cannot be reproduced is misleading and wastes time and money.

Reproducible models can be tested by others or used as building
blocks to advance the field when they pass muster.
3.2. A Simple The modeling code. For an introduction to JSim we use a two-
Compartmental Model compartment closed system with passive exchange between the
Implemented in JSim compartments, and a conversion reaction of solute A to solute B
in either or both compartments. This model has analytic solutions
which could be used either to show the solutions or to provide
verification of the accuracy of the numerical solution, but since
these solutions run no faster than the numerical solutions, they
will not be used here. Detailed instruction in JSim use is available at
http://www.physiome.org/jsim/. This model is #246.
Many model programs are available at http://www.physiome.
org/Models. One can search to find a model similar to what one
might like to construct, e.g., from a tutorial list of compartmental
models: http://www.physiome.org/jsim/models/webmodel/NSR/
TUTORIAL/COMPARTMENTAL/index.html.
Open model #246, Comp2ExchReact, and you will be asked to
allow the display on your computer, wait a moment for it to be
compiled, and then click on “Source” at the bottom of the JSim
page to show the source code, Table 1 for the model shown in
Fig. 1. All the models on the Web site are archived to keep track of
model changes (previous versions can be found under the “Model
History” section on each model Web page). (Models edited over
the Web cannot be saved on the Physiome Web server, so simply
save what you want to your own directory. The JSim system can
likewise be downloaded directly and the model worked on from
your own computer.) The model for the code in Fig. 1 is dia-
grammed at the bottom of the figure; it is a two-compartment
model for two substances, A and B. Both substances can passively
move from one compartment to the other. A is irreversibly con-
verted to B in either or both compartments. After the title a short
description is provided. (Text enclosed by /* . . .*/ is ignored by
the compiler, as is comment text following // on any line.)
An important JSim feature is invoked next, the reading of a
units file, nsrunit, which allows automated unit balance checking
and automated unit conversion to common units during the
compilation phase; MKS, CGS, and English units can all be
used. (The unit conversion can be turned off, a feature useful
Fig. 1. JSim code for a two-compartment model in which a solute A can exchange across
a membrane between volume V1 and V2 and can react irreversibly to form solute B in
either compartment. This is available for download at www.physiome.org/jsim/models/
webmodel/NSR/Comp2ExchangeReaction/index.html and is model #246.
Table 1
Code for two-compartmental model with passive bidirectional exchange
when importing models from CellML or SBML which sometimes

have unit conversion factors hidden as dimensionless factors in their
archived models, whereas they should be dimensioned, e.g., as 60 s/
min.) The word “math” and the curly bracket designate the start of
the model code, written in JSim’s MML, which will be seen to be
merely the equations for the model. The “realDomain” defines the
independent variable as t, time in seconds, along with a starting and
ending time and a time interval for graphing the solutions.
The parameters of the model are assigned units. We have put
the units in physiological form in order to represent those for a
perfused tissue, so much per gram of tissue. In general it is practical
to avoid mistakes by using the same units as used for the
experimental studies in which the data are acquired. The volumes of

the two compartments are V1 and V2 ml/g. In this context, flow
per gram of tissue (F), clearances (Cl), and permeability-surface
area products (PS) all have the same units, ml/(g min). (The PS is
defined as permeability, P cm/min, times membrane surface area
per gram of tissue, S cm2/g.) The reaction is described as a clear-
ance, but here is given the symbol G, ml/(g min) for a gulosity or
consumption. The solute A is converted from A to B. (The termi-
nology used here is that used by the American Physiological Society
(13), where an extensive set of terms for transport, flow, and
electrophysiology are given.)
The model’s variable are functions of the independent variable
time, for example the concentration of A is defined as a real number
and as a function of time by “real A(t).” In the equations A(t) can
be written without the “(t),” but in JSim’s MML the use of the (t)
initially is required to establish it as a function of time. Other
languages such as Madonna and XPPaut do not demand this, but
JSim users find that this reduces errors. The MML is fairly similar to
that in those languages; the basic intention is just to write the
equations directly.
The ODEs, used to describe compartmental models, need to
be supplied with initial conditions, here provided as A10, A20,
B10, and B20 where the first subscript digit refers the compartment
and the second is “0” to refer to the initial time t ¼ t min, which
can be negative or positive, but is usually t ¼ 0.
In the ODEs, the derivative of A(t), which you would expect to
be dA/dt is written A:t. The second derivative, d2A/dt2, would be
A:t:t. The equations are written in state variable form in this model,
for tidiness and simplicity, but this is not FORTRAN and one could
equally well write the equation as “V1 A1:t ¼ A2 PSa A1
(PSa + G1);”. In this latter form, the left-hand side, LHS, of the
equation is the rate of change of mass of solute A, (mol/g)/min ¼
mM (ml/(g min)) for A1 in V1.
3.2.1. Unit Balance Checking The acute observer will have noticed that the independent variable, t,
is in seconds: the differential equation therefore looks to have unbal-
anced units. It actually computes correctly since the phrase “unit
conversion on” at the top of the program preceding the model code
enlists the automated unit conversions so that in the compiled code
the multiplier of the left side, “60 s/min,” is inserted. Without unit
conversion on, the best way to handle this is to put the independent
variable t in minutes. In languages like Matlab there is no unit
checking. In huge projects like the Mars Climate Orbiter mission
(http://www.spaceref.com/news/viewpr.html?pid¼2937), mixing
units from the European and American programs led to the crash
of the space vehicle and the termination of the billion dollar mission.
There would have been no problem in JSim as long as the units were
stated and unit conversion “on” (11). The nsrunit file may be viewed
by clicking “Debug” (left, bottom) to drop down a menu including

“View system units file.” The four ODEs in Fig. 1 are a complete
description of the system behavior after the initial moment. Follow-
ing them is an algebraic equation for TotalC, also a function of time,
to check the total mass at each point in time divided by the total
volume and giving the average concentration. With the initial con-
ditions and volumes given, the result is 0.33333333 mM for every
time step. Run the model and then check the numbers for the
plotted variables by clicking on “Text” instead of “Graph” at the
bottom of the plot page, as in Fig. 2. The program end is marked by
a right curly bracket, beyond which one can put notes, comments,
key words, diagrams, references, etc.
3.2.2. Graphical User The JSim GUI for simulation control is shown in Fig. 2. The left
Interface panel has overlays for project contents: one or more models, data
sets, parameters sets, setups for solvers for ODEs and PDEs,
sensitivity analysis, optimization, and confidence limits. To start a
simulation run one clicks “RUN” at the top of the run time
window. On the right page one clicks on a “Message” panel for
error messages, or the plot pages (1_Conc or 2_ReacSite, names
chosen by the user), and at the page bottom choose “Graph” for
displaying the output graphically or “Text” for seeing the numerical
listings of the experimental data and the model solutions at each
time point.
There is an extensive introduction to JSim at www.physiome.
org/jsim giving precise detail to supplement this outline of usage.
The JSim MML code is pretty easy for a beginner to use since
it contains just the parameters with their units, the variables
followed by (t) to indicate their dependence on the independent
variable t for time, the initial conditions, and the equations for the
model. The most common mistakes are to misapply the uses of
commas and semicolons. Each of a sequence of events is usually
comma separated, and a string of them is closed with a semicolon.
The model code itself begins with the left curly bracket, “{” after
the word “math” and ends with the right curly bracket,“},”
after all of the equations have been written. Comments are pre-
ceded by a double slash, “//,” or alternatively can be preceded by
a “/*” followed by an “*/,” without the quote signs. Equations
end with a semicolon, including those in the initial condition
statement.
3.2.3. Exploring Parameter In loop mode, the user can choose to enter a sequence of values
Influences Using the “LOOP” (under Other Values) to explore model behavior widely, using
Mode of Operation comma separation, e.g., 2, 3, 5, 8, etc., and in the “auto” mode
will do as many runs as there are values entered. One can also enter
arithmetic changes such as @*2 or @ + 3 or more complicated
expressions to indicate automatic changes in the starting value by
Fig. 2. Standard JSim Input/Output control and plot pages. Left panel: Runtime control: “Domain” is time, t, with tmin
starting time, tmax ending, and tdelta the time step for plotting. “Model Inputs” gives parameter values and initial conditions.
Model Outputs shows values of variables at the end of the run, at t ¼ 60 s. The mass balance check, TotalC, is exact to
eight decimals. Right panel: The time courses of concentrations and A and B in compartments 1 and 2 are shown for the
parameters and initial conditions shown in the left panel and in the code in Fig. 1. The user chooses the variables to plot,
and the colors and line type or point type. The title and labeling are user written and retained in the JSim project file. The
legend for the graphics output is automated.
multiplication by 2 on each run, or the automatic addition of 3 on

each run, for a chosen number of runs (Fig. 3).
3.2.4. Sensitivity Analysis Sensitivity analysis is available to provide a quantitative measure of

the effect of any chosen parameter on a particular variable. This
extends the information gained by “looping.” The sensitivity func-
tion S(t) of a variable to a parameter is calculated by calculating the
change in a variable such as A(t) per fractional change in a parame-
ter value, P. The linear sensitivity function is
SðtÞ ¼ @AðtÞ=@P; (4)
or alternatively, the log sensitivity function is
Slog ðtÞ ¼ ð@AðtÞ=AðtÞÞ=ð@P=PÞ ¼ logðAðtÞÞ=@ log P: (5)
Fig. 3. Loop Mode Operation. The control panel for the looping operator is shown on the left. The starting values for G1 and
G2 are those shown in the code (Fig. 1) and their solutions are given by the solid lines, and are the same as in Fig. 2, with
conversion of A ! B only in compartment 1. On the second run (dashed lines of same colors) the values entered by the
user under Other Values (left panel, top right ) are used, in this case setting G1 ¼ 2 and G2 ¼ 0, so that the conversion of
A ! B occurs now only in compartment 1 instead of in compartment 2.
Sensitivities are local linear approximations, and are calculated

by computing two solutions, A1(t) and A2(t), the second having P
changed by a small fraction, 0.001P or 0.01P, so that
SðtÞ ¼ ðA1 ðtÞ A2 ðtÞÞ=DP; (6)
and using the definition of the derivative as DP goes to zero gives
the formal value.
3.2.5. Optimization in Data Optimization, either by manual parameter adjustment or by

Fitting and Analysis automated methods, is the procedure of adjusting parameters of the
equations so that the solution provides a good fit to the data. The
evaluation of the goodness of fit by minimizing the sum of squares of
the differences between the model solution and the experimental data
is based on an implicit assumption that the differences are Gaussian
random. This is seldom correct, and it is important to appreciate that
the choice of the distance function is personal, i.e., it is up to the
investigator to characterize the noise in the data and to weight the
influence of individual data points. See next section.
4. Applications
4.1. A Compartmental Aspirin is a very old drug used to reduce fever and inflammation;
Approach to Aspirin only recently has its mechanism of action begun to be understood,
Kinetics the first being its action as a blocker of prostaglandin formation. Its
kinetics have not been thoroughly worked out, so we present here
our analysis of three sets of representative data from three different
studies. Aspirin is acetylsalicylic acid, and its reactions are as follows:
Acetylsalicylic acid ! salicylic acid ! salicyluric acid: (7)
Aspirin, acetylsalicylic acid, is hydrolyzed quickly via a plasma
esterase to salicylate. Salicylate is pharmacologically active. Salicy-
late kinetics dominate the clearance. The modeling examines only
salicylate’s enzymatic conversion to product, where product is con-
sidered to be equivalent to excretion into the urine, a saturable
process. This may produce salicyluric acid or a glucuronate. The
model captures the kinetics of salicylate clearance over a 100-fold
range of concentrations through consideration of one enzymatic
reaction with parameters optimized to fit three very different data
sets taken from the referenced papers.
The second reaction to form the excretable product is slow and
is enzymatically facilitated. At high doses the enzyme becomes
saturated, i.e., the reaction is limited by the fact that the enzyme
is all in the form of the bound enzyme–substrate complex and
raising the salicylate concentration does not accelerate the reaction.
At low dosage the clearance is rapid; at medium or high therapeutic
dosage the clearances are slower; at near-lethal toxic levels clearance
is very slow and often requires treatment by infusion of alkaline salts
and sometimes dialysis. For this example we have taken the data
from three research studies. (Low dose data are from Fig. 1 Right of
Benedek (14) Dose Period 1 (squares). Medium dose data are from
Fig. 4 of Aarons (15) oral dosage, last nine points. High dose data
are from Fig. 1 of Prescott (16) Control.) These particular data
were chosen because the chemical methods and procedures
appeared to be excellent, the data covered many hours, and, while
we do not have the original data, conversion from the symbols in
the figures to numerical representation was accomplished with
good accuracy.
A reason for choosing aspirin as the subject for compartmental
analysis is that the time course of the clearance, mainly by loss into the
urine, is long compared to circulatory mixing times, so that the bio-
chemical processes appeared to limit the clearance, and they could
therefore be characterized. If the circulation and distribution times
were long compared to the reaction processes, the latter would not be
meaningfully determined from the observations. Another reason was
that a comparison among the different data sets suggested that the
clearance was enzymatically mediated: at high concentration the dimi-

nution was almost linear, at low concentrations it was almost exponen-
tial, and at intermediate concentrations the rate of clearance appeared
to speed up with time. This fits the expectations for a saturable enzy-
matic process: at low concentrations well below the dissociation con-
stant for the enzyme substrate complex almost all of the enzyme is free
and available for the reaction and therefore the fractional transforma-
tion of substrate is at its highest. With all the enzyme free this is a first-
order process, a single exponential. In contrast at high concentrations
the enzyme is almost totally saturated and the conversion is at a
maximum rate independent of the concentration; this is zero order
kinetics, giving a linear diminution in concentration. Thus the shift
from zero order kinetics to what occurs at intermediate concentration,
a gradually increasing fractional rate of reaction, is what was suspected
by looking at the middle-level concentrations.
Thus we hypothesize an enzyme conversion model for salicylic
acid clearance, and then test the surmise by attempting to simulta-
neously fit the three independent data sets using one set of para-
meters. To do this in one program and to optimize the fit to the
data for all three salicylate levels we coded three identical models in
the one program (Table 2). These are computed simultaneously in
order to fit the three independent data sets from the three research
reports using one common set of parameters, and automated optimi-
zation was used to minimize the set of differences between data and
model solutions.
The reactions, both the binding to the enzyme and the product
formation, are reversible:
kon1 ! koff 2 !
SA þ E ! SAE ! E þ P: (8)
koff 1 kon2
The equations and parameters are identical for the three dosage
levels. Here is the model code for the low dose, where the prefix L
distinguishes this model equations from those for the medium
dose, prefixed M and the high dose, prefixed H, neither of which
are shown in Table 2.
In undertaking an analysis on a single enzymatic reaction we
lack knowledge of the exact mechanism. In addition the compart-
mental approximation is certainly questionable for whole body
studies. The hypothesis that a single enzymatic reaction dominates
the clearance would be strengthened if the model provides good fits
to the three data sets. There is no guidance from the literature on
the dissociation constant, KD, for our presumed enzymatic reac-
tion, so that we are neither constrained nor aided.
We chose to use a simple enzymatic reaction, one that allowed
characterizing the rates of binding and unbinding of substrate and
enzyme, and a rate of the forward reaction to yield the product. We
allow also a backward reaction, on the basis that all reactions are
Table 2
Model code for salicylate clearance (model downloadable
from www.physiome.org/Models: search for model 280)
/* MODEL NUMBER 280

MODEL NAME: Aspirin
SHORT DESCRIPTION: Salicylic acid (SA) clearance for three different dose ranges is modeled as
an enzyme reaction. This table is abbreviated by omitting the code for the Mid and High dose
reactions.
*/
import nsrunit; unit conversion on;
math Aspirin {
// INDEPENDENT VARIABLE
realDomain t hour; t.min ¼ 0; t.max ¼ 16.0; t.delta ¼ 0.05;
// PARAMETERS (SAME FOR ALL THREE MODELS)
real kon1 ¼ 0.174 L*mg^(-1)/hour, // On rate for SA + enzyme
KD1 ¼ 6.3 mg/L, // Dissociation constant for SA enzyme complex
koff1 ¼ KD1*kon1, // Off rate for SA enzyme complex
kon2 ¼ 0.003 L/(mg*hour), // On rate for Product + enzyme
KD2 ¼ 250 mg/L, // Dissociation constant for Product enzyme complex
koff2 ¼ KD2*kon2, // Forward rate to form Product from complex
Gp ¼ 0.03 1/hour, // Clearance rate from plasma
Etot ¼ 10 mg/L; // Total enzyme concn
// LOW DOSE MODEL: Data from Benedek (1995) Fig 1 Right Dose Period 1 (squares)
// LOW DOSE PARAMETER
real LSAtot ¼ 8.2 mg/L, // Total Low Dose concentration
// LOW DOSE MODEL VARIABLES
LSA(t) mg/L, // Low dose SA
LSAE(t) mg/L, // Low dose SA-enzyme complex
LE(t) mg/L, // Low dose free enzyme
LP(t) mg/L; // Low dose product
// LOW DOSE INITIAL CONDITIONS
when(t ¼ t.min) {LSA ¼ LSAtot; LSAE ¼ 0; LP ¼ 0;}
// LOW DOSE ORDINARY DIFFERENTIAL AND MASS BALANCE EQUATIONS
LSA:t ¼ -kon1*LSA*LE + koff1*LSAE;
LSAE:t ¼ kon1*LSA*LE-koff1*LSAE-koff2*LSAE + kon2*LE*LP;
LP:t ¼ koff2*LSAE - kon2*LE*LP - Gp*LP;
LE ¼ Etot - LSAE; // LSAtot ¼ LSA + LSAE + LP; for an overall mass balance
accounting
} //End of Low Dose Model. The omitted code for the medium and high dose models is identical
// but L is replaced with M or H. Copy and paste the Low Dose Model into the space preceding the
// right curly bracket, twice, and change the L to M for one copy, and L to H for the other; recompile.
thermodynamically reversible, at least in principle. This reaction

setup allows reduction to the simpler, commonly used Michaelis–-
Menten reaction; this is accomplished by speeding up the binding
and unbinding reactions of salicylate to enzyme and eliminating the
Fig. 4. Data on salicylate clearance from the three laboratories are fitted simultaneously with a single enzyme model to
describe the clearance. Parameters are given in the model code in Table 2. Data are from Benedek (14) in the left panel,
Aaron (15) in the middle panel, and Prescott (16) in the right panel. Note that the concentration ranges are markedly
different. (http://www.physiome.org/jsim/models/webmodel/NSR/Aspirin/ model 280).
reverse reaction from the product back to salicylate. So for the first
level of testing all parameters are considered as open and adjustable,
including the initial values of the concentrations in the system. In
this analysis the system is considered to be a single well-stirred tank,
as if the circulation were instantaneously mixed. This gross over-
simplification also makes the assumption that the product either
goes directly into the urine or to some other location in the body
from which it does not return. Enzymatic conversion in the liver
followed by excretion into the bile would be equivalent kinetically
to conversion to a glucuronate followed by the circulation through
the blood and clearance in the kidney. Ideally one would measure
the urinary excretory rates simultaneously and model both the
plasma clearance rates and the urinary clearance. This was not
done here.
The results are shown in Fig. 4 where the curves at the low,
middle, and high concentration ranges are shown to be fitted
reasonably well by the model. The high concentration data (right
panel) are fitted less well but do illustrate that the slope is
Fig. 5. Semilog plots of the data (symbols) fitted with the model solutions for the salicylic acid (LSA, MSA, and HSA, solids
lines). All the data in both Figure 5 were fitted simultaneously with one parameter set for the enzyme as given in Table 2.
The Product concentrations (dashed line for LP, MP, and HP) are merely predicted product concentrations. Since there are
no data on product concentrations, the assumption that there is no degradation of Product must lead to some overestima-
tion of its influence on the backward reaction.
approximately linear, as expected, whereas at the low concentra-

tions (left panel) the curve is nearly exponential, as expected.
A different view is provided by Fig. 5, a semilog plot of the
three sets of data each fitted by the model. The high concentration
curve (open circles, purple line) appears to be almost linear on this
plot, but the slope is shallow, and judgment based on Fig. 4 is
better. The high dose concentrations are very much higher than the
KD1 for substrate binding, estimated at 6.3 mg/l, so that there is
no doubt that the enzyme was almost saturated. In fact, with the
High Dose the concentrations are above apparent KD2 of 250 mg/l
for product binding, so there is significant reversal of the reaction.
At Mid Dose level (triangles, blue line) the slope diminishes as time
progresses, that is to say the fractional clearance is increasing as the
enzyme becomes less saturated. The Low Dose data (diamonds,
green line) are fitted well and have considerable curvature on the
semilog plot, the slope at late times being much diminished: this
leads to the idea that there is some tendency for retention at the low
concentrations, which could be either due to recirculation from
other parts of the body where the concentrations were initially
Fig. 6. Three sets of linear sensitivity functions versus time. Solid lines are sensitivities to the initial zero-time
concentrations resulting from the doses. Long dashed lines are sensitivities to KD1: the sensitivities are all positive.
Dotted lines are sensitivities to KD2, the dissociation constant for the reverse reaction: the sensitivities are all negative.
higher or due to “product pressure” to form salicylate by the

reverse reaction. The reversibility is governed by the KD2 for the
product, which here was about 250 mg/l. There must be even
more reverse flux with the middle and higher level concentrations,
HP and MP, a feature of product inhibition, and the estimate of this
KD2 is determined from both of these, even though the effect is
most evident from the curvature of the low salicylate dose, LSA.
4.1.1. Sensitivity Analysis In Fig. 6 the linear or absolute sensitivity functions are shown
for initial concentrations and for the two dissociation constants.
The solid lines are sensitivities to the initial zero-time concentra-
tions resulting from the doses; most of the sensitivity is at the
earliest points. With the high dose the fractional clearance is so
low that the high sensitivity extends throughout the 16 h of the
study. For the middle dose, MSAtot, the sensitivity diminishes most
steeply as a function of time at around 10 h when the concentration
is close to KD1, the dissociation constant for substrate binding.
The long dashed lines are the sensitivities to KD1; these are all
positive, meaning that if KD1 were increased (decreasing the affinity
of the enzyme for the substrate SA) the model solutions for all
three doses would be at higher levels and the rate of disappearance
would be diminished. Note that the time of peak sensitivity to KD1
is at early times for the low dose, at 10 h for the middle dose, and at
late times for the high dose. The dotted lines are sensitivities to KD2:
the sensitivities are all negative, meaning that if KD2 were increased
(decreasing the affinity of the enzyme for the product P) the model
solutions would be at lower levels and the rate of disappearance
would be increased because of reduced rates of reverse flux from
product to salicylate.
Technically the sensitivity calculations are set up by a special
mechanism: at the bottom of the left page is a button labeled
“Sensitivity.” Clicking on it takes one to the “Sensitivity Analysis
Configurator.” There in the leftmost column of the configurator
table one types in, or chooses from the drag down menu, the
parameter for which one wants to find the sensitivity. By clicking
on the down arrowheads you bring up the choices. In the setup
provided on the Web site at www.physiome.org, etc. the three start-
ing values for the initial concentrations at t ¼ 0 are listed: LSAtot,
MSAtot, and HSAtot. Next on the list are the dissociation constants
KD1 and KD2.Their current values are automatically displayed under
“value.” The calculations of each S(t) are made on the basis of the
parameter change of 1% set under “delta” at 0.01. The tick marks in
the OK column indicate that the calculation will be made as
described earlier, namely, that the standard solution will be calcu-
lated and then another solution calculated for each of the five para-
meters listed, with this 1% change in parameter value. The S(t) is the
difference in the solution at each time point from the standard
solution divided by the 1% change in the parameter value, Eq. 4.
4.1.2. Optimization This is the process of fitting the model solutions as closely as
possible to the data in order to guide one’s thinking about and
one’s use of the model. When the fit is very close, then one has a
descriptor of the fitted data sets, that is, the model and its parameter
set provide a record of that description. Descriptions of many
different studies, patients studies or experiments, allow compari-
sons and possible classification into categories having specific dis-
tinctions. Descriptive models are useful for diagnosis and possibly
for prognosis or choosing modes of therapy. (If the model
“explains” the data by defining the physical and chemical mechan-
isms, that is even better.)
When the fit is poor, then more exploration is needed. Was
automated optimization used? A typical set of trajectories of param-
eter values during an optimization run using SENSOP is shown in
Fig. 7. The values do not range widely; to assure one that they have
not settled into local minima we also used other optimizers that
search widely, e.g. simulated annealing. Try weighting the data
differently: a simple sum of squares minimization may not be
appropriate; it is almost never appropriate if the data range over one

order of magnitude. For example when fitting a decay process that
is exponential, one can use a reciprocal weighting so that points are
weighted more evenly, or a weighting adjusted to the individual
result such as “1/exp(t/(2 h))” for the low dose data set, and “1/
exp(t/(8 h))” for the middle dose data. These choices of weight-
ing are to be typed, without the quotes, into the appropriate line
under Pwgt (point weighting) in the Data to Match Table on the
Optimizer Configuration Page.
In the same table there is opportunity to even up the weightings
for the three quite different curves by using Curve Weighting, Cwgt.
In this case we weighted Low Dose data with 140 times the weight of
the High Dose Data since the latter were about 140 times higher
concentrations. This evens up the contributions to the sum of
squares. Likewise the Middle dose data were weighted at 14, that
is, about ten times higher than the High Dose data. This combina-
tion of Point weighting within each data curve and Curve weighting
among the three data curves gives each point in the triple data set
about the same weight. This makes the calculation of a root mean
square error provided by the minimization of the sum of squares a
reasonable strategy. This is probably close to what one would do if
using “eyeball best fitting,” examining the fitting to all the points.
Both Figs. 1 and 2 plots are valuable in this regard since they weight
the curves differently from the “eyeball” view. While arguing that the
eyeball fitting is just as valid as that from automated optimization,
the latter has the virtue that the weighting scheme is explicitly stated,
and that the weighted sum of squares can be reproducibly reported
when the weighting scheme is reported also.
There are situations where an optimizer fails to reach a good fit
because the sum of squares has settled in a local minimum. Opti-
mizers like Gridsearch and Simulated Annealing are designed to
cover a wide range in state space so that even the corners get
explored (Fig. 7). Others like SENSOP (17), GGOPT (18), and
NL2SOL (19) make excursions, jumping out of the locale to test
other regions for a better fit. To switch optimizers, simply pick
another one from the drag down menu.
4.2. Multiple Sequential Test of clinical functions evolves in a variety of ways, usually being
Dose Administration: designed long after the function has been clearly understood. Here
Hepatic Function is a counterexample, a case in which the idea of the clinical evalua-
tion from the administration of a drug was evident right at
the beginning. There was a coalescence of features that brought
this about.
For the estimation of cardiac output using the indicator dilu-
tion technique Mayo Clinic’s Earl Wood needed a dye that
absorbed light at a wavelength of 800 nm, the isosbestic point at
which oxyhemoglobin and reduced hemoglobin absorbed equally.
This would allow optical detection and quantitation of the dye
Fig. 7. Optimization: Trajectories of values for parameters being optimized, in this case the three initial concentrations and
the two dissociation constants during 100 trials of fitting the model solutions to the salicylate data. The optimizer used was
SENSOP (17); the staircase nature of the plotted values is that SENSOP reports the previous estimate of the parameter
vector as it calculates each new sensitivity function, one for each parameter being optimized, and then reports and plots
the new value of each parameter.
concentration independent of the blood oxygen level. The dye,

indocyanine green, was found by Fox in Kodak’s repository (20);
it was not toxic; its absorbance peaked at 805 nm; it bound to
albumin and so stayed in the circulation and did not color the body.
The densitometer measuring the absorbance was developed (21),
and a spin-off from it was an earpiece densitometer that could be
used to detect blood concentration noninvasively. In early experi-
ments to test the accuracy and reproducibility of the estimates of
cardiac output, dye was injected repeatedly as a bolus into a vein
and the arterial concentration–time curve recorded each time. The
repeated injections led to a rise in the background arterial concen-
tration, raising the question of whether the detector could
be calibrated accurately over a large range of concentrations
(22, 23). In the first studies, we found right away that the dye
CInject
V1, C1, Flow1 = C.O.
CRecirc
Lungs and Heart
V2, C2 Flow2
Kidney and Head
V3, C3 Flow3
G3
Liver and Gut
V4, C4 Flow4
Muscle
Fig. 8. Indocyanine Green Injection and Clearance. Upper panel: Blood concentrations in a
14 kg dog with 22 successive intravenous injections of 2.5 mg ICG (vertical pips). Dashed
lines represent estimated single exponential decay after each series of injections. Data
are from Edwards et al. (22). Lower panel: Circulatory model for hepatic clearance of ICG,
a mammillary compartmental model. Model code, with all parameters and equations, is in
Table 3. Modeling results are in Fig. 9.
was excreted via the bile: the feces turned green! Within a few years
the Indocyanine Green clearance Test was used as a liver function
test; it was a valued test even before the mechanisms of its hepatic
excretion became known (24).
The data from one dog (Edwards, 1960), shown in Fig. 8
(top panel), invite analysis. After each series of injections a set of
blood samples were obtained as the concentration diminished: the
diminutions appeared as straight lines on the semilog plot shown in
Fig. 8, top. This suggests a first-order clearance, i.e., a constant
fraction of the dye was being removed per unit time. Now let us
develop a simple model, and then test it against the data.
The anatomy and physiology provide the framework for the
model. The dye distributes throughout the whole blood volume.
We hypothesize that it is removed at one point only, by the liver,

and is excreted into the bile. We have no data on biliary concentra-
tion versus time, except that the dye ends up in the lower bowel
within the duration of the experiment, 4 h. We model the quanti-
tative data: the doses and their time of injection, the blood con-
centrations at the particular times, and use the dog’s weight, 14 kg,
as a constraint. The cardiac output is known from the areas of the
dye–dilution curves: C.O. ¼ Dose injected (mg)/Area under dye
curve (min. mg/l). These averaged 1.5 l/min, with a standard
deviation of about 10% (25).
A simple model of the whole body circulation is chosen as a
compromise, three compartments to represent fast, moderate, and
slow flows throughout the body and a single lumped compartment to
represent blood in the heart and lungs, as in Fig. 8, lower panel. A
mammillary model of this sort is far too crude to represent the
indicator dilution curves used to measure the cardiac output; this
requires higher temporal resolution (26, 27) and more precise mod-
els (28, 29). In this model there are up to three parameters for each of
the four compartments, flow and volume, and in compartment 3, a
consumption term G3. The shape of the transport function is fixed by
the assumption of a mixing chamber, so each has the impulse response
h(t) ¼ (1/t) exp(t/t), where the time constant t ¼ volume/flow.
The consumption influences the fraction leaving the compartment. If
we were to try to use the observed concentrations in Fig. 8 upper as
the input to the model we would see at once that we do not have
enough data: the points sampled are too sparse.
But the doses and their times of injection were defined pre-
cisely, so we use the sequence of doses as the input. We represent
the process of flow distribution and dilution of the injected ICG by
the set of four first-order differential equations, stirred tanks, as in
traditional compartmental analysis. Within each tank the concen-
tration, Cin, is given by
dCi =dt ¼ ðFi ðCini Ci Þ Gi Ci Þ=Vi ; (9)
where the index subscript i denotes the compartment, F is its flow, V
is its volume, and G is the rate of consumption within the volume.
Cini is the concentration in the blood entering the volume, and
equals the concentration in the compartment just upstream. The
concentration Ci is the same as the concentration flowing out
because of the basic compartmental assumption that the tank is
stirred instantaneously from entrance to exit. Via the test of fitting
the model to the data we can assess whether or not this assumption is
refuted by the data; if the data are fitted, then we do not prove that
the assumption is correct, but only that it was not invalidated by the
data. (This philosophical point underlies every assumption in for-
mulating a model: a model can never be proven correct. It can only
be demonstrated as adequate for the situation, and then can be used
as a “working hypothesis,” yet to be disproven. As T. H. Huxley put
it: “The great tragedy of science: the slaying of a beautiful hypothesis

by an ugly fact.” The downfall leads to the advancement!)
The JSim code is given in Table 3. The equations for each
compartment are closely similar. The recirculated dye and newly
injected dye enter the heart and lung compartment with central
blood volume V1, and mix with the total cardiac output. The only
clearance, G3, is during passage through the liver where ICG is
highly extracted and secreted into the bile; G2 and G4 are available
for exploration, but are set to zero.
The modeling results are shown in Fig. 9. Each injection
resulted in a sharp rise in concentration; when closely spaced the
mean concentration rises in spite of the rapid clearance. In the
absence of injections the concentration diminishes almost mono-
exponentially as was suggested by the straight lines on the semilog
plot in Fig. 8 (upper). If the system were a perfect single mixing
chamber with first-order clearance the time constant for washout
would be the volume divided by the clearance, Vtot/G3, and in fact
this is very close. The test is to plot the theoretical concentration
diminution against the model and the data. The equation,
Test ¼ 10 expððt 58Þ=ðVtot =G3 ÞÞ;
is the red dashed line in Fig. 9, and it does fit the phase where there
are no injections after 64 min and before the next injection at
100 min. The multiplier, 10, and the delay, 58 min, in the equation
simply position the curve. Test curves with longer delays (130 and
195 min for the purple and blue dashes) fit the data for the later
phases without injections. The closeness of these fits suggests that
approximating the whole blood volume as a single mixing chamber
would give a fairly good fit too. But the exponential Test curves
decay a little too rapidly compared to the model function and the
data at the lowest concentrations. This systematically better fit of
the model compared to the single exponential curve emphasizes
that the circulation is not really instantaneously wholly mixed, and
that a complete washout curve must be multi-exponential.
With respect to the physiological state of the animal, the fact
that the same value for G3 fits the data throughout the 4-h study
says that the hepatic clearance was stable over this long period of
anesthesia. We can also conclude that there is no evidence for the
saturation of the clearance process over the range of ICG concen-
trations seen here, fairly high concentrations, nearing 10 mg/l.
This implies that the processes for ICG clearance are not only
very effective but must also have binding constants high compared
to the concentrations found here. This dye soon found use as a
clinical test of liver function (24), and is widely used clinically today
(30–32). The hepatocyte’s apical transporter has a very high capac-
ity so that in the normal liver the transhepatic extraction is nearly
100% and the dye clearance can be used to estimate hepatic blood
flow. The reason the extraction is so high is that there is an active
ATP-supported extrusion of the dye from the hepatocyte into
Table 3
Compartmental model for hepatic ICG clearance
/* MODEL NUMBER: 0103

MODEL NAME: Comp4ICG
SHORT DESCRIPTION: Four-compartment whole body model with recirculation:
Repeated injections and first-order hepatic clearance of Indocyanine Green dye. */
import nsrunit; unit conversion on;
math Comp4ICG {
// INDEPENDENT VARIABLE
realDomain t min; t.min ¼ 0; t.max ¼ 30; t.delta ¼ 0.1;
// PARAMETERS
real Flow1 ¼ 1.5 L/min, // Blood Flow through Heart/Lung ¼ Cardiac Output
Dose ¼ 2.5 mg, // Amt injected
Vtot ¼ 1. 256 L, // Total Blood Volume
// First Compartmental unit (Heart/lung)
Vfr1 ¼ 0.25 dimensionless, // fraction of Vtot in Comp 1
V1 ¼ Vfr1*Vtot, // Volume of blood 1
C10 ¼ 0 mg/L, // Initial concentration
// Second Compartmental unit (Kidney and Head)
Fr2 ¼ 0.33, // Fraction of flow through comp 2
Flow2 ¼ Fr2*Flow1, // Flow through Comp 2
Vfr2 ¼ 0.40, // fraction of Vtot in Comp 2
G2 ¼ 0 ml/min, // consumption 2
// Third Compartmental unit (Liver and Gut)
Fr3 ¼ 0.25, // Fraction of Flow through Comp 3 (Liver)
Vfr3 ¼ 0.25, // fraction of Vtot in Comp 3
G3 ¼ 63.7 ml/min, // consumption 3
// Fourth Compartmental unit (Muscle)
Fr4 ¼ 1 - Fr2 - Fr3, // Fraction of Flow through Comp 4
Vfr4 ¼ 1 -Vfr1 -Vfr2 -Vfr3, // fraction of Vtot in Comp 4
G4 ¼ 0 ml/min; // Consumption 4
// Total flow relationship: Flow4 ¼ Flow1 - Flow2 - Flow3;
// EXTERNAL VARIABLE: The series of ICG injections, 2,5 mg each
extern real Cin1(t) 1/min; // 60 second Pulse injection @ 0.1 min
extern real Cin2(t) 1/min; // 60 second Pulse injection @ 30 min
(continued)
Table 3
(continued)

// DEPENDENT VARIABLES
real Crecirc(t) mg/L, //Inflow + recirculation
Qtot(t) mg, // Total amt in 4 compartments
Cinsum(t) 1/min, // Sum of the string of pulse injections
// Comp1 Comp2 Comp3 Comp4
C1(t) mg/L, C2(t) mg/L, C3(t) mg/L, C4(t) mg/L, // concn in each region
Q1(t) mg, Q2(t) mg, Q3(t) mg, Q4(t) mg; // Quantity in each regions
// INITIAL CONDITIONS
when(t ¼ t.min) {C1 ¼ C10; C2 ¼ C20; C3 ¼ C30; C4 ¼ C40;}
// INPUT CONCENTRATION FUNCTION
real Qinjrate(t) mg/min;
Cinsum ¼ (Cin1+ Cin2 + Cin3 + Cin4 + Cin5 + Cin6 + Cin7 + Cin8+
Cin9+ Cin10+ Cin11+ Cin12+ Cin13+ Cin14+ Cin15+
Cin16+ Cin17+ Cin18+ Cin19+ Cin20+ Cin21+ Cin22);
Qinjrate ¼ Dose*Cinsum;
Crecirc ¼ (C2*Flow2+ C3*Flow3 + C4*Flow4)/Flow1;
// QUANTITY Retained at time t:
Qtot ¼ V1*C1 + V2*C2 + V3*C3 + V4*C4 // Qtot is the total amount of ICG in the
body at time t.
real Area ¼ Dose*(integral(t ¼ t.min to t.max, Cinsum)); // check amt injected
// ORDINARY DIFFERENTIAL EQUATIONS
V1*C1:t ¼ Flow1*(Crecirc-C1) + Cinject; // Note that V1 can be left of ¼ sign
C2:t ¼ Flow2*(C1-C2)/V2 - G2*C2/V2;
C3:t ¼ Flow3*(C1-C3)/V3 - G3*C3/V3;
C4:t ¼ Flow4*(C1-C4)/V4 - G4*C4/V4;
} // END of program
/*
DETAILED DESCRIPTION:
This is a whole body model composed of a central blood volume from which flows the whole cardiac
output. The C.O. is distributed into three organs labeled "kidney" for kidney and head, "liver" for
liver and intestines, and "muscle" for the rest of the body.
(continued)
Table 3
(continued)
The four-compartment model has recirculation. The input function to compartment 1 (Heart and Lung)
is the sum of the recirculated indicator plus the series of injection pulses into V1, Qinjrate, each
injection at x mg/min for a short duration. Each injection, Cin#, is defined at run time using a separate
function generator. Clearance of the injected dye, indocyanine green, is hepatic extraction from the
blood via a saturable transporter on the hepatocyte sinusoidal membrane followed by ATP-dependent
excretion across the hapatocyte apical membrane into the bile, but is represented here by a passive first-
order loss, G3. This is adequate kinetically only at low concentrations of ICG, where the transporter is
mainly uncomplexed. */
KEY WORDS: compartment, flow and exchange, mixing chamber, hepatic clearance. first-order
consumption, washout, organ, multi-organ, recirculation
REFERENCES:
Edwards AWT, Bassingthwaighte JB, Sutterer WF, and Wood EH. Blood level of indocyanine green in
the dog during multiple dye curves and its effect on instrumental calibration. Proc S M Mayo Clin 35:
747-751, 1960.
Edwards AWT, Isaacson J, Sutterer WF, Bassingthwaighte JB, and Wood EH. Indocyanine green
densitometry in flowing blood compensated for background dye. J Appl Physiol 18: 1294-1304,
1963.
REVISION HISTORY:
Author: BEJ 06jan11
Revised by: JBB 09jan11 to combine the function generators, fgens, to speed computation
COPYRIGHT AND REQUEST FOR ACKNOWLEDGMENT OF USE:
Copyright (C) 1999-2011 University of Washington. From the National Simulation Resource,
Director J. B. Bassingthwaighte, Department of Bioengineering, University of Washington, Seattle WA
98195-5061.
Academic use is unrestricted. Software may be copied as long as this copyright notice is included.
This software was developed with support from NIH grant HL073598.
Please cite this grant in any publication for which this software is used and send an e-mail with the
citation and, if possible, a PDF file of the paper to [email protected].
*/ Table 3 is an example of a complete set of MML code in the format used in general for JSim project files, having the
same sequence and format as those used in the repository of models at www.physiome.org/Models.
the bile, thus keeping the intra-hepatocyte concentration low.

These features cannot be demonstrated at the low concentrations
seen in Figs. 8 and 9, but one would expect that the model would
have to be revised to include a saturable transporter if the concen-
trations were a lot higher.
5. Summary:
The Processes
Undertaken in
Pharmacokinetics In Subheading 4 we covered a standard approach to the steps in a
modeling analysis of data. The order of the steps depends a little on
the nature of the task. In the first model we performed no
Fig. 9. Model solution to fit the ICG data in a 14 kg dog. The triangles are the data shown in Fig. 8 (upper). The parameters
and initial conditions for the model are those given in the code in Table 3. The parameters are the result of optimization
using NL2SOL from Dennis and Schnabel (19); the total blood volume, Vtot, and the hepatic clearance, G3, were the only
free parameters. The 2.5 mg injection pulses are shown along the abscissa. The dashed lines are mono-exponentials with
time constant Vtot/G3. The model is #103 at http://www.physiome.org/jsim/models/webmodel/NSR/Comp4ICG/.
verification steps, and in the second the verification was done after
the fitting of the model to the data. This is clearly in the wrong
order: there is no point fitting a model to the data until it has been
demonstrated to be computed correctly, so in the list that follows,
the verification is done as soon as the draft model is constructed.
One cannot argue that the verification is not needed until the
model fits the data: a mathematically incorrect model might fit
the data, and after all that work the effort would be shown to be a
waste if the code had an error. Our failure to precede the data fitting
by a formal verification in these two cases is based on prior observa-
tions that JSim’s solutions for compartmental models provide four-
digit accuracy compared to analytical solutions.
Taking a listing of the steps to a detailed level:
l Ideally, design the experimental protocol to be the best test of
the model.
l Gathering supporting data, assess experimental accuracy.
l Obtain information on necessary parameters, a priori, and on
possible constraints.
l Complete development of the model. List all the assumptions.
l Conduct simulations, compare the results with other methods

of solution, perform verification tests, and find limiting cases
for which there are analytical solutions.
l Determine how one should weight the data, for use in mini-
mizing the sum of squares of differences between model solu-
tions and data.
l Validate that the model is reasonable, that it provides good fits
and that the, residual errors are not systematic or localized.
l Comparing simulation results to experiments and/or results of
other methods.
l Post-processing: Analyzing the results for consistency with
respect to physical and chemical and physiological expectations.
l Interpret the results scientifically. What does the model predict?
l Rethink the process. What are the weakest assumptions in the
model? How might it be improved? The model is always wrong.
Figure out new tests of it. Where might its predictions fail?
This overall process can be considered a success, if:
l Observations hitherto unexplained now fit a rational working
hypothesis.
l The essence of the phenomenology has been captured. (A
descriptive success only, perhaps.)
l Diagrams and schema of interactions aid understanding the
model, the kinetic relationships, and the code.
6. Model
Alternatives
and Modifications:
Interactive When the fit is not precise, outside of the limits of expectation
Hypothesis Revising relative to the noise in the data, despite all the attempts, then
maybe the model is just wrong. Certainly it is not nicely descriptive,
let alone explanatory! Given the philosophical premise that all
models are wrong, in the sense of being incomplete, or incorrect
mechanistically, every failure is a stimulus to find an alternative
model. A most rewarding and successful strategy is that of Platt
(12): he proposes that right from the outset one should
have alternative hypotheses in mind, and that the experiment
should be designed to distinguish between these hypotheses.
“Strong Inference” is the title of his paper. We advocate that each
hypothesis be expressed in terms of a computational model, since
that means that it is described explicitly and is therefore testable.
The strategy pays off because at least one of the hypotheses
is proven wrong when the model fails to fit the data. Sometimes
both are wrong! Regroup, rethink!
The least squares approach, minimization of the sums of

squares of the differences between the data points and the model
function, is a blunt tool, without specific information with respect
to any misfitting. Displaying the residuals, plotting the point-by-
point differences as a function of time (or in general, of the
independent variable), is an excellent way to get insight into
what to do next: a series of points above or below zero reveals a
systematic misfitting. Comparing the time course of the deviations
of the data from the model with the sensitivity functions of the
various parameters could suggest that a particular parameter has
not been optimized well, but this is most unlikely when the
various options for optimization have been explored. It is more
likely that the model lacks a feature that is needed. Back to the
drawing board!
This is the usual iterative process: hypothesis in general terms,
develop the model that represents a clear precise hypothesis, if
possible design the model for an alternative hypothesis, design
the distinguishing experiment while taking into account the accu-
racy of the data to be acquired, and disprove one or both hypoth-
eses. The next step is to improve the model and repeat the series of
steps until a satisfactory level of synchrony between model and
experiment is achieved. This version of the model is usually then
designated the working hypothesis. No working hypothesis is to be
regarded as “the truth,” though it is useful for practical purposes.
And in fact it serves as the standing target to be disproved in order
to advance the science.
7. What to Do When
the Compartmental
Representation Is
not so Good? What follows is a common example that applies in biophysics,
physiology, and pharmacology. It is usual that drugs and substrates
for metabolism and signaling molecules of molecular weight less
than 1,000 Da are partially extracted during single passage through
a capillary. Since capillaries are about a 1,000 mm long, but only
5 mm in diameter (as in Fig. 10), there is no possibility that they are
instantaneously stirred tanks with uniform concentration from end
to end: there must be gradients between capillary entrance and exit
for any solute that is exchanging between blood and tissue. If the
extraction is less than 5%, the gradient might be ignored, but for
solutes of interest to us here the steady-state extractions are
30–90%, and so affect the estimates of the permeabilities and
consumption rates. Consequently we now consider the computa-
tional differences between a stirred tank and an axially distributed
capillary–tissues exchange unit.
Fig. 10. A venule and capillaries on the epicardium of a dog heart casted with microfil.
Capillary diameters are 5.6 1.3 mm, average intercapillary distances are 17–19 mm,
and lengths are 800–1,000 mm. The distance between the long calibration lines is
100 mm. Modified from Bassingthwaighte et al. (33).
7.1. Capillary–Tissue A two-compartment system is here modified to incorporate flow, as

Exchange: Convention, shown in Fig. 11 (left panel), thus identifying V1 as the vascular
Permeation, Reaction, region, the membrane as the capillary barrier, and V2 as the tissue.
and Diffusion The system is considered as a homogeneously perfused organ with
constant volumes and steady flow, F, in and out. Now, in order to
put it into the context of substrate delivery and metabolism, we
switch to standard physiological representation of the units, defin-
ing them per gram of organ mass. F, PS, and the consumption G
have units ml/(g min), and the volumes have units ml/g. This
notation normalizes flows, substrate use, etc. to be independent
of organ mass. (The model is www.physiome.org/jsim/models/
webmodel/NSR/Comp2FlowExchange/, #247.)
To keep the system simple so as to focus on the blood–tissue
exchange, the intratissue consumption is considered to be a first-
order process, as if the substrate concentration is far below the KD
for any enzymatic reaction.
In the right panel of Fig. 11 is the equivalent axially distributed
model that accounts for gradients in concentration along the length
of the capillary. It is more general, but reduces to an exactly analo-
gous model when the PS’s are set to zero, so there is no entry into
the cells. Other parameters for cellular permeation and reaction are
as follows: for passage across endothelial cell luminal membrane
(PSecl); endothelial cell abluminal membrane (PSeca); and paren-
chymal cell membrane (PSpc); G, intracellular consumption ml/
(g min) or metabolism of solute by endothelial cells (Gec) or by
parenchymal cells (Gpc). The Vs, ml/g, are volumes of distribution
in plasma (Vpl), endothelial cell (Vec), interstitial (Visf), and paren-
chymal cell (Vpc). With the cell permeabilities, PSecl, PSeca, and
Fig. 11. Compartmental versus axially distributed models for capillary–tissue exchange. Exchange between capillary
plasma and interstitial fluid regions can be regarded as providing two conceptually similar but mathematical distinguish-
able methods of representation. Left panel: Two compartment stirred tank for the exchange of solute C between flowing
blood and surrounding stagnant tissue. Flow F, ml/(g min), carries in solute at concentration Cin, mM, and carries out a
concentration Coutt ¼ Cp, mM. The flux from capillary to ISF (compartment 2) is limited by the conductance PS ml/(g min),
the permeability-surface area product of the membrane separating the two chambers, allowing bidirectional flux. G2 ml/
(g min) is a reaction rate for a transformation flux forming product at a rate G2·C2 mmol/(g min). Right panel: Axially
distributed model equivalent to the two-compartmental model when the solute does not enter the endothelial or
parenchymal cells. PSg, the conductance for permeation through the interendothelial clefts, is equivalent to the PS of
the compartmental model in the left panel. Right panel from Gorman et al. (18) with permission from the American
Physiological Society.
PSpc, set to zero, PSg is equivalent to the compartmental PS, and

Vpl and Visf are equivalent to V1 and V2. In each of the four regions
there is an axial dispersion coefficient, equivalent to a diffusion
coefficient, setting the rate of random spreading along the length.
The full version of this model, GENTEX, a general multispecies
model for blood-tissue exchange is available at: www.physiome.
org/jsim/models/webmodel/NSR/gentex.proj.
Chinard (34) developed a technique to distinguish individual
processes involved in blood–tissue exchange and reaction: the
Multiple-Indicator Dilution (MID) technique (Fig. 12). He first
used it for the purpose of estimating the volumes of distribution for
sets of tracers of differing characteristics: the mean transit time vol-
ume, Vmtt ¼ F t where t is the mean transit time through the
system. He did not estimate permeabilities as his studies were on
highly permeable solutes. Crone (35) analyzed the technique,
showing how it could be used to estimate PS from the outflow curves
for a simultaneous injection into the inflow of a solute and imperme-
able reference intravascular tracer as shown in Fig. 12. The figure dia-
grams an experimental setup for examining the uptake of D-glucose in
an isolated perfused heart as by Kuikka et al. (36). L-Glucose, the
stereoisomer, serves as an extracellular, non-metabolized reference.
A more realistic diagram of a capillary–tissue exchange includes the
endothelial cells and interstitial fluid (ISF), as shown in Fig. 11.
To determine capillary permeability, the relevant reference sol-
ute is one that does not escape from the capillary blood during
single transcapillary passage; for example, albumin is the relevant
reference solute to determine the capillary permeability to glucose.
In this situation the albumin dispersion along the vascular space
Fig. 12. Schematic overview of experimental procedures underlying the application of the multiple-indicator dilution
technique to the investigation of multiple substrates passing through an isolated organ without recirculation of tracer. The
approach naturally extends also to their metabolites.
may be assumed to be the same as that of the glucose; thus

the shape of the albumin impulse response, halb(t), accounts for
the intravascular transport of all the solutes. (L-Glucose, an extra-
cellular reference tracer with the same molecular weight and diffu-
sivity as D-glucose, is the extracellular reference for D-glucose,
having the same capillary PSg and the same interstitial volume of
distribution, Visf. Having simultaneous data on such reference tra-
cers greatly reduces the degrees of freedom in estimating the para-
meters of interest for D-glucose.)
7.2. Model Equations The two diagrams in Fig. 11 look quite different, but the second
for Tracer can be reduced to the compartmental model, as we will show
below. The essential difference is that the distributed model
accounts for concentration gradients along the capillary length.
Capillaries are about 1 mm long, and are 5 mm in diameter, an
aspect ratio of 200. Diffusional relaxation times thus differ by a
factor of 200 between radial and axial directions. Consequently,
considering the capillary as a stirred tank is unreasonable.
The stirred tank expressions account for the flow through
compartment 1, the permeation, and consumption terms G2 ml/
(g min) in the second compartment:
dC1 PS F
¼ ðC1 C2 Þ ðCin C1 Þ;
dt V1 V1
dC2 PS G2
¼ þ ðC1 C2 Þ C2 : (10)
dt V2 V2
The use of these ODEs implies and builds into the calculations
a discontinuity between the concentration of solute in the inflow
and that in V1. Because V1 is assumed instantly mixed, there is no

gradient along the capillary and the tracer entering the tank is
immediately available to be washed out with the same probability
as any molecule dwelling in there for a longer time.
Alternatively the capillary–tissue unit of Fig. 11 right, can be
reduced to two regions represented by PDEs that allow a continu-
ous gradient along the length of the capillary between entrance and
exit, like those of Krogh and Erlang (37), Sangren and Sheppard
(38), and Renkin (39). It improved upon these by incorporating
axial dispersion in both the capillary and extravascular regions.
Using the spatially distributed analogs for plasma, Cpl, or blood,
and extravascular tissue, Cisf, instead of the lumped variables C1 and
C2:
@Cp ðx; tÞ Fpl L @Cpl PSg @ 2 C1
¼ Cpl Cisf þ Dx1 ;
@t Vpl @x Vpl @x 2
@Cisf ðx; tÞ PSg Gisf @ 2 C2
¼ Cpl Cisf Cisf þ Dx2 ; (11)
@t Visf Visf @x 2
where Cpl and Cisf are spatially distributed functions of both x and t,
not just t. The axial position is denoted by x, where 0 < x < L, the
capillary length, cm. The analogy between this model, Eq. 11, and
the compartmental version in Eq. 10 is Fp ¼ F, Vp ¼ V1, Visf ¼ V2,
and PSc ¼ PS, the permeability-surface area of the capillary wall,
but we retain the two sets of names in order to allow comparisons
between the estimated parameter values. The capillary length, L, is
arbitrarily set to an average value such as 0.1 cm; in the computa-
tions what is being used is the dimensionless fractional length, x/L.
PDEs require boundary conditions. At the capillary entrance, in
contrast to the compartmental model there is no discontinuity in the
concentration profile, but there is a requirement for matching the
diffusional and convective terms so that the influx is just the con-
vected mass, FCp(x ¼ 0, t). The boundary conditions are written:
when ðx ¼ xmin ÞfðFp L=Vp Þ ðCp Cin Þ þ Dp dCp =dx ¼ 0; g
at inlet to capillary,
(12)

when ðx ¼ xmax Þ dCp =dx ¼ 0; Cout ¼ Cp ;
at exit from capillary: (13)
The form of the inlet condition is important when the diffusion is
large, and we use it here for conceptual and practical accuracy. The
outflow concentration Cout is set equal to the concentration just inside
the exit, Cp(x ¼ L, t), the same condition described by the ODEs.
The last term in each equation is the diffusion along the length
of the capillary–tissue regions; the use of an anatomically correct
length then makes using observed diffusion coefficients for Dp and
Dtiss, cm2/s, practical and meaningful. Gross exaggeration of the
diffusion coefficients can be used in the equations to turn the

distributed model into a de facto well-mixed, compartmental
model. With both diffusion coefficients set to infinity this model
behaves identically to the compartmental model of Eq. 10.
The negative sign on the flow term, the first term on the right
of the equal sign in Eq. 11, merits further explanation. Consider
the inflow to contain a bolus of solute: as it enters, the concentra-
tion at the capillary entrance rises. At this time, the slope of the
curve of concentration versus position x, ∂C/∂x, is negative as
illustrated in the lower panel of Fig. 13 by the slope of C(x, t) for
the bolus shape at t ¼ 1.5 s, at the capillary mid point, x ¼ 0.5 L.
The spatial slope has always the sign opposite to the temporal
derivative ∂C/∂t at the same point, thus the negative sign on the
term.
Functionally, therefore Eq. 11 is analogous to Eq. 10. But
using the PDEs avoids the unrealistic discontinuity in the compart-
mental model at the entrance. Obnoxious, unrealistic discontinu-
ities at the entrances are the consequence of the instantaneous mixing
within a compartment. Using the PDE allows continuity in concen-
trations and concentration gradients along the capillary, and not
only in concentration but also in the properties of the system such
as axial gradients in transporter and enzyme densities that are
evident in the liver sinusoid. For the following analysis, all para-
meters are assumed spatially uniform so as to minimize the differ-
ence from the compartmental models. There are many ways of
representing axially distributed convecting systems, and two are
shown in Fig. 13. One solves PDEs using a PDE solver (there are
several choices within JSim), and here we used a Lagrangian
method (40–42). A compartmental type of alternative is to approx-
imate the capillary as a series of stirred tanks, each with the same
volume and PS, as diagrammed in the upper part of the figure. With
a large number of serial stirred tanks, the longitudinal concentra-
tion gradient is approximated increasingly accurately as the steps
from one to the next are small. The intravascular transport process
with serial stirred tanks is a Poisson process. In modeling, serial
stirred tanks are convenient because the number of tanks, Ntanks,
can be used as a free parameter. The relative dispersion RD over the
length of the tube is determined by Ntanks such that the RD of the
outflow curve (RD equals the coefficient of variation, the standard
deviation divided by the mean transit time), induced during transit,
is the reciprocal of the square root of Ntanks, so that with 100 tanks
the RD is 10%.
Figure 13 (upper) shows that curves for the PDE solution and
for the Poisson process are essentially similar, so that the dispersion
coefficient Dp sufficed to create the same dispersion as occurred
with the Poisson process using 109 tanks. The choice of 109 tanks
is arbitrary, large that a plot of C(Ntanks) versus N would appear
smooth. Because the capillary PS > 0, there is loss of solute as the
Fig. 13. Pulse responses in axially distributed models. The input function, Cin, is a pulse of duration 1.4 s. Upper panel:
Outflow concentration–time curves for a partial differential equation solution using a Lagrangian sliding fluid element
method and an intravascular dispersion coefficient, Dp ¼ 2.6 105 cm2/s (gray curve), and for a serial stirred tank
algorithm representing a Poisson process with 109 stirred tanks (black curve almost superimposed on the gray one). Lower
panel: Intracapillary spatial profiles in the distributed model (using the PDEs) at a succession of times, 1.5, 2.0, 2.5, and
3.0 s. The pulse slides and disperses due to the diffusion while some solute is lost from the vascular space by permeation
of the capillary wall. Parameters were the same for the compartmental 109 tank Poisson model and the PDE: Fp ¼ 1 ml/
(g min), PSC ¼ 2 ml/(g min), and tissue volume Vtiss was set to 10 ml/g so that there was negligible tracer flux from tissue
back into the plasma space. (The model: http://www.physiome.org/jsim/models/webmodel/NSR/Anderson_JC_2007/FIG-
URES/Anderson_JC_2007_fig11/index.html) Figure from Anderson (6).
bolus progresses along the capillary. The permeative loss is the same
for both methods, with the result that the peak outflow concentra-
tions are similar. Figure 13 (lower) shows the shape of the bolus as a
function of position as it deforms continuously from its initial
square pulse at the entrance to the capillary. The diminution in
peak height is therefore due not only to the spreading but also to
the loss. This loss is reflected of course in the reduction in the areas.
Fig. 14. Effect of reduction of Ntanks on the output C(t ). Responses of the Nth order
Poisson operator with Ntanks varied from 109 tanks in series down to 50, 20, 10, 5, 2,
and finally to a single mixing chamber, Ntanks ¼ 1. The gray curve is the Lagrangian
solution to the PDEs as in Fig. 13. All of the Poisson operator outflow curves (black) have
the same mean transit time, and the same parameters: Vp ¼ 0.05 and Vtiss ¼ 0.15 ml/g;
Fp ¼ 1 ml/(g min), PSg ¼ 1 ml/(g min). Figure from Anderson (6). Model is #46, running
loop mode to change Ntanks (http://www.physiome.org/jsim/models/webmodel/NSR/
Anderson_JC_2007/FIGURES/Anderson_JC_2007_fig12/index.html).
The grey curve touching to the top of each spatial profile is the
theoretical curve from Crone (35) and Renkin (39), as in the model
of Sangren and Sheppard (38):
PSg x

CðxÞ ¼ 1 e Fp L ; (14)
where x/L is the fractional distance along the capillary, the abscissa
in the lower panel.
Now that we know that the multi-compartmental serial tank
representation can give results approximating the normal PDE
representation, and that they differ basically only in the numerical
method used, the question becomes: “Which of the methods pro-
duces the correct assessment of the parameters with the greatest
efficiency?” The serial stirred tank model has the disadvantage that
the waveforms are seriously distorted by reducing the number of
stirred tanks, as is shown in Fig. 14.
While reducing the number of tanks in the stirred tank method
has a dramatic effect on the shapes of the outflow curves, the
problem is much less severe with the PDE representation, as
shown in Fig. 15. Solutions are shown for Ngrid ¼ 109, 51, 21,
11, and 7 for two methods of solving PDEs, one using a robust
solver TOMS731 (43) and the other using a Lagrangian sliding
fluid element algorithm (42).
Fig. 15. Effect of reduction of Ngrid on the output C(t ) with two PDE solvers. The green curve is the serial compartmental
model with 109 tanks. The red rectangle is the input function divided by 4. Black curves are the solutions to the PDE
at varied resolutions. Upper panel: Lagrangian sliding fluid element algorithm, LSFEA, with Ngrid ¼ 121 black solid line,
51 short dashes superimposed on the black solid line, 21 long dashes, and 11 dotted. This algorithm, while computation-
ally fast, describes the outflow dilution curve as a series of square pulses; with a large number of segments the curves
appear smooth. With fewer segments, e.g., Ngrid ¼ 21 (long dashes) or 11 (dotted line), the steps are obvious but the
approximation is reasonably good. Lower panel: TOMS731 solver. This slow, robust solver, like most PDE solvers,
broadens the solution somewhat with reduced Ngrid but with Ngrid ¼ 11 (dotted line) there is obvious spreading and
oscillation in the solutions (www.physiome.org Model # 126).
Fig. 16. Outflow dilution curves for D-glucose, 2-deoxy-D-glucose, and albumin (dog expt
4048-6) fitted with the model. The deoxyglucose curve has been shifted downward by half
a logarithmic decade (ordinate values divided by 10 in.) to display it separately from D-
glucose curve. Parameter estimates for D- and deoxyglucose were PSc, 0.97 and 1.0;
PSpc ¼ 0.7 and 0.5; Gpc ¼ 0.01 and 0.05 ml/(g min); the volume ratio Visf/Vp ¼ 6.5 and
Vpc/Vp ¼ 13.3 ml/g, the same for both glucoses. Coefficients of variation were 0.19 and
0.09. From Kuikka (36) with permission from the American Physiological Society (Model is
at www.physiome.org: Search on Kuikka, model 126).
Computation times for the PDE solvers are almost propor-

tional to Ngrid, the number of spatial elements chosen for the
computation. LSFEA is very much faster than TOMS731, and is
also faster than the serial tanks model since fewer segments are
needed. The key thing is that when the number of grid points is
reduced the PDE algorithms do not change shape so much as the
serial compartmental models. The result is that the PDEs are more
efficient for fitting the data and provide a realistic solution.
Accounting correctly for smooth intracapillary gradients is
important in the analysis of indicator dilution data for low-molecu-
lar-weight solutes using high-resolution techniques. Figure 16 is an
example of modeling of the uptake of glucoses in the dog heart.
The data from Kuikka (36) are outflow dilution curves from actively
contracting blood-perfused dog hearts. Albumin is the intravascu-
lar reference data curve defining the dispersion and delay from the
coronary inflow to the effluent veins for all of the tracers. The study
was on D-glucose, a normal substrate, and 2-deoxy-D-glucose, an
abnormal glucose that is transported into the cardiomyocyte, phos-
phorylated by hexokinase, and then cannot be further metabolized.
The model is that illustrated in Fig. 11 (right panel) but without
considering the negligible uptake by the endothelial cells. Para-
meters for the capillary wall permeation through the clefts and for
the cellular uptake are given in the legend. The quality of the data is
shown by the smoothness of the curves over time; the goodness of
fit is a testimonial to the quality of the model defined through
the PDEs.
8. Commentary
This presentation has emphasized issues important for new users of

compartmental systems to contemplate. This reduced coverage
neglects many important issues such as the identifiability of the
models, distinguishing one option from another, the numbers
and accuracy of the data points and their timing in obtaining
parameter estimates, the consequences of model verification test-
ing, and the basic principles of numerical methods and optimiza-
tion techniques. These are all covered in the more general
references listed in the Introduction.
We did however carefully express the models in physiological
terms, trying to emphasize recognition of the nature of the pro-
cesses of exchange and reaction. Thus Eq. 10 was deliberately not
written in the common parlance used in compartmental modeling:
dC1
¼ k1 ðC1 C2 Þ k2 ðCin C1 Þ;
dt
dC2
¼ þk3 ðC1 C2 Þ k4 C2 ; (15)
dt
where k1 ¼ PS/V1, k2 ¼ F/V1, k3 ¼ PS/V2, and k4 ¼ G2/V2.
These equations look simpler, for it appears that there are fewer
parameters, namely, four, while the original equations appear to
have five: F, PS, V1 and V2, and G2. But they are deceptively simple,
quietly masking their true identities behind the fact that there are
really only four independent parameters in the original Eq. 10.
Putting t ¼ V1/F, g ¼ V2/V1, d ¼ G2/F, and e ¼ PS/F, we
have k1 ¼ PS/V1 ¼ d/t; k2 ¼ F/V1 ¼ 1/t; k3 ¼ PS/V2 ¼ e/
(gt), and k4 ¼ G2/V2 ¼ d/(gt); this illustrates that there were
only four parameters in the original expressions, as expected,
remembering that the exponential response of a single compart-
mental system has just one parameter t ¼ V1/F, its mean transit
time. The point is that PS, G, F, and Vs have real anatomic and
physiological meaning, and reminding us that when the flow is
known, V1 is the relevant unknown, and it is usually constrained
by knowledge of the anatomy.
In general, “data” should include more than just a set of con-
centration–time curves, for one must recognize the underlying
anatomy as data. Anatomic data constrain the estimates of water
space or physical volumes as shown by Yipintsoi (44) and Vinnakota
(45). Further, large comprehensive models can be used to great

advantage when built upon anatomic and physicochemical con-
straints combined into the physiological kinetics (46).
In this essay we have not emphasized the distinction between
tracer and non-tracer kinetics except in the Introduction in
Eqs. 1–3. We have however considered the nonlinearity of reactions
in Application 4A, Aspirin Kinetics, where the dependence of the
reaction flux on the concentration is evident. If the salicylate con-
centration were held constant at any particular value, then for a
tracer, the tracer flux would be determined entirely by the concen-
tration of the ambient non-tracer salicylate, and the system would
be reduced to a linear one with constant coefficients (47). Likewise,
when examining tracer fluxes in a situation where mother substance
is varying, this should be taken into account by computing the dual
model, for mother and tracer simultaneously, a consideration all
too often forgotten in the common usage of compartmental analy-
sis. The warning is: when using tracers, measure the background
concentrations of mother substance often enough to assure its
constancy.
The traditional compartmental analysis therefore has its great-
est strength in situations where the overall system is in steady state
for all substances related to the tracer substance, for then the power
of linear systems analysis applies, convolution integration and sta-
tionarity apply, and matrix inversion can be used on sets of linear
ODEs. Then all is well set up for linear systems analysis.
References
1. Berman M (1963) The formulation and testing 7. Knopp TJ, Anderson DU, Bassingthwaighte JB
of models. Ann N Y Acad Sci 108:182–194 (1970) SIMCON–Simulation control to opti-
2. Jacquez JA (1972) Compartmental analysis in mize man-machine interaction. Simulation
biology and medicine. Kinetics of distribution 14:81–86
of tracer-labeled materials. Elsevier Publishing 8. Sauro HM, Fell DA (1991) SCAMP: a meta-
Co, Amsterdam, 237 pp bolic simulator and control analysis program.
3. Jacquez JA (1996) Compartmental analysis in Math Comput Model 15:15–28
biology and medicine, 3rd edn. BioMedware, 9. Sauro HM, Hucka M, Finney A, Bolouri H
Ann Arbor, MI, 514 pp (2001) The systems biology workbench con-
4. Cobelli C, Foster D, Toffolo G (2000) Tracer cept demonstrator: design and implementa-
kinetics in biomedical research. From data to tion. Available via the World Wide Web at
model. Kluwer Academic, New York http://www.cds.caltech.edu/erato/sbw/
5. Zierler KL (1981) A critique of compartmental docs/detailed-design/
analysis. Annu Rev Biophys Bioeng 10. Raymond GM, Butterworth E, Bas-
10:531–562 singthwaighte JB (2003) JSIM: free software
6. Anderson JC, Bassingthwaighte JB (2007) package for teaching physiological modeling
Tracers in physiological systems modeling. and research. Exp Biol 280.5:102. (www.phy-
Chapter 8 Mathematical modeling in nutrition siome.org/jsim)
and agriculture. In: Mark D. Hanigan JN, 11. Chizeck HJ, Butterworth E, Bassingthwaighte
Casey L Marsteller. Proceedings of the ninth JB (2009) Error detection and unit conversion.
international conference on mathematical Automated unit balancing in modeling inter-
modeling in nutrition, Roanoke, VA, 14–17 face systems. IEEE Eng Med Biol 28(3):50–58
August 2006, Virginia Polytechnic Institute 12. Platt JR (1964) Strong inference. Science
and State University Blacksburg, VA, 146:347–353
pp 125–159
13. Bassingthwaighte JB, Chinard FP, Crone C, 26. Stewart GN (1897) Researches on the circula-
Goresky CA, Lassen NA, Reneman RS, Zierler tion time and on the influences which affect it:
KL (1986) Terminology for mass transport and IV. The output of the heart. J Physiol
exchange. Am J Physiol Heart Circ Physiol 22:159–183
250:H539–H545 27. Hamilton WF, Moore JW, Kinsman JM, Spur-
14. Benedek IH, Joshi AS, Pieniazek JH, King ling RG (1932) Studies on the circulation. IV.
S-YP, Kornhauser DM (1995) Variability in Further analysis of the injection method, and of
the pharmacokinetics and pharmacodynamics changes in hemodynamics under physiological
of low dose aspirin in healthy male volunteers. and pathological conditions. Am J Physiol
J Clin Pharmacol 35:1181–1186 99:534–551
15. Aarons L, Hopkins K, Rowland M, Brossel S, 28. Thompson HK, Starmer CF, Whalen RE,
Thiercelin JF (1989) Route of administration McIntosh HD (1964) Indicator transit time
and sex differences in the pharmacokinetics of considered as a gamma variate. Circ Res
aspirin, administered as its lysine salt. Pharm 14:502–515
Res 6:660–666 29. Bassingthwaighte JB, Ackerman FH, Wood
16. Prescott LF, Balali-Mood M, Critchley JAJH, EH (1966) Applications of the lagged normal
Johnstone AF, Proudfoot AT (1982) Diuresis density curve as a model for arterial dilution
or urinary alkalinisation for salicylate poison- curves. Circ Res 18:398–415
ing? Br Med J 285:1383–1386 30. Krenn CG, Krafft P, Schaefer B, Pokorny H,
17. Chan IS, Goldstein AA, Bassingthwaighte JB Schneider B, Pinsky MR, Steltzer H (2000)
(1993) SENSOP: a derivative-free solver for Effects of positive end-expiratory pressure on
non-linear least squares with sensitivity scaling. hemodynamics and indocyanine green kinetics
Ann Biomed Eng 21:621–631 in patients after orthotopic liver transplanta-
18. Glad T, Goldstein A (1977) Optimization of tion. Crit Care Med 28:1760–1765
functions whose values are subject to small 31. Krenn CG, Pokorny H, Hoerauf K, Stark J,
errors. BIT 17:160–169 Roth E, Steltzer H, Druml W (2008) Non-
19. Dennis JE, Schnabel RB (1983) Numerical isotopic tyrosine kinetics using an alanyl-
methods for unconstrained optimization and tyrosine dipeptide to assess graft function in
nonlinear equation. Prentice-Hall, New York liver transplant recipients - a pilot study. Wien
20. Fox IJ, Brooker LGS, Heseltine DW, Essex Klin Wochenschr 120(1–2):19–24
HE, Wood EH (1957) A tricarbocyanine dye 32. Kortgen A, Paxian M, Werth M, Recknagel P,
for continuous recording of dilution curves in Rauschfusz F, Lupp A, Krenn C, Muller D,
whole blood independent of variations in Claus RA, Reinhart K, Settmacher U, Bauer
blood oxygen saturation. Proc Staff Meet M (2009) Prospective assessment of hepatic
Mayo Clin 32:478 function and mechanisms of dysfunction in
21. Edwards AWT, Isaacson J, Sutterer WF, Bas- the critically ill. Shock 32(4):358–365
singthwaighte JB, Wood EH (1963) Indocya- 33. Bassingthwaighte JB, Yipintsoi T, Harvey RB
nine green densitometry in flowing blood (1974) Microvasculature of the dog left ventric-
compensated for background dye. J Appl ular myocardium. Microvasc Res 7:229–249
Physiol 18:1294–1304 34. Chinard FP, Vosburgh GJ, Enns T (1955)
22. Edwards AWT, Bassingthwaighte JB, Sutterer Transcapillary exchange of water and of other
WF, Wood EH (1960) Blood level of indocya- substances in certain organs of the dog. Am J
nine green in the dog during multiple dye Physiol 183:221–234
curves and its effect on instrumental calibra- 35. Crone C (1963) The permeability of capillaries
tion. Proc Staff Meet Mayo Clin 35:747–751 in various organs as determined by the use of
23. Bassingthwaighte JB (1966) Plasma indicator the indicator diffusion method. Acta Physiol
dispersion in arteries of the human leg. Circ Scand 58:292–305
Res 19:332–346 36. Kuikka J, Levin M, Bassingthwaighte JB
24. Hunton DB, Bollman JL, Hoffman HN (1986) Multiple tracer dilution estimates of D-
(1961) The plasma removal of indocyanine and 2-deoxy-D-glucose uptake by the heart.
green and sulfobromophthalein: effect of dos- Am J Physiol Heart Circ Physiol 250:
age and blocking agents. J Clin Invest 30 H29–H42
(9):1648–1655 (PMCID PMC290858) 37. Krogh A (1919) The number and distribution
25. Bassingthwaighte JB, Edwards AWT, Wood of capillaries in muscles with calculations of the
EH (1962) Areas of dye-dilution curves sam- oxygen pressure head necessary for supplying
pled simultaneously from central and periph- the tissue. J Physiol (Lond) 52:409–415
eral sites. J Appl Physiol 17:91–98
38. Sangren WC, Sheppard CW (1953) A mathe- 44. Yipintsoi T, Scanlon PD, Bassingthwaighte JB
matical derivation of the exchange of a labeled (1972) Density and water content of dog
substance between a liquid flowing in a vessel ventricular myocardium. Proc Soc Exp Biol
and an external compartment. Bull Math Bio- Med 141:1032–1035
phys 15:387–394 45. Vinnakota K, Bassingthwaighte JB (2004)
39. Renkin EM (1959) Transport of potassium-42 Myocardial density and composition: a basis
from blood to tissue in isolated mammalian for calculating intracellular metabolite concen-
skeletal muscles. Am J Physiol 197:1205–1210 trations. Am J Physiol Heart Circ Physiol 286:
40. Bassingthwaighte JB (1974) A concurrent flow H1742–H1749
model for extraction during transcapillary 46. Bassingthwaighte JB, Raymond GR, Ploger
passage. Circ Res 35:483–503 JD, Schwartz LM, Bukowski TR (2006)
41. Bassingthwaighte JB, Wang CY, Chan IS GENTEX, a general multiscale model for
(1989) Blood-tissue exchange via transport in vivo tissue exchanges and intraorgan
and transformation by endothelial cells. Circ metabolism. Phil Trans Roy Soc A Math Phys
Res 65:997–1020 Eng Sci 364(1843):1423–1442. doi:
42. Bassingthwaighte JB, Chan IS, Wang CY 10.1098/rsta.2006.1779
(1992) Computationally efficient algorithms 47. Bassingthwaighte JB, Goresky CA, Linehan JH
for capillary convection-permeation-diffusion (1998) Ch. 1 Modeling in the analysis of the
models for blood-tissue exchange. Ann processes of uptake and metabolism in the
Biomed Eng 20:687–725 whole organ. In: Bassingthwaighte JB, Goresky
43. TOMS./TOMS. Association of computing CA, Linehan JH (eds) Whole organ approaches
machinery: transactions on mathematical soft- to cellular metabolism. Springer, New York,
ware. http://www.netlib.org/toms/index. pp 3–27
html
View publication stats

Compartmental Modeling in The Analysis of Biological Systems

Uploaded by

Copyright:

Available Formats

Compartmental Modeling in The Analysis of Biological Systems

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Compartmental Modeling in The Analysis of Biological Systems

Uploaded by

Copyright:

Available Formats

See discussions, stats, and author profiles for this publication at: https://www.researchgate.

Compartmental modeling in the analysis of biological systems

Article in Methods in molecular biology (Clifton, N.J.) · January 2012

James Bassingthwaighte Bartholomew E Jardine

SEE PROFILE SEE PROFILE

Hepatic Clearance of Indocyanine green View project

The user has requested enhancement of the downloaded file.

Compartmental Modeling in the Analysis

Key words: Physiological and pharmacologic modeling, PKPD, Pharmacokinetics–pharmacody-

both administration and distribution, a relatively new term, PBPK,

distribution sites of positron-labeled drug may be shown with high

1.2.2. In Terms 1. Catenary system: Two or more compartments arranged in

mixing throughout the compartment, and therefore uniform inter-

1.3. Distribution Distribution, the D of ADME, is by convection, diffusion, and

necessarily well characterized for the particular drug. Heterogeneities

Cellular Open Resource (COR) is for a Windows environment:

Thus JSim, since it is designed around the analysis of experimen-

cannot be reproduced is misleading and wastes time and money.

when importing models from CellML or SBML which sometimes

experimental studies in which the data are acquired. The volumes of

by clicking “Debug” (left, bottom) to drop down a menu including

multiplication by 2 on each run, or the automatic addition of 3 on

3.2.4. Sensitivity Analysis Sensitivity analysis is available to provide a quantitative measure of

Sensitivities are local linear approximations, and are calculated

3.2.5. Optimization in Data Optimization, either by manual parameter adjustment or by

clearance was enzymatically mediated: at high concentration the dimi-

/* MODEL NUMBER 280

thermodynamically reversible, at least in principle. This reaction

approximately linear, as expected, whereas at the low concentra-

higher or due to “product pressure” to form salicylate by the

appropriate; it is almost never appropriate if the data range over one

concentration independent of the blood oxygen level. The dye,

We hypothesize that it is removed at one point only, by the liver,

it: “The great tragedy of science: the slaying of a beautiful hypothesis

/* MODEL NUMBER: 0103

extern real Cin9(t) 1/min; // 60 second Pulse injection @ 99 min

the bile, thus keeping the intra-hepatocyte concentration low.

l Conduct simulations, compare the results with other methods

The least squares approach, minimization of the sums of

7.1. Capillary–Tissue A two-compartment system is here modified to incorporate flow, as

PSpc, set to zero, PSg is equivalent to the compartmental PS, and

may be assumed to be the same as that of the glucose; thus

and that in V1. Because V1 is assumed instantly mixed, there is no

diffusion coefficients can be used in the equations to turn the

Computation times for the PDE solvers are almost propor-

This presentation has emphasized issues important for new users of

(45). Further, large comprehensive models can be used to great

View publication stats

You might also like