A Flexible Nonlinear Feedback System That Captures Diverse Patterns of Adaptation and Rebound
 462 Downloads
Abstract
An important approach to modeling tolerance and adaptation employs back mechanisms in which the response to the drug generates a counterregulating action which affects the response. In this paper we analyze a family of nonlinear back models which has recently proved effective in modeling tolerance phenomena such as have been observed with SSRI’s. We use dynamical systems methods to exhibit typical properties of the responsetime course of these nonlinear models, such as overshoot and rebound, establish quantitive bounds and explore how these properties depend on the system and drug parameters. Our analysis is anchored in three specific in vivo data sets which involve different levels of pharmacokinetic complexity. Initial estimates for system (k _{in}, k _{out}, k _{tol} ) and drug (EC _{50}/IC _{50}, E _{max}/I _{max}, n ) parameters are obtained on the basis of specific properties of the responsetime course, identified in the context of exploratory (graphical) data analysis. Our analysis and the application of its results to the three concrete examples demonstrates the flexibility and potential of this family of back models.
Key words
adaptation back pharmacokineticpharmacodynamic modeling rebound system biology system dynamics toleranceINTRODUCTION
Clinical tolerance may be associated with the effects of many drugs especially opioids, various central nervous stimulants, and organic nitrates (cf. Goodman and Gilman 1996 (1)). Acute tolerance is also seen in an experimental setting in laboratory animals. Factors like dose rates and frequencies, and diurnal oscillations (Sällström et al. 2005 (2) and Visser et al. 2006 (3)) often manifest themselves as tolerance to the pharmacological activity of the test compound. In our own experience we then see how the rate and extent of an intravenous infusion may impact the onset, intensity and duration of e.g., a cardiovascular response. Therefore, we believe that a parametric characterization of the rate and extent of tolerance development has a practical bearing on design and analyses of pharmacodynamic studies.
Over the years there have been several different approaches to modeling tolerance. We mention timedependent attenuation of drug constants (cf. Colburn et al. 1994 (4)), effect compartment models (Porchet et al. 1988 (5) and Ouellet and Pollack 1997 (6)), systems analysis (Urquhart and Li 1968 (7) and 1969 (8), VengPedersen et al. 1993 (9)), pool/precursor models (Licko and Ekblad 1992 (10), Bauer and Fung 1994 (11) and Sharma et al. 1998 (12), and back turnover models (Ackerman et al. 1964 (13), Resigno and Segre 1961, 1966 (14), Holford et al. 1990 (15), Wakelkamp et al. 1996 (16), Agersø et al. 2001 (17), Zuideveld et al. 2001 (18), Lima et al. 2004 (19), Sällström et al. 2005 (2), Visser et al. 2006 (3), and de Winter et al. 2006 (20)).

Rapid drug induced changes in the pharmacodynamic steady state yield pronounced overshoot and rebound.

Inhibition of the loss term which raises the steady state R _{ss}, causes an overshoot which is larger than the rebound upon the return to the baseline after washout.

Stimulation of the loss term which lowers the steady state R _{ss}, causes an overshoot which is smaller than the rebound upon the return to the baseline after washout.
Similar behavior is found when the production term is stimulated or inhibited.
It is the characteristic shape of these response curves: a large overshoot with respect to the new pharmacodynamic steady state when R _{ss} > R _{0} and a large rebound upon return to the baseline when R _{ss} < R _{0}, and the difference in magnitude of these two, which we want to study in this paper. An important factor in determining the size of overshoot or rebound is the rate of input or removal of the drug in relation to the intrinsic rates of the system. We show how a rapid input or removal tends to cause a large overshoot or rebound, whilst a gradual input or removal tends to suppress overshoot or rebound.
 (a)
A data set displaying a pronounced rebound effect of the response when the exposure to the test compound is abruptly terminated. The pharmacokinetics of this test compound has a 1–3 min halflife (Example I).
 (b)
A pharmacodynamic data set derived from constant rate intravenous infusion of an experimental compound displaying both tolerance, drift in the baseline, and no rebound upon return towards baseline (Example II).
 (c)
A recently applied in vivo model for serotonin (5HT) turnover in rat brains during intravenous administration of selective serotonin reuptake inhibitors (SSRI’s) (Bundgaard et al. (21)) (Example III).
It will prove very illuminating to study the evolution of the turnover back sytem in State Space, which for these models is twodimensional, the dimensions being the response R and the moderator M. Viewing the evolution of the system as an Orbit in the (R, M)plane gives valuable insight into the properties of the response versus time curve. Also, we shall see that the structure of the state space immediately reveals some of these properties, such as the tendency of the system towards overshoot when the response jumps to a higher value, whether at onset, when R _{ss} > R _{0} or at washout when R _{ss} < R _{0}.
The plan of the paper is the following. We begin with an analysis of the main class of back models and exhibit a series of characteristic properties such as bounds for R _{max}, for overshoot and for rebound, and analyze how these characteristic properties depend on the different pharmacodynamic rate constants in the model (system properties). Important mathematical ingredients here will be (i) comparison with simpler systems and (ii) asymptotic methods. Then we use these results to discuss successively the data of the three examples exhibited in Fig. 2. We also apply them to the data of Example II to demonstrate how the pharmacodynamic and kinetic constants in the model can be estimated from the response versus time course. In order to gain a deeper understanding of the system, we then return to the back model and present a geometric analysis of the dynamics in state space. We conclude with a discussion of the main results we obtain and the methods we use.
METHODS
In this section we introduce the class of back systems and derive analytically a series of characteristic properties of its response versus time curves. In particular we discuss the propensity of these models to show overshoot and rebound. We first give a detailed description of the class of models. Then we analyze the onset and termination of a constant rate infusion with fast kinetics.
Introduction of the Models
Thus, the production of the response R is here inversely proportional to the strength M of the moderator. In addition, production of response tends to infinity when M drops down to zero and tends to zero when M tends to infinity.
In Bundgaard et al. (21) and in the study of Compound X (cf. Example II) this model was used to fit the data while the drug was assumed to inhibit the loss term (H _{1} = 1 and H _{2} = I). And in Gabrielsson and Peletier (22) this model was fitted to the data from Urquhart and Li (7), also assuming that loss was inhibited by the drug (H _{1} = 1 and H _{2} = I).
 (a)
Overshoot and rebound are biased towards larger response values in that when a drug raises the pharmacodynamic steady state response R _{ss}, then the overshoot tends to be larger than when the drug depresses the steady state response R _{ss}. Similarly, when a return to baseline causes the steady state response R _{ss} to rise, then the overshoot will be larger than when it drops. We shall see that this property can be traced back to the nonlinearity of the system, i.e., to the fact that the moderator acts nonlinearly on the production term. We shall return to this important issue when we give a geometric analysis of the system at the end of the paper.
 (b)
It is possible to obtain precise bounds and estimates for overshoot and rebound. In particular the accuracy of these bounds increases as the dimensionless ratio κ=k _{tol}/k _{out} tends to zero.
Let the departure of the response from the baseline value R _{0} reach its maximum value R _{max} at T _{max}, i.e., R _{max} = R(T _{max}). Then we show that
 (c)R _{max} is either a decreasing or an increasing function of κ. Specifically:

If R _{max} > R _{0} then R _{max} decreases as κ increases.

If R _{max} < R _{0} then R _{max} increases as κ increases.

In the sections which deal with the three examples, the drug regimen and the pharmacokinetics are much more complex, and we shall see how the properties found for the simple plasma concentration profile given by Eq. 9 can still be identified and help in fitting the model to the data.
In the following two subsections we study the impact of the above drug regimen and discuss in succession: (i) onset of a constant rate infusion, assuming that t _{1} = ∞, and (ii) washout at t _{1} < ∞.
Onset of Test Compound
We assume that initially, the system is at baseline, i.e., R(0) = R _{0} and M(0) = M _{0}, and we follow the response R and the moderator M as they evolve with time.
In order to demonstrate the effect of back we follow the response versus time profile as k _{tol} increases from k _{tol} = 0 (no back) to values of k _{tol} which are large compared to k _{out} (strong back).
 Case(i) k _{ tol } « k _{ out } : In this case the moderator changes slowly compared to the response. Thus for an initial period M(t) ≈ M _{0} and hence \(R{\left( t \right)} \approx \overline{R} {\left( t \right)} \). This means that the response quickly rises to a level that is almost equal to R _{top}, which is defined as the root of the equation$$ k_{{in}} \cdot \frac{1} {{M_{0} }}  k_{{out}} \cdot R \cdot I{\left( {C_{{ss}} } \right)} = 0 $$Then, as time progresses, M(t) increases because R > M. However, because of the relatively fast dynamics of the response, any departures from the equalityare quickly corrected. Thus, we may assume that Eq. 15 holds approximately once R has reached its maximum value R _{max} at time T _{max}. When we use this relation to eliminate R from the equation for M in Eq. 8 we obtain$$ k_{{in}} \cdot \frac{1} {{M{\left( t \right)}}}  k_{{out}} \cdot R{\left( t \right)} \cdot I{\left( {C_{{ss}} } \right)} = 0 $$(15)$$ \frac{{dM}} {{dt}} = k_{{tol}} \cdot {\left( {\frac{{R^{2}_{0} }} {{I{\left( {C_{{ss}} } \right)}}} \cdot \frac{1} {M}  M} \right)} $$(16)Since the time of maximal response T _{max} is small and k _{tol} is small, it follows that approximately M(T _{max}) = M _{0}. The solution of equation Eq. 16 with this value at t = T _{max} is given byso that we obtain for the response,$$ M{\left( t \right)} = M_{{ss}} \cdot {\left\{ {1  {\left[ {1  I{\left( {C_{{ss}} } \right)}} \right]} \cdot e^{{  2k_{{tol}} {\left( {t  T_{{\max }} } \right)}}} } \right\}}^{{1/2}} \quad \quad t > T_{{\max }} $$(17)Therefore,$$ R{\left( t \right)} = R_{{ss}} \cdot {\left\{ {1  {\left[ {1  I{\left( {C_{{ss}} } \right)}} \right]} \cdot e^{{  2k_{{tol}} {\left( {t  T_{{\max }} } \right)}}} } \right\}}^{{  1/2}} \quad \quad t > T_{{\max }} $$(18)i.e., in this extreme case, t _{1/2} = ln 2/(2 k _{tol}). Thus, as k _{tol} → 0 it takes longer for the response to settle on its final value R _{ss}. In Fig. 3 this gradual descent is clearly manifest.$$ {\left {R{\left( t \right)}  R_{{ss}} } \right} = O{\left( {e^{{  2k_{{tol}} t}} } \right)}\quad \quad as\quad \quad t \to \infty $$
 Case(ii) k _{ tol } » k _{ out } : Now, the dynamics of the moderator is relatively fast and any departures from the equalityhave soon disappeared. Since R(0) = M(0), we may now assume that Eq. 19 holds approximately for all t ≥ 0. Using Eq. 19 to eliminate M, we find that the equation for R becomes, approximately,$$ k_{{tol}} {\left\{ {R{\left( t \right)}  M{\left( t \right)}} \right\}} = 0 $$(19)$$ \frac{{dR}} {{dt}} = k_{{in}} \frac{1} {R}  k_{{out}} R \cdot I{\left( {C_{{ss}} } \right)} = k_{{out}} I{\left( {C_{{ss}} } \right)} \cdot {\left( {\frac{{R^{2}_{0} }} {{I{\left( {C_{{ss}} } \right)}}} \cdot \frac{1} {R}  R} \right)} $$(20)This equation is of the same form as Eq. 16; its solution, starting from the baseline R _{0}, is given by$$ R{\left( t \right)} = R_{{ss}} \cdot {\left\{ {1  {\left[ {1  I{\left( {C_{{ss}} } \right)}} \right]} \cdot e^{{  2k_{{out}} I{\left( {C_{{ss}} } \right)}t}} } \right\}}^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern\nulldelimiterspace} 2}} $$(21)
