Impact of saturable distribution in compartmental PK models: dynamics and practical use

We explore the impact of saturable distribution over the central and the peripheral compartment in pharmacokinetic models, whilst assuming that back flow into the central compartiment is linear. Using simulations and analytical methods we demonstrate characteristic tell-tale differences in plasma concentration profiles of saturable versus linear distribution models, which can serve as a guide to their practical applicability. For two extreme cases, relating to (i) the size of the peripheral compartment with respect to the central compartment and (ii) the magnitude of the back flow as related to direct elimination from the central compartment, we derive explicit approximations which make it possible to give quantitative estimates of parameters. In three appendices we give detailed explanations of how these estimates are derived. They demonstrate how singular perturbation methods can be successfully employed to gain insight in the dynamics of multi-compartment pharmacokinetic models. These appendices are also intended to serve as an introductory tutorial to these ideas.


Introduction
In practical applications, population pharmacokinetic modellers are regularly confronted with data suggesting nonlinear kinetics of the investigational compound. This may include disproportionate increases in C max in single ascending dose (SAD) data or disproportionate accumulation in multiple ascending dose (MAD) data. Such nonlinearities may be difficult to account for using the standard linear compartmental pharmacokinetic (PK) model, even when nonlinear elimination is employed. Here we investigate a class of compartmental PK models which can be characterized as saturable distribution models, which we feel can provide an additional tool enabling pharmacometric modelers to tackle observed nonlinearities in their data.
Compartmental PK models usually combine a central or plasma compartment, which represents the site at which pharmacokinetic sampling takes place, with one or more peripheral or tissue compartments. Such multi-compartmental models typically assume that drug enters the blood stream in the central compartment, is distributed from there via linear first order processes to the peripheral compartments, and finally is eliminated again from the central compartment via either a linear first order process or a saturable Michaelis-Menten process (see e.g. Wagner et al. [1] and more recently, Wu et al. [2], Brocks et al. [3] and Scheerens et al. [4]). While linear distribution from central to peripheral may often provide an adequate description of the observed PK, very few processes in biology are truly linear. Most, if not all biological processes are saturable and may only appear linear because their maximum capacity has not been approached in the observed data. It follows that the standard multi-compartmental PK model with linear distribution can be seen as a special case of a & Lambertus A. Peletier peletier@math.leidenuniv.nl 1 more general class of multi-compartmental PK models with saturable distribution. Snoeck et al. [5] first developed a population PK model with saturable distribution to account for the nonlinear PK of draflazine. This nonlinearity was found to be related to a capacity-limited, high-affinity binding of draflazine to nucleoside transporters located on erythrocytes and endothelial tissue, and could not be accounted for by conventional, linear distribution PK models. In the model developed by Snoeck et al., draflazine was distributed from a central compartment with linear elimination to three peripheral compartments, two of which were capacity-limited with different capacities but similar affinity and were thought to represent the specific binding of draflazine to its receptors on erythrocytes and tissue, respectively. This model was found to satisfactorily predict the nonlinear, dose-dependent PK of draflazine and its disposition in whole blood and plasma.
In an unpublished study, the approach developed by Snoeck et al. was used to model the PK of compound X, which also showed a markedly nonlinear PK and was also known to bind specifically to receptors on the erythrocytes.
Starting from a conventional three compartment PK model, transformation of one of the two peripheral compartments to a low capacity, high affinity compartment with saturable distribution resulted in a highly significant improvement of the model fit. This compartment was thought to represent specific binding to the receptors on the erythrocytes, and addressed a nonlinear dose-dependent increase of C max observable in single ascending dose (SAD) studies. However, Fig. 1 shows that this 1-receptor model still failed to address nonlinear dose-dependencies in both accumulation and time to steady-state in multiple ascending dose (MAD) studies. Transformation of the second peripheral compartment to a very high capacity, low affinity compartment with saturable distribution addressed this problem and yielded a further, highly significant improvement of the model fit. This 2-receptor saturable distribution model was used to develop a successful individual dose titration protocol, and was mathematically analysed by Peletier et al. [6].
What kind of non-linearities in the observed PK can be addressed by saturable distribution models, when and how should we apply them? In the following we address such questions by exploring the dynamics of a two-compartmental model with a saturable, Michaelis-Menten type rate function for the distribution of drug from the central to the peripheral compartment. We do this for two opposing variants of saturable distribution: first, we explore the dynamics of a model with a low affinity, high capacity distribution process, and then discuss the dynamics of a model with high affinity, low capacity distribution. In order to assess the impact of saturation, we analyse the dynamics of two classes of models: one with linear and one with saturable distribution.
The objectives of this paper are (i) To identify characteristic properties of the time courses in the central compartment, and identify differences between linear and saturable models which may serve as handles to determine which class of models should be used to fit a given set of data. (ii) To study the dynamics of the nonlinear model incorporating saturation with a view to understand the impact of the relative capacities and the rate constants of the system and identify the characteristic time-scales. (iii) To identify the impact of saturable distribution in practical applications, such as the exposure resulting from SAD and MAD regimens.
The mathematical analysis that is used to prove the results in this paper is presented in three appendices. The first two are devoted to the large capacity case, the linear model and the saturable model, and the third appendix is devoted to the small capacity case. The analysis relies strongly on applications of singular perturbation theory (cf. [7,8]). The appendices are written so that they can be used as a tutorial for applications of this method in pharmacokinetics and pharmacodynamics.