Note that R(t) is now a strictly increasing function. Hence, in the limit as k _{tol} → ∞, there is no tolerance. We see this behavior demonstrated in Fig. 4.
We also see thati.e., in this extreme case, t _{1/2} = ln 2/(2 k _{out} I(C _{ ss })).$$ {\left {R{\left( t \right)}  R_{{ss}} } \right} = O{\left( {e^{{  2k_{{out}} I{\left( {C_{{ss}} } \right)}t}} } \right)}\quad \quad as\quad t \to \infty $$It is interesting to note that since R(t)  M(t) → 0 as k _{tol} → ∞, the graph of M(t) is also increasing in the limit as k _{tol} → ∞. However, we see clearly in Fig. 4 that M overshoots its limiting value M _{ss} for k _{tol} = 0.1. In Gabrielsson and Peletier (22) it is shown that there is a range of values of k _{tol} for which R(t) and M(t) tend to, respectively, R _{ss} and M _{ss} through a damped oscillation, whilst outside this range, the approach to R _{ss} and M _{ss} is monotone.
Washout of Test Compound
This equality implies that when k _{tol} is small compared to k _{out}, then overshoot is much larger than rebound.
As with the overshoot at onset, we find that at washout, rebound diminishes monotonically as κ = k _{tol}/k _{out} increases and eventually disappears when κ is large enough.
Conclusion
 1.Uniform upper and lower bounds for the response: If k _{tol} > 0, then$$\begin{array}{*{20}l} {{R_{0} < R{\left( t \right)} < R_{{top}} \quad \quad after{\text{ }}onset} \hfill} \\ {{R_{{bottom}} < R{\left( t \right)} < R_{{ss}} \quad after{\text{ }}washout} \hfill} \\ \end{array} $$(28)
 2.Bounds for overshoot and rebound of the response: If k _{tol} > 0, thenso that for κ small,$$ \begin{aligned} & \left. {\Delta R} \right_{{overshoot}} < \frac{{R_{0} }} {{I(C_{{ss}} )}}{\left[ {1  {\sqrt {I{\left( {C_{{ss}} } \right)}} }} \right]}\quad and\quad \\ & \left. {\Delta R} \right_{{rebound}} < R_{0} {\left[ {1  {\sqrt {I{\left( {C_{{ss}} } \right)}} }} \right]} \\ \end{aligned} $$(29)$$ \left. {\Delta R} \right_{{overshoot}} \approx \frac{1} {{I{\left( {C_{{ss}} } \right)}}}\left. { \cdot \Delta R} \right_{{rebound}} $$(30)
 3.
Pharmacodynamic steady state and baseline are “global attractors”: If k _{tol} > 0 and the infusion is constant for all time, then R(t) → R _{ss} as t → ∞. Similarly, after washout when C = 0 for all time, then R(t) → R _{0} as t → ∞.
 4.Influence of the ratio κ = k _{tol}/k _{out} : We have exhibited the pivotal importance of the ratio κ of the rate constants of the two equations: k _{tol} and k _{out}. For instance,
 (a)
Overshoot and rebound are greatest when κ is small.
 (b)
As κ increases, both overshoot and rebound become smaller and eventually disappear.
 (a)
For further properties we refer to Gabrielsson and Peletier (22).
GENERALIZATIONS
The analysis we have presented for the back system 8 in which the loss term is inhibited by the drug, and the results we have obtained for that system, easily carry over to the other three types of systems defined in Eq. 2, i.e., the systems in which the loss term is stimulated (H _{1} = 1 and H _{2} = S) or the production term is either stimulated (H _{1} = S and H _{2} = 1) or inhibited (H _{1} = I and H _{2} = 1).
Many of the properties of the response versus time curves which we have sketched above can be explained in a geometrical way by viewing the state (R,M) of the system as moving along a trajectory in state space. We present such an explanation in a special section towards the end of the paper.
EXAMPLE I: INHIBITION OF PRODUCTION OF RESPONSE
Parameter estimates for Example I
R _{0}  k _{out}  k _{tol}  n  IC_{50}  κ 