Methods
In order to study the impact of saturation we compare the dynamics of two distribution models, one with linear and one with nonlinear distribution that involves saturation. In both models a test compound or drug is supplied to an absorption space (1). The drug is then discharged into a central compartment (2), distributed over a peripheral compartment (3), as well as eliminated from the central compartment.

Linear distribution model
This is the standard linear two-compartment distribution model in which drug flows between the central compartment (1) and the peripheral compartment by diffusion in which the flux is proportional to the difference of the concentrations in the two compartments. The amount of drug in the absorption space is denoted by A 1 and the concentrations in the central and the peripheral compartment are denoted by, respectively, C 2 and C 3 . These quantities satisfy the following system of differential equations: Here q denotes the infusion rate, k a a first order rate constant, Cl the non-specific clearance, Cl d the intercompartmental distribution and V 2 and V 3 the volumes of the central-and the peripheral compartment.
In comparing this linear distribution model to the nonlinear model involving saturation below, it is convenient to use the amount of drug in the central compartment (A 2 ¼ V 2 Â C 2 ) and in the peripheral compartment (A 3 ¼ V 3 Â C 3 ). Introducing these amounts into the system (1) then results in the following system of differential equations: where and H is a dimensionless constant which can be viewed as a measure of the ''relative capacity'' of the central and the peripheral compartment.

Nonlinear or saturable distribution model
In this model the transfer from the central compartment to the peripheral compartment is saturable, whilst that from the peripheral back to the central compartment is linear.
Specifically we study the model where q; k a ; k 20 and k p are as in the linear problem. Here B max is referred to as the capacity of the peripheral compartment and K M the Michaelis-Menten constant. Both B max and K M have the dimension of an amount. Thus, saturation is modelled by a Michaelis-Menten term which involves two new parameters, the capacity B max and K M . This model has five parameters whereas the linear model has four.
Remark For values of A 2 which are small relative to K M , the Michaelis-Menten term in the nonlinear system may be approximated by ðB max =K M Þk p A 2 . Thus the relative capacity H in the linear system may be compared to the quotient B max =K M in the nonlinear system.
In the large capacity case, the infusion rate q is assumed to be constant, and initially the system is assumed to be empty, i.e., the amounts in the compartments are all assumed to be zero: In the small capacity case, the infusion rate q is assumed to be zero, and the initial conditions after an iv dose D are given by

Steady state
For reference we give here the steady state values of A 1 ; A 2 and A 3 when A 1 is supplied to the absorption space at a constant rate k f ðtÞ q. Equating the temporal derivatives in Eqs. (2) and (3) to zero we obtain the following expressions for the steady state amounts A i;ss (i ¼ 1; 2; 3): Thus, we can write A 3;ss in terms of A 2;ss : ;ss ðLinearÞ and We conclude that in both models A 1;ss and A 2;ss are the same and increase linearly with the infusion rate q. In the linear model the amount A 3;ss in the peripheral compartment also increases linearly with q, but in the nonlinear model it increases nonlinearly and converges to the capacity B max as the infusion rate tends to infinity: We shall see however that whereas in the linear model the time needed for A 2 ðtÞ to reach steady state is independent of q, in the nonlinear model it varies with the infusion rate. Evidently, in the absence of an infusion rate, i.e., when q ¼ 0, the steady state is given by We contrast the dynamics of models with large capacity peripheral compartment, combined with slow transfer with models with small capacity peripheral compartments endowed with rapid transfer.

Large capacity and slow distribution
We assume, A.1 The capacity of the peripheral compartment is large compared to that of the central compartment.
A.2 The drug flows back from the peripheral compartment into the central compartment at a much smaller rate than it is eliminated from the central compartment. Specifically, in terms of the rate constants we assume that: Small capacity and rapid distribution We assume, A.3 The capacity of the peripheral compartment is small compared to that of the central compartment.
A.4 Elimination from the central compartment is much slower than the rate with which the drug flows back into the central compartment.