0.58  0.59  0.058  1.4  0.13  0.1 
3%  36%  25%  7%  26%  – 
The exact start and stop of drug input (infusions) and the exact sampling times of the pharmacological response are particularly important for this system with rapid turnover of both test compound and response. Since the half life of kinetics falls in the range of 1–2 min and half life of response is about 1 min, a 15–30 sec error in sampling times at pivotal occasions will not only impact parameter estimates (accuracy) but also their precision. We will expand on this particular system (response) and test compound in a subsequent series of papers.
EXAMPLE II: INHIBITION OF LOSS OR STIMULATION OF PRODUCTION OF RESPONSE
Parameter estimates for models (A) and (B)
 Estimate (A)  CV % (A)  Estimate (B)  CV % (B) 

k _{tol}  0.057  10  0.055  10 
k _{out}  0.26  5  0.13  6 
IC _{50}/SC _{50}  9.1  6  15  6 
I _{max}/S _{max}  0.51  2  1.0  3 
n  1.5  6  1.5  6 
k _{in}(0)  5,773  5  2200  5 
α  0.00039  17  0.00040  13 
κ  0.22  –  0.42  – 
This suggests that more doses and eventually an extended input (infusion) period are necessary challenges for model discrimination. It was necessary to model the drift in baseline as a linear function, as part of the full model. We would not advocate modeling baseline subtracted response data because the meaning of different parameters would become less transparent.
In future experiments we would suggest a separate control group that captures the unperturbed baseline drift.
An ordinary turnover model (by inhibition of loss), or a pool/precursor model (Stimulation of the loss from the pool compartment) did not adequately capture the upswing, overshoot, roof and decline to the baseline of this data set. For the turnover model, this can be explained by analytical arguments (cf. (25)).
EXAMPLE III: INHIBITION OF LOSS OF RESPONSE
Parameter estimates for the model Eq. 38
R _{0}  k _{out}  k _{tol}  I _{max}  IC _{50}  n  κ 

105  0.18  0.003  0.84  4.1  0.87  0.02 
5%  11%  18%  2%  25%  11%  — 
Critical response levels in Examples I–III
Example  R _{0}  R _{ss}  R _{top}  R _{bottom} 

I  \({\sqrt {\frac{{k_{{in}} }} {{k_{{out}} }}} } \)  \(R_{0} {\sqrt {I{\left( {C_{{ss}} } \right)}} } \)  \(\frac{{R_{0} }} {{{\sqrt {I{\left( {C_{{ss}} } \right)}} }}} \)  \(R_{0} I{\left( {C_{{ss}} } \right)} \) 
II(A)  \( {\sqrt {\frac{{k_{{in}} }} {{k_{{out}} }}} } \)  \( \frac{{R_{0} }} {{{\sqrt {I{\left( {C_{{ss}} } \right)}} }}} \)  \( \frac{{R_{0} }} {{I{\left( {C_{{ss}} } \right)}}} \)  \( R_{0} {\sqrt {I{\left( {C_{{ss}} } \right)}} } \) 
II(B)  \( {\sqrt {\frac{{k_{{in}} }} {{k_{{out}} }}} } \)  \( R_{0} {\sqrt {S{\left( {C_{{ss}} } \right)}} } \)  \( R_{0} S{\left( {C_{{ss}} } \right)} \)  \( \frac{{R_{0} }} {{{\sqrt {S{\left( {C_{{ss}} } \right)}} }}} \) 
III  \( {\sqrt {\frac{{k_{{in}} }} {{k_{{out}} }}} } \)  \( \frac{{R_{0} }} {{{\sqrt {I{\left( {C_{{ss}} } \right)}} }}} \)  \(\frac{{R_{0} }} {{I{\left( {C_{{ss}} } \right)}}} \)  \( R_{0} {\sqrt {I{\left( {C_{{ss}} } \right)}} } \) 
By extending the original Bundgaard et al. (21) analysis to also include a transduction step for the moderator compartment, the parameter precision of k _{tol} was improved. Increasing k _{tol} also resulted in a shortening of its corresponding halflife (ln 2/k _{tol}). All other parameters and their precision remained the same.
This demonstrated the flexibility of this back turnover model without increasing the number of model parameters. We will expand on this topic in a subsequent series of papers.
It is not surprising that the model predicts an uptake halflife (t _{1/2kout} ) of about 3–4 min in light of the rate limiting distribution across the dialysis probe membrane. The true uptake halflife into presynaptic vesicles does probably lie in the microsecond to second range, but the equilibrium process across the dialysis probe confounds the more rapid reuptake process.
INITIAL PARAMETER ESTIMATES
 (a)The evolution of the baseline over time, modeled by$$ R_{0} {\left( t \right)} = R_{0} {\left( 0 \right)}{\left( {1 + \alpha t} \right)},\quad \quad R_{0} {\left( 0 \right)} = 149,\quad \quad R_{0} {\left( {250} \right)} = 162 $$(40)
 (b)The evolution of the pharmacodynamic steady state modeled by$$ R_{{ss}} {\left( t \right)} = R_{{ss}} {\left( 0 \right)}{\left( {1 + \alpha t} \right)},\quad \quad R_{{ss}} {\left( 0 \right)} = 220 $$(41)
 (c)The slope at onset of infusion (t = t _{0} = 32):$$ slope{\mathop = \limits^{def} }\left. {\frac{{dR}} {{dt}}} \right_{{t = t_{0} }} = 17 $$(42)
For convenience we collect some pertinent formulae in Table IV.
THE STATE SPACE—A GEOMETRIC DESCRIPTION OF THE DYNAMICS
Below we list a few observations about these orbits.
Onset of response
Since I(C _{ss}) < 1, the null cline Γ _{ R } shifts upward.
k _{tol} = 0: Since M(t) = M _{0} along this orbit, it traces out a vertical line until it reaches Γ _{ R } at R=R _{top}.
k _{tol} << k _{out}: Here M(t) grows slowly and the orbit is almost vertical until it reaches Γ _{ R }, crosses it at time T _{max}, where the response reaches its maximal value R _{max}. The orbit then follows Γ _{ R } down until it reaches the pharmacodynamic steady state (R _{ss}, M _{ss}). Note that R _{max} ≈ R _{top}.
k _{tol} ≈ k _{out}: Here R(t) and M(t) develop more in tandem and the overshoot is much diminished, i.e. R _{max} is much lower. Notice the spiraltype behavior near (R _{ss},M _{ss}).
Washout of response
Since I(0) = 1, the null cline Γ _{ R } drops back to the starting position.
k _{tol} = 0: Since M(t) = M _{ss} along this orbit, it traces out a vertical line until it reaches Γ _{ R } at R = R _{bottom}
k _{tol} << k _{out}: The orbit drops down and quickly reaches Γ _{ R }. Thereafter, it slowly follows Γ _{ R } back towards the baseline (R _{0},M _{0}).
k _{tol} ≈ k _{out}: The rebound is now less pronounced and the baseline point (R _{0},M _{0}) is approached via a small spiral.
 (i)
R _{max} cannot exceed R _{top} and that R _{min} cannot drop below R _{bottom}.
 (ii)