Simulations
In order to acquire a qualitative understanding of the structure of the dynamics of both models, given the relative magnitudes of the rate constants k a , k 20 and k p , and the capacity of the peripheral compartment of the linear model (H) and the nonlinear model (B max ), we perform a series of simulations. We do this separately for the large and the small capacity peripheral compartment.

Large capacity and slow distribution
We select a series of different values of the infusion rate q in order to demonstrate the differences between the linear and the nonlinear model. These simulations will then be done for the following parameter values: Because of the large value of k a , the compound in the absorption space very quickly reaches a quasi-steady state so that we may put A 1 ðtÞ ¼ A 1;ss ¼ q=k a for t [ 0. Thus, the dynamics of the system is effectively determined by the interaction between the central and the peripheral compartment.
In Figs. 2 and 3 we show how in the linear and the nonlinear model the amount of compound in the central compartment (A 2 ) evolves with time for the different infusion rates. The simulations for the linear and the nonlinear system look similar. Both exhibit a clear two-phase structure, which can be divided into: A brief initial phase in which A 2 climbs to what appears to be a plateau. We shall refer to this value of the amount of compound as the Plateau value and denote it by A 2 .  (2) graphs of A 2 ðtÞ for increasing infusion rates q ¼ 1, 2, 3, 4, 5 mg h À1 when k a ¼ 10, k 20 ¼ 10 À2 , k p ¼ 10 À4 h À1 and H ¼ 100 A second, much longer phase in which the final plateau value A 2 of the first phase serves as a starting point of a slow rise towards the final limit which, as expected, is the steady-state value A 2;ss .
However, Figs. 2 and 3 demonstrate that the impact of the infusion rate q is very different. Here we focus on how the infusion rate q affects the following characteristics of the dynamics: (1) The plateau value A 2 after the first phase.
(2) The half-life of the convergence to the plateau value A 2 as well as the half-life of the convergence to the final steady state value A 2;ss .
As can be expected from a linear problem, we see in Fig. 2 and Eq. (6) that A 2 and A 2;ss depend linearly on q and that the half lives in the two phases are independent of the infusion rate. The simulations in Fig. 3 demonstrate that for the nonlinear model the influence of the infusion rate q is more complex. However, the terminal state A 2;ss is the same as for the linear problem (cf. (6)) and hence depends linearly on the infusion rate: Thus, in comparing the two models one needs to focus on the complete temporal profile i.e., the concentration versus time profile for all time. We make the following observations: • The plateau value, A 2 , increases with increasing q. For the linear model A 2 is seen to increase linearly with q (cf. Fig. 2) whilst for the nonlinear model the dependence on q appears to be super-linear, i.e., A 2 appears to grow faster than linearly with q (cf. Fig. 3). • The half-life in the two phases. As the infusion rate q increases, the half-life in the first phase appears to increase whilst the half-life in the second phase appears to decrease.

Small capacity and rapid distribution
In Fig. 4 we present a series of simulations for nonlinear, saturable distribution model which exhibit the impact of an iv bolus dose on the initial peak of A 2 ðtÞ. The doses and the parameter values are given in Table 2 in Appendix 2:  (3) graphs of A 2 versus time for the parameter values k a ¼ 10, It is seen that in this case the disposition also has a twophase structure: soon after administration, A 2 ðtÞ jumps up to a high Peak value A 2;max quickly drops thereafter (left figure) and then, in a second phase slowly returns to zero (right figure). This Peak value is seen to increase rapidly with D in a super-linear manner: when D increases from 10 to 20 mg then A 2;max rises by about 4 mg and when D increases from 60 to 70 the rise is about 7 mg, i.e., almost double the low-dose increase.
Thus, as in the large capacity case, the graphs of A 2 versus time exhibit a two-phase structure, albeit with a completely different shape. A brief initial phase, say for 0\t\t 0 , in which A 2 ðtÞ exhibits a violent up-and down swing which ends with A 2 at an intermediate plateau value A 2 , followed by a much longer elimination phase.

Results
Many of the observations made in the simulations can be explained through mathematical analysis of the linear twocompartment model (2) and the nonlinear model (3). Below we present a series of results from such analysis. We discuss the large capacity and the small capacity case in succession.

Large capacity and slow distribution
At first sight the simulations in Figs. 2 and 3 for the two models are qualitatively similar: a rapid rise of A 2 towards an intermediate plateau A 2 , the plateau value, followed by a slow rise towards the final steady state A 2;ss given in Eq. (6). In order to discriminate between the dynamics of the linear and the nonlinear model it is therefore important to obtain detailed and quantitative information about characteristics of the dynamics over time. We focus here on two such characteristic properties: -The intermediate plateau value A 2 , and -The half-life of the convergence as A 2 tends to A 2 , and as A 2 tends to A 2;ss and the way these quantities depend on the infusion rate, the capacity and the different rate constants.
For both models we present such quantitative estimates of the plateau value and the half-life in the first and the second phase. Their proofs are given in the mathematical analysis presented in Appendices 1 and 2.

Plateau value
The existence of a plateau value is a result of the two-phase structure of the dynamics of this system in which two different time scales can be distinguished: 1 In light of the basic assumption (9) there is a significant difference between these two time scales. For the parameter values of Table 1 the half-life of the first phase is about a factor 100 shorter than that of the second phase. During the first phase, return flow from the peripheral compartment is still negligible because k p is very small and A 3 is still building up. Therefore, during this phase the term k p A 3 modelling the back flow from the peripheral compartment into the central compartment may be omitted. Removing this term from the equation for A 2 in the systems (2) and (3) yields a single differential equation involving A 2 only.
-Linear model: In the absence of back flow from the peripheral compartment, the amount of compound in the central compartment is governed by the equation In this equation the input term k a A 1 has been replaced by the infusion rate q because, thanks to the large value of k a , within a very short time we have k a A 1 ðtÞ % q.
The right hand side of Eq. (13) has a unique zero, the plateau value A 2 , and it can be shown that Observe that i.e., the plateau value is smaller than the steady state value. Thus, the plateau value can be seen as the starting value of the second phase in which A 2 ðtÞ climbs further towards the final value A 2;ss .
Remark Because the system (2) is linear, the amounts A 1 ; A 2 and A 3 will depend linearly on the infusion rate q. This is indeed seen in the expression for the plateau value. Thus, -Nonlinear model: Without back-flow from the peripheral compartment, the dynamics in the central compartment is now governed by the equation and where A 2 is the unique positive zero of the right hand side of Eq. (17), or of the quadratic equation Therefore In Fig. 5 we show how in the nonlinear model, the plateau value A 2 and the plateau value normalised with respect to the infusion rate A 2 =q vary with q.
In contrast to the linear model, where this quotient is constant, in the nonlinear problem the normalised plateau is seen to be an increasing function of q, which connects two asymptotes. Expanding the expression for A 2 =q in (19) for small and large values of q, we find that The limits ' AE reflect the fact that (4) The small infusion limit in Eq. (21) demonstrates the sensitivity of the plateau value to changes in B max .

Conclusion
The simulations shown in Fig. 5, together with the analytical estimates derived from the model equations provide valuable Fig. 5 Variation of the plateau value A 2 ðqÞ (left) and the normalised plateau value A 2 ðqÞ=q (right) for the nonlinear model as they vary with q, when the data are k p ¼ 10 À4 h À1 , k 20 ¼ 10 À2 h À1 , the capacity takes the values: B max ¼ 10 4 (blue) 3 Â 10 4 (red) and 6 Â 10 4 (green) mg, and K M ¼ 100 mg diagnostic tools for identifying saturable elimination.
Increasing the infusion rate we observe (i) An increasing plateau value which, when normalised by the infusion rate q, is still increasing and is uniformly bounded above and below by positive limits ' AE . (ii) Simple explicit expressions for ' AE which yield quantitative information about k 20 and k p B max =K M . (iii) Additional estimates for B max , K M and k p can be obtained from the value of q at the transition from ' À to ' þ .

Terminal slope
In both models, the amount of compound A 2 ðtÞ in the central compartment converges, in the first phase towards the plateau value A 2 and then in the second phase towards the steady state A 2;ss . The rate of convergence towards these limits is characterised by the half-life (t 1=2 ) or the terminal slope k z . We obtain accurate approximations for the terminal slope for each of the models, which we denote by k ð1Þ z for the first phase and k ð2Þ z for the second phase, and discuss how k ð1Þ z and k where A 2 ðq; B max Þ is the plateau value. We deduce the following properties: (1) As we have seen in Fig. 5, the plateau value A 2 increases when the infusion rate q increases. Hence, it follows from (25) that k z ðq; B max Þ is a decreasing function of q.
(2) When q ! 1, then A 2 ðq; B max Þ ! 1 and hence, by (25), (3) When q ! 0, then A 2 ðq; B max Þ ! 0 and hence, by , (4) The terminal slope in the first phase k ð1Þ z ðq; B max Þ increases as B max increases. To see this note that according to Fig. 5, the plateau value A 2 ðq; B max Þ decreases when the capacity B max increases.
For the Second phase the terminal slope is well approximated by the formula The right hand side suggests the following properties: (1) k ð2Þ z ðq; B max Þ is an increasing function of q and a decreasing function of B max .
(2) By expanding the expression for k ð2Þ z ðq; B max Þ in (28) for small and large values of q we obtain Note that as q ! 0, the terminal slope k ð2Þ z ðq; B max Þ of the nonlinear model approaches that of the linear problem given by (23) with H ¼ B max =K M . Figure 6 illustrates and confirms the analytical properties presented above. For the linear model they will be proved in Appendix 1 and for the nonlinear model in Appendix 2.

Conclusion
The simulations displayed in Fig. 6, together with analytical expressions for the dependence on q of the terminal slope in the first and the second phase are a rich source of information for estimating the different parameters in the models. For both phases, the terminal slope depends monotonically-first down and then up-on q and tends to finite non-zero limits as q ! 0 and q ! 1 which can be computed explicitly.

Impact of slow leakage from the peripheral compartment
In many practical situations, data are only available for the first phase, and only predictions can be made about the second phase [6]. Clearly, during the long second phase, with its slow dynamics, the influence of leakage from the peripheral compartment may well be relevant. In light of the large capacity of the peripheral compartment this may result in significant losses.
In order to assess the impact of leakage, we modify the nonlinear model and increase the first order loss term in the equation for the peripheral compartment by a factor ð1 þ aÞ, where a [ 1. The equation for A 3 in the nonlinear system (3) then becomes whilst the equation for A 2 , which does not involve a, remains the same.
Because it is assumed that k p ( k 20 , the two-phase structure is not affected by moderate leakage. And because during the first phase the elimination term in the equation for A 3 is small and may be neglected, the first phase will hardly change when some leakage takes place from the peripheral compartment.
On the other hand, during the second phase the impact of leakage will be felt. For instance, leakage has an impact on the steady state values of A 2 and A 3 . They now become: where A 2;ss ðaÞ is the root of the quadratic equation Note that this equation is the same as Eq. (19) for the plateau value A 2 , except for the factor a=ð1 þ aÞ which multiplies B max . An elementary computation shows that Thus, when there is little leakage ða ( 1Þ, then A 2;ss ðaÞ is close to the steady state value A 2;ss given by (6) and when leakage is substantial ða ) 1Þ, the steady state value drops down to the plateau value A 2 given by (20). In Fig. 7 we show how the temporal behaviour of A 2 changes as the elimination from the peripheral compartment increases beyond the original back-flow into the central compartment. The rate of infusion is kept constant (q ¼ 5) and the elimination is increased from the original value (a ¼ 0) in four steps to a ¼ 0:5; 1; 2 and 4.
The simulations confirm the analysis: the two-phase structure remains intact, and in the first phase ( Fig. 7 left panel) the the additional elimination does not show up in Fig. 6 Terminal slopes: k ð1Þ z ðq; B max Þ (left) in the first phase and k ð2Þ z ðq; B max Þ (right) in the second phase versus the infusion rate q for the nonlinear model for two values of the capacity: B max ¼ 10 4 (red) and B max ¼ 3 Â 10 4 mg (blue) and the rate constants k a ¼ 10, k 20 ¼ 10 À2 , k p ¼ 10 À4 h À1 , and K M ¼ 10 2 mg the graphs. In the second phase (right panel) elimination does have an impact, and shows a drop in final steady state, starting from the original value (a ¼ 0) and approaching a value close to the plateau value when a ¼ 4. Evidently, the half-life in the second phase decreases as elimination from the peripheral compartment increases.

Conclusion
Elimination is a long-term phenomenon, as is to be expected since it takes place from the peripheral compartment which fills up slowly since k p is small. Nonetheless, the impact on the central compartment can be significant and, even for moderate elimination rates, can obliterate most of the growth beyond the first phase.

Small capacity and rapid distribution
To fully appreciate the effect of a large capacity of the peripheral compartment combined with a slow exchange between the two compartments, we conclude with a brief discussion of the dynamics of the nonlinear model for the converse situation: small capacity of the peripheral compartment combined with a fast exchange between the two compartments. Thus, we here assume that Since in this case the peripheral compartment has small capacity and direct elimination is relatively small, one expects that an iv bolus administration will lead to a large peak in concentration in the central compartment. In practical situations the height of this peak can be critical. Thus, to gain insight into this feature we focus here on dynamics after an iv bolus dose. As expected, soon after administration, A 2 ðtÞ jumps up to a high peak value A 2;max . This peak value is seen to increase rapidly with D 0 in a super-linear manner: when D 0 increases from 10 to 20 mg then A 2;max rises by about 4 mg and when D 0 increases from 60 to 70 the rise is about 7 mg, i.e., almost double the low-dose increase.
Thus, as in the large capacity case, the graphs of A 2 versus time exhibit a two-phase structure, albeit with a completely different shape. A brief initial phase, say for 0\t\t 0 , in which A 2 ðtÞ exhibits a violent up-and down swing which ends with A 2 at an intermediate plateau value A 2 , followed by a much longer elimination phase.
In order to analyse the dynamics of this system for the parameer values constrained by the conditions (34) and obtain an estimate for A 2;max it is necessary to transform the system to dimensionless variables. This analysis, carried out in Appendix 3, yields the following estimates for A 2;max .
Because the initial phase is short and the elimination rate k 20 is small, the total amount of drug in both compartments is conserved during this initial phase. i.e., Because of the larger value of k p the two compartments are quickly in quasi-steady state, so that after a brief initial adjustment, we may put Because Eqs. (36) and (37) both hold at t 0 , we may use Eq. (36) to eliminate A 3 from Eq. (37) to obtain from which we can compute the value of A 2 , right after the initial peak. For small and large dose D we find (cf. Appendix 3), which clearly demonstrates the super-linear behaviour of A 2 ðDÞ. For the terminal slope of the first phase k ð1Þ z ðB max Þ we find In order to determine the long time behaviour of A 2 ðtÞ, we add the equations for A 2 and A 3 from the system (3). This yields the equation because q ¼ 0. We now use the expression for A 3 given by Eq. (37), which is valid in the second phase to eliminate A 3 from Eq. (41) to obtain Using the expression for uðA 2 Þ this equation can be written as where u 0 ðA 2 Þ denotes the derivative of the function uðA 2 Þ. The terminal slope k ð2Þ z of the graph of A 2 as it approaches its steady-state value A 2;ss is given by Since u 0 ðA 2 Þ is a decreasing function of A 2 it follows that the terminal slope increases as A 2;ss increases, i.e, as D increases. In particular, since u 0 ðA 2 Þ ! B max =K M as A 2 ! 0, it follows that so that t 1=2 ! f1 þ ðB max =K M Þg lnð2Þ=k 20 % 38 h for the parameter values used in Fig. 4. We see that this estimate is confirmed in Fig. 4.

Conclusion
We find that for small capacity and rapid exchange between central and peripheral compartment the dynamics has a brief initial phase followed by a long terminal phase, with an appropriately defined plateau value in between. As in the previous case the terminal slopes yield sensitive markers that can be used to identify the impact of saturation on drug distribution. The plateau value informs about the capacity B max and K M , whilst the terminals slope yields estimates for B max , K M , and about k a when ðB max =K M Þ k p [ k a and about k p when ðB max =K M Þ k p \k a .

Discussion
We have compared the dynamics of two types of models for the distribution of a compound over a central and a peripheral compartment. In one type the elimination of compound from the central compartment into the peripheral compartment is linear, and the other it is saturable and hence nonlinear. In both models, the return flow from the peripheral compartment to the central compartment is linear. We have focussed on two contrasting extreme cases: (i) In one case, the capacity of the peripheral compartment is large and the back-flow is slow, and (ii) In the other case capacity and back-flow are respectively, small and fast. These cases can be viewed as bench marks in parameter space since they exhibit very different dynamics, each being endowed with its own characteristic ligand versus time graphs.
Both types of graphs exhibit a two-phase structure. However, within these two phases each case has its own characteristic behaviour: the large capacity peripheral compartment retaining ligand for a long time, whilst in the small capacity compartment the presence of ligand, though large, is short-lived.
It is demonstrated that saturable distribution can lead to disproportionately higher steady-state exposures. Specifically: • In the large capacity/slow distribution case, multiple ascending doses (MAD) yield disproportionately higher steady state exposures. • In the small capacity/fast distribution case, SAD yield disproportionately higher C max .
Thus saturable distribution models merit a careful analysis in light of the impact saturation may have on exposure.
In analysing these models subject to the conditions (46-49) listed in Methods, a mathematical framework has been created which can be used to analyse comparable models, which involve additional processes such as (i) leakage, or (ii) binding of the ligand to proteins, lipids and receptors in the central or the peripheral compartment, such as discussed in [9], or (iii) when the model involves additional compartments. This analytical machinery makes it possible to give quantitative estimates of the impact of these processes on the drug distribution between compartments and over time.
As an application of the methods developed in this paper, we show that leakage from the peripheral compartment may have considerable impact over a period of time. If this period extends beyond the period over which measurements are available the need for accurate quantitive predictions is evident.
Distribution over two compartments in which the peripheral compartment has a limited capacity, has much in common with tissue-binding. Here it is the tissue, viewed as a separate compartment, which can become saturated when maximal occupancy is reached. Thanks to this similarity in structure many of the results established in this paper can easily be transposed to the dynamics of tissue-binding.
The mathematical analysis is presented in a series of appendices. They offer an introduction to the use of such methods as (i) the use of dimensionless variables and parameters and (ii) multi-scale analysis. Dimensionless parameters are often a numerical measure of the relative importance of different processes involved in the model, such as direct elimination from the central compartment and distributional transfer between the compartments. Different time-scales are a common occurrence in pharmacokinetics and pharmacodynamics, often due to large differences in concentrations, in rate constants or in binding constants. They make it possible to simplify the often complex systems by means of singular perturbation theory (cf. [7,8]). The appendices demonstrate the practical usefulness of this theory for the study of complex pharmacokinetic and pharmacodynamic systems, and can serve as an introductory tutorial.
In summary, we have demonstrated a number of interesting dynamic properties of saturable distribution models which can be of value in practical modelling applications. In particular, we have shown that such models can account for disproportional accumulation evident in MAD data as well as disproportional increase in C max in SAD data. This is achieved by relaxing the assumption of linear distribution in the standard model at the cost of only one extra parameter per peripheral compartment. Saturable distribution models share many properties with models for tissue and receptor binding, which provides another attractive mechanistic underpinning for this class of models. For these reasons, we feel that the saturable distribution model deserves a more prominent place in the pharmacometrician's toolbox than it currently has. Here we try to promote this by providing a guide to its dynamics and, hence, applicability.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendixes Appendix 1: mathematical analysis of the linear problem
In Appendices 1 and 2 we focus on the large capacityslow distribution case (cf. 46 and 47). Thus, and initially A i ð0Þ ¼ 0 for i ¼ 1; 2; 3. In order to compare the ''weight'' of different terms in the system (2) we introduce dimensionless variables. Thus, we use the steady state values as reference value and put where A 1;ss , A 2;ss and A 3;ss are given by (6). Introducing these dimensionless variables we obtain the system dx dt ¼ k a ð1 À xÞ dy dt ¼ k 20 ðx À yÞ À H Â k p ðy À zÞ dz dt ¼ k p ðy À zÞ This is a linear system, so that on any bounded time interval the solution is bounded. It remains to introduce a dimensionless time variable. For this purpose we choose three different time scales, each corresponding to one of the rate constants in the system (47).

A very short time scale
We choose k À1 a h as a reference time and define the dimensionless time s 0 ¼ k a t. Introducing this time variable into (47) yields the system By assumption (9) the parameters m and e are small. For the parameter values of Table 1 for the linear problem which are used in Fig. 2, they are: e ¼ 10 À3 and m ¼ 10 À2 . It follows that given any finite time interval 0 s 0 \T 0 , then as m ! 0, the solution of the system (48) converges uniformly to that of the Reduced system: Remembering that initially, ðx; y; zÞ ¼ ð0; 0; 0Þ, it follows that for 0 s 0 T 0 , the solution of the original system (48) is well approximated by Evidently, xðs 0 Þ ! 1 as s 0 ! 1. The half-life s 0;1=2 -i.e., the time it takes for xðs 0 Þ to reach half of its final value, is equal to lnð2Þ. Thus, in the original time variable t, the half-life is given by t 1=2 ¼ lnð2Þ=k a ¼ 0:07 h, which amounts to about 4 min.

A short time scale
We choose 1=k 20 h as a reference time and define the dimensionless time s 1 ¼ k 20 t. Introducing this new time variable into the system (47) and putting x ¼ 1, we obtain the reduced system dy ds 1 ¼ ð1 À yÞ À e Â Hðy À zÞ dz ds 1 ¼ eðy À zÞ Using phase-plane arguments (cf. [11]) it can be shown that 0\yðs 1 Þ\1 for all s 1 [ 0. Therefore, 0\zðs 1 Þ\e for all s 1 [ 0 and hence, since e ( 1, we may put zðs 1 Þ ¼ 0, and approximate the first equation of (50) by where we have retained the product e Â H, because it may not be small. The steady state value in this first phase is y ¼ ð1 þ e Â HÞ À1 ; it is the plateau value y. It follows that and the half-life s 1;1=2 in this initial phase are given by In terms of the original variables we thus obtain When H ¼ 100 and k 20 ¼ 0:01, then A 2 ¼ 50 Â q mg and t 1=2 ¼ 34:7 h (cf. Fig. 2).

A large time scale
We choose k À1 p h as a reference time and define the dimensionless time s 2 ¼ k p t. Putting this time variable into the system (48) and setting x ¼ 1 then yields the system e dy ds 2 ¼ 1 À y À e Â H Á ðy À zÞ Again, by standard singular perturbation theory (cf. [7,8,10]) it follows that after a very short time 1 À y À e Â H Á ðy À zÞ ¼ 0 Using this expression for y in the second equation of the system (54), we obtain Remembering that zð0Þ ¼ 0, it follows that and the half-life s 2;1=2 in the second phase is given by and the terminal slope k z is given by Comparing the expressions (55) and (59) for the half-life in, respectively, the first and the second phase, we see that as the capacity increases, the half-life decreases in the first phase and increases in the second phase.

Appendix 2: mathematical analysis of the nonlinear model
Here we define the dimensionless variables in which for A 1 and A 2 we have chosen the same reference values, whilst for A 3 we have chosen the fixed capacity B max which serves as a uniform bound for A 3 (cf. (7)). In addition we define the dimensionless constants Substitution of these variables into the system (3) yields As with the linear model, within a very short time xðtÞ % 1 so that we may put xðtÞ ¼ 1 and reduce the system (63) to A short time scale As with the linear model we use the dimensionless time s 1 ¼ k 20 t. Introducing this variable into the system (63) and setting x ¼ 1, we obtain the reduced system Since e ( 1 it follows that for s 1 ¼ Oð1Þ, we may put zðs 1 Þ ¼ 0 and write the first equation in (65) as an equation for y only: Since f(y) is strictly decreasing, f ð0Þ ¼ 1 and f ð1Þ\0, it follows that f(y) has a unique zero y between 0 and 1. Plainly, y is one of the two roots of the quadratic equation One of these roots is negative and the other is positive. Plainly y is the positive one which is given by Because yð0Þ ¼ 0 and f ðyÞ [ 0 for 0 y\y, it follows from (66) that In order to determine the rate of convergence towards y we linearise equation (66) at y. Writing y ¼ y þ g, we can write equation (66) as in which Omitting the small rest term qðgÞ we conclude that the terminal slope is approximately given by Returning to the original variables, this translates into (cf. equation (25)): A graph of the terminal slope in the first phase, as it varies with q and B max , is shown in Fig. 6.
Example When q ¼ 5 we obtain for the parameter values given by Table 1, that y ¼ 0:558, and hence A 2 ¼ 279 mg, and f 0 ðA 2 Þ ¼ À0:0223 which is consistent with what we observe in the simulation shown in Fig. 3 for q ¼ 5.

A large time scale
We now use the dimensionless time s 2 ¼ k p t, which yields the system e dy ds 2 On the basis of singular perturbation theory (cf. [7,8,10]), we may put which yields an expression for z in terms of y: When we substitute this expression into the equation for z in Eq. (74) we obtain a single equation for the variable y only.
u 0 ðyÞ Â dy or, when we divide by u 0 ðyÞ, The terminal slope (in terms of s 2 ) is now given by In terms of the original time t this translates into We conclude from this expression that the terminal slope increases as the infusion rate increases, and that its limiting values values are given by k z ðqÞ ! k p as q ! 1 ð82Þ and A graph of the terminal slope in the second phase, as it varies with q and B max , is shown in Fig. 4.

Appendix 3: mathematical analysis: small capacity and rapid exchange
In Appendix 3 we focus on the small capacity -fast distribution case (cf. 47 and 48). Thus, We introduce dimensionless variables using K M as a reference value for the amounts and 1=k a as a reference time. Thus we put In terms of these new variables the system (3) becomes for q ¼ 0, dx ds ¼ Àx Note that for the parameter values of Table 2, used in Fig. 4, we obtain b ¼ 10, e ¼ 0:2 and l ¼ 0:04. In order to obtain a first estimate, we put l ¼ 0 and also e ¼ 0, except when e is multiplied by b since the product b e ¼ 2 is clearly not small. Solving the first equation we An elementary computation shows that When y ) 1 we scale Eq. (87) and write y ¼ dz. The resulting equation may then be approximated by Its solution which starts at the origin is given by Returning to y, we then deduce that y max ¼d 1 À be d À be ln d be ¼ d À be ln d be À be as d ! 1: In terms of the original variables, the expressions (93) and (93) yield the limits given in Eq. (39).