Overshoot (R _{max}  R _{ss}) as well as rebound (R _{0}  R _{min}) become smaller as the ratio κ = k _{tol}/ k _{out} becomes larger.
The null clines in the (R,M) plane are now given by \(\Gamma _{R} = {\left\{ {{\left( {R,M} \right)}:M = {M_{0} } \mathord{\left/ {\vphantom {{M_{0} } {\text{I}}}} \right. \kern\nulldelimiterspace} {\text{I}}{\left( {C_{{ss}} } \right)}\,} \right\}}\,\,\,\,and\,\,\,\,\Gamma _{M} = {\left\{ {{\left( {R,M} \right)}:R = M} \right\}}\)
After onset, Γ _{ R } shifts to the right since I(C _{ss}) < 1 if C _{ss} > 0, and it shifts back at washout. We see how the orbit, which starts at the baseline (R _{0}, M _{0}) = (1,1) spirals towards (R _{ss}, M _{ss}) = (2, 2), as predicted by Eq. 52 and after washout spirals back to the baseline.
DISCUSSION
Background
From a drug discovery point of view many biological variables including the ones affected by drugs, are subject to adaptation and tolerance. In the description of the time course of the response to drugs, it may therefore be necessary to consider the endogenous back control as an intrinsic part of the pharmacodynamics. With still the same system parameters and properties, different drug parameters will result in different responsetime courses and consequently different levels of onset, intensity, duration and adaptation.
Does tolerance confound the estimability of system and drug parameters? Yes, particularly with slow drug kinetics relative to kinetics of response, and with a substantial baseline drift, unless data are rich and wellspaced such as in Examples II and III. Since tolerance/adaptation may develop on a totally different (slower) time scale than the primary response (k _{out}), this may result in unexpected (and unpredictable from acute dosing) responses at pharmacodynamic steadystate. We believe that not only tolerance per se, but tolerance coupled to a substantial rebound effect may jeopardize the clinical outcome of unscheduled drug holidays (missed doses) and temporarily result in clinically harmful rebound effects. Therefore, multiple dosing and continuous tracking of the pharmacodynamic effect will be necessary at two or more dose levels. The kappa parameter (k _{tol}/k _{out}) may also serve as an early indicator variable of how slow tolerance develops and therefore how long extended exposure is needed.
In quantitative pharmacology, the physiology and the system need to be (i) condensed into a conceptual (picture) model that captures the main features, including rate limiting step(s). When the mechanism of action has been presented as a conceptual model, the latter then needs to be translated into (ii) equations, i.e., the mathematical model, e.g., of the rate process (quantitative pharmacological), which is implemented into the regression software package. This is the second condensation, because model parameters also need to be identifiable and estimable.
The next step (iii) will be to add the statistical model (error and/or covariate models) to the regression model. This is, of course, of importance in most instances, but it also needs to be balanced against mechanism of action (including the mathematical structure) and supported by a thorough experimental design. We envision more of a focus on mechanistic/mathematical model building including emphasis on experimental design. However, pattern recognition and exploratory data analysis are more related to the art of successful model building than the pure mathematical/statistical regression.
Examples
The analysis in this paper has been anchored in data sets obtained from three examples.
Example I covers rapid test compound kinetics, multiple provocations, washout dynamics, and a value of κ = k _{tol}/k _{out} which is less than 1. These combinations reveal the system properties and result in accurate and precise system and drug parameters. The pharmacological response exhibits adaptation and particularly extensive rebound.
Example II covers rapid system kinetics, relatively slow washout of test compound and a substantial drift in the baseline. Again, the value of κ is less than 1. The pharmacological response shows a pronounced overshoot, adaptation and a slow return towards baseline without any indication of a rebound. However, simulations also demonstrate rebound of this system provided rapid kinetics and stable baseline.
Example III covers intravenous infusion at a 4fold dose (exposure) range, slow washout of test compound and a stable baseline. The pharmacological response shows a rapid onset, superimposing responsetime courses for the two highest doses and lack of rebound upon return towards the baseline. This time the value of κ was much less than 1, but the slow pharmacokinetics of the test compound confounded the possibility of a rebound.
The intravenous way of administration avoids confounding absorption processes and firstpass phenomena. In neither of the three examples, pharmacologically active metabolites are known to be involved.
Initial Estimates
Using the data set of Example II we have shown how initial estimates for system and drug parameters can be obtained from just a few critical properties of the response versus time curve, such as baseline, pharmacodynamic steady state and initial slope.
Qualitative Analysis of the Model
Using qualitative methods we have obtained apriori estimates for the response versus time course, such as upper and lower bounds, rate of convergence to the pharmacodynamic steady state, and to the baseline and insight into the crucial importance of the ratio κ of the rate constants for the moderator k _{tol} and the response k _{out}. Viewing the dynamics of the system as an orbit in the state space represented by the (R, M)—plane has yielded valuable insight into the specific properties of the back model such as a tendency towards overshoot and rebound, and an “upward bias” in that the tendency towards overshoot tends to be larger if the response jumps up than when it jumps down.
With the apriori knowledge we have obtained for this class of back models it will be easier to assess before hand whether or not a given data set can be fruitfully modeled by a model of this type.
Design
What are study designs generally lacking in situations where we also have back, tolerance and adaptation? We advocate several dose levels coupled to different routes and rates of input, multiple dosing, redosing during overshoot, pharmacodynamic steadystate or during rebound. Multiple consecutive intravenous infusions were given to the system presented in Example I. This is an example of multiple perturbations of a system, which reveals the system properties nicely. Since adaptation is a timephenomenon, experimental design also requires repeated provocations (repeated dosing, different rates and routes of administration) of the system that is being studied. There is no generic schedule for this, but multiple provocations for at least 3–4 t _{1/2ktol}, coupled to washout dynamics with sufficient baseline data, are probably informative. In a subsequent series of papers we will demonstrate how multiple rates and durations of intravenous dosing, coupled to observations of pharmacology during and after the provocations (washout), may help to elucidate the rate limiting steps behind a physiological turnover system exhibiting rapid adaptation.
We also advocate ex vivo plasma protein binding of the test compound for appropriate estimation of unbound exposure, since changes in plasma protein binding have been shown to affect the exponent (Hill coefficient) in the sigmoid E _{max} model when converting total to unbound concentrations. Ex vivo plasma protein binding has, in our experience, given different results compared to pooled in vitro plasma protein binding experiments. One may also need to consider dosing (infuse) potentially active metabolites if their potencies and/or exposure levels fall in the neighborhood of parent compound. Making model simulations prior to executing the full experiment also helps to set up a strategy for experimental design.
A constant (stable) drug input is probably the least satisfactory design for a back experiment. For tolerant physiological systems we believe that the constant 24 hr exposure paradigm may be deleterious in that most physiological systems need to be reset to some extent at some point during each dosing interval for the system to be susceptible to a new dose.
Final Conclusion
We have demonstrated the flexibility of a class of turnover models in their ability to mimic overshoot, pharmacodynamic steady state, return to baseline and rebound. We anchored this analysis to three experimental data sets of varying complexity. Finally, we discuss analytical tools such as phase plane plots and show a strategy as to how to derive initial parameter estimates.
Notes
Acknowledgement
The second author wishes to record his gratitude to the Center for BioDynamics at Boston University for the hospitality enjoyed during the writing of this review.
References
 1.L. S. Goodman, and A. G. Gilman. The Pharmacological Basis of Therapeutics, 9th Ed., McGrawHill, New York, 1996.
 2.B. Sällström, S. A. G. Visser, T. Forsberg, L. A. Peletier, A. C. Ericson, and J. Gabrielsson. A novel pharmacodynamic turnover model capturing asymmetric circadian baselines of body temperature, heart rate and blood pressure in rats: challenges in terms of tolerance and animal handling effects. Journal of Pharmacokinetics and Pharmacodynamics. 32:835–859 (2005).CrossRef
 3.S. A. G. Visser, B. Sällström, T. Forsberg, L. A. Peletier, and J. Gabrielsson. Modeling drug and system related changes in body temperature: application to Clometiazoleinduced hypotyhermia, longlasting tolerance developments and circadian rhythm in rats. J. Pharmacol. Exp. Ther. 317:209–219 (2006).CrossRef
 4.W. A. Colburn, and M. A. Eldon. Simultaneous Pharmacokinetic/Pharmacodynamic modeling. In N. R. Butler, J. J. Sramek and P. K. Narang (eds.), Pharmacodynamics and Drug Development: Perspectives in Clinical Pharmacology, Wiley, New York 1994.
 5.H. C. Porchet, N. L. Benowitz, and L. B. Sheiner. Pharmacodynamic model of tolerance: Application to nicotine. J. Pharmacol. Exp. Ther. 244:231–236 (1988).
 6.D. M. C. Ouellet, and G. M. Pollack. Pharmacodynamics and tolerance development during multiple administrations in rats. J. Pharmacol. Exp. Ther. 281:713–720 (1997).
 7.J. Urquhart, and C. C. Li. The dynamics of adrenocortical secretion. Am. J. Physiol. 214:73–85 (1968).
 8.J. Urquhart, and C. C. Li. Dynamic testing and modeling of andrecortical secretory function. Ann. NY Acad. Sci. 156:756–778 (1969).CrossRef
 9.P. VengPederson, and N. B. Modi. A system approach to pharmacodynamics. Inputeffect control system analysis of central nervous system effect of alfentanil. J. Pharm. Sci. 82:266–272 (1993).CrossRef
 10.V. Licko, and E. B. Ekblad. Dynamics of a metabolic system: what singleaction agents reveal about acid secretion. Am. J. Physiol. Gastrointest. Liver Physiol. 262:G581–G592 (1992).
 11.J. A. Bauer, and H. L. Fung. Pharmacodynamic models of nitroglycerininduced hemodynamic tolerance in experimental heart failure. Pharm. Res. 11:816–823 (1994).CrossRef
 12.A. Sharma, W. F. Ebling, and W. J. Jusko. Precursordependent indirect pharmacodynamic response model for tolerance and rebound phenomena. J. Pharm. Sci. 87:1577–1584 (1998).CrossRef
 13.E. Ackerman, J. W. Rosevear, and W. F. McGuckin. A mathematical model of the glucosetolerance test. Phys. Med. Biol. 9:203–213 (1964).CrossRef
 14.A. Rescigno, and G. Segre. Drug and tracer kinetics. Blaisdell Publishing Company, London p. 1961 (1966).
 15.N. H. G. Holford, J. Gabrielsson, L. B. Sheiner, N. Benowitz, and R. Jones. A Physiological Pharmacodynamic Model for Tolerance to Cocaine Effects on Systolic Blood Pressure, Heart Rate and Euphoria in Human Volunteers. (Poster presented at the Symposium on Measurement and Kinetics of In Vivo Drug Effects: Advances in Simultaneous Pharmacokinetic/Pharmacodynamic Modeling. Noordwijk, The Netherlands, 28–30 June, 1990).
 16.M. Wakelkamp, G. Alvan, J. Gabrielsson, and G. Paintaud. Pharmacodynamic modeling of furosemide tolerance after multiple intravenous administration. Clin. Pharmacol. Ther. 60:75–88 (1996).CrossRef
 17.H. Agersø, L. Ynddal, B. Søgaard, and M. Zdravkovic. Pharmacokinetic and pharmacodynamic modeling of NN703, a growth hormone Secretagogue, after a single po dose to human volunteers. J. Clin. Pharmacol. 41:163–169 (2001).CrossRef
 18.K. P. Zuideveld, H. J. Maas, N. Treijtel, J. Hulshof, P. H. Van der Graaf, L. A. Peletier, and M. Danhof. A setpoint model with oscillatory behaviour predicts the time course of 8OHDPATinduced hypothermia. Am. J. Physiol. Regul. Integr. Comp. Physiol. 281:R2059–R2071 (2001).
 19.J. J. Lima, N. Matsushima, N. Kissoon, J. Wang, J. E. Sylvester, and W. J. Jusko. Modeling the metabolic effects of terbutaline in β_{2} adrenergic receptor diplotypes. Clin. Pharmacol. Ther. 76:27–37 (2004).CrossRef
 20.W. De Winter, J. De Jongh, T. Post, B. Ploeger, R. Urquhart, I. Moules, D. Eckland, and M. Danhof. A mechanismbased disease progression model for comparison of longterm effects of Pioglitazone, Metformin and Gliclazide in disease processes underlying Type 2 Diabetes Mellitus. Journal of Pharmacokinetics and Pharmacodynamics. 33:313–343 (2006).CrossRef
 21.C. Bundgaard, F. Larsen, M. Jørgensen, and J. Gabrielsson. Mechanistic turnover model including autoinhibitory back regulation of acute 5HT release in rat brain after administration of selective serotonoin reuptake inhibitors (SSRI’s). Eur. J. Pharm. Sci. 29:394–404 (2006).CrossRef
 22.J. Gabrielsson, and L. A. Peletier. A nonlinear back model capturing different patterns of tolerance and rebound. Eur. J. Pharm. Sci. 32:85–104 (2007).CrossRef
 23.N. L. Dayneka, V. Garg, and W. J. Jusko. Comparison of four basic models of indirect pharmacodynamic responses. J. Pharmacokinet. Biopharm. 21:457–478 (1993).CrossRef
 24.P. Blanchard, R. L. Devaney, and G. R. Hall. Differential Equations, Brooks/Cole Publishing Co, Pacific Grove, CA, 1998.
 25.L. A. Peletier, J. Gabrielsson, and J. Den Haag. A dynamical systems analysis of the indirect response model with special emphasis on time to peak response. Journal of Pharmacokinetics and Pharmacodynamics. 32:607–654 (2005).CrossRef
 26.J. Gabrielsson, and Weiner, D. Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications, 4th Ed., Swedish Pharmaceutical Press, Stockholm. ISBN 13 978 91 9765 1004, 2006.