Impact of mathematical pharmacology on practice and theory: four case studies

Drug-discovery has become a complex discipline in which the amount of knowledge about human biology, physiology, and biochemistry have increased. In order to harness this complex body of knowledge mathematics can play a critical role, and has actually already been doing so. We demonstrate through four case studies, taken from previously published data and analyses, what we can gain from mathematical/analytical techniques when nonlinear concentration-time courses have to be transformed into their equilibrium concentration-response (target or complex) relationships and new structures of drug potency have to be deciphered; when pattern recognition needs to be carried out for an unconventional response-time dataset; when what-if? predictions beyond the observational concentration-time range need to be made; or when the behaviour of a semi-mechanistic model needs to be elucidated or challenged. These four examples are typical situations when standard approaches known to the general community of pharmacokineticists prove to be inadequate.


Introduction
In recent years application of mathematics in drug development has gained momentum. Even the FDA is considering approval of compounds in part on the basis of arguments based on modelling and simulation (cf. [1]). But there is a great variety of ways in which mathematics can play a role in drug discovery and development. On the one hand, the industrial scientist is often faced with the problem to make reliable predictions about such issues as optimal dose or assessment of safety, on the basis of data about onset, intensity and duration of response, when quantitative information about the underlying mechanism of action is limited. The challenge is then to combine available physiological knowledge, well-designed experiments and mathematical analysis to develop a model which can be used to make such reliable predictions. In addition, with expanding knowledge about biological and physiological processes, more systems-based studies are being carried out in which mathematical ideas about dynamical systems are used, for instance, to model complex regulatory networks, or gain insight in the behaviour of such networks, i.e. locate sensitive spots (cf. [2,3]).
We demonstrate the role mathematics can play in various aspects of pharmacology, such as (i) analysing complex data sets; (ii) using mathematical reasoning for dissecting model structure and acquiring quantitive information out of unconventional response-time courses; (iii) predicting the effect of chronic drug exposure on the basis of relatively short-time data sets in the context of disease & Lambertus A. Peletier peletier@math.leidenuniv.nl Large molecule compounds, such as monoclonal antibodies, exhibit interesting nontrivial interactions with their target, involving binding, saturation and target turnover. This results in complex ligand-concentration versus time courses. In Fig. 1 we show a typical such data set. Their dynamics is often referred to as Target-Mediated Drug Disposition and has been the subject of many studies (e.g. Mager et al [9][10][11], Gibiansky et al [12], Krippendorff et al [13], Peletier et al [14,15], Ma [16] and Dua et al [17]). Characteristically, the concentration-time courses display a series of phases. Initially they exhibit a quick drop, which may easily be missed if the first plasma sample is 12-24 h post dose. Then the curves display a concave bend towards a slower decline. Here, at higher exposure of ligand, one often has first-order linear (dose-proportional) kinetics.
The third typical phase is a convex bend downward with a shorter apparent half-life as we approach lower concentrations. This is where elimination of ligand through target internalisation starts to contribute. The kinetics is now nonlinear. Interestingly, the downward bend occurs at a ligand concentration which is independent of the dose. Throughout this phase the target route of elimination is more or less saturable, but less and less so as concentrations decrease.
Finally, the ligand concentration enters a slower terminal phase, again after a concave bend, with a longer apparent half-life.

The model
The disposition of the antibody (ligand, L) is described by a two-compartment disposition model, involving a plasmaand a tissue compartment, coupled to a target pool (R) with zero-order production and first-order loss. Ligand and target form a target-ligand complex RL via a second-order process. The complex can either be degraded into ligand and target via the first-order k off process or be irreversibly lost via k eðRLÞ (first-order internalisation). (Fig. 2) The combination of the second order formation and firstorder loss of complex results in the following nonlinear system involving the concentrations of ligand in the central compartment (L) and in the tissue compartment (L t ), and of target (R) and ligand-target complex (RL): with parameters defined by in which Input denotes the zeroth order input flux, Cl L the first-order non-specific clearance, Cl d the inter- compartmental distribution, and V c and V t the volume of the central-and the tissue compartment.

Steady states
In Fig. 3 we show how the steady-state concentrations of ligand, target and ligand-target complex vary as the infusion rate k infus changes when the parameter values are given by Table 1. Note that at steady state the ligand concentrations in the central and the peripheral compartment are the same.
For the values of Table 1, the dissociation constant is K d ¼ k off =k on ¼ 0:011 mg/L and the Michaelis-Menten constant K m ¼ ðk off þ k eðRLÞ Þ=k on ¼ 0:044 mg/L.
The graphs in Fig. 3 are complex and offer unique diagnostic material to asses the strength of the different processes and the values of the parameters under very general conditions. Besides, whether or not there is a tissue compartment makes no difference, the steady states are the same.
The curves for L ss , R ss and RL ss , have nontrivial shapes. For instance, in the figure on the left in which they are plotted linearly, they have the following properties: (i) They each consist of two segments that are approximately linear. (ii) The two approximately linear segments are joined at a narrow interface located at approximately k infus ¼ k syn ¼ 0:11 (mg/L)/h.

Mathematical analysis
The slopes of these segments are well-defined and can be computed explicitly. Hence, their dependence on the parameters of the system (1) is transparent. Specifically, for the ligand-receptor complex RL one can prove for the slope (A) at low infusion rates: At a critical value, when k infus % k syn , the growth of RL stops abruptly and the graph becomes flat. The level (R Ã ) at which this happens is given by: Thus, the complex increases more or less linearly up to a plateau where it abruptly levels off. The height of this plateau depends on two parameters only: the synthesis rate of receptor k syn and the elimination rate k eðRLÞ of ligandreceptor complex.
Similarly, for the ligand L ss versus k infus curve we find for small infusion rates: For the data of Table 1, K m ( R 0 so that the initial slope of L ss is much smaller than that of RL ss . This is also evident in Fig. 3. For large infusion rates it is possible to prove the following limit This shows that L ss climbs as k infus increases more or less along a straight line with slope 1=k eðLÞ and shifted to the right by an amount equal to k syn . Because L ss depends monotonically on the infusion rate, one can also express the steady state values of RL and R in terms of the steady state ligand concentration. The resulting formula's are in which According to (7) the graphs of R ss and RL ss have a familiar sigmoidal shape with different limits at large and small ligand concentrations. They are shown in Fig. 4: In particular, the limits at small and large ligand concentrations are given by Fig. 2 Schematic description of target-mediated drug disposition. Ligand L is distributed over a central-and a tissue compartment with respective volumes V c and V t , is eliminated again via a first order process (k eðLÞ ¼ Cl ðLÞ =V c Þ, and binds reversibly (k on =k off ) to the target R to form a ligand-target complex RL, which then is irreversibly removed via a first order rate process (k eðRLÞ ) J Pharmacokinet Pharmacodyn (2018) 45:3-21 5 and for both curves the ligand values for which they reach their maximum value is given by Detailed information about the asymptotic formulae (3)-(6) and the expressions (7) and (8)   Target suppression R ss and complex formation RL ss at equilibrium versus the steady state ligand concentration L ss . The fractional targetturnover rate k deg is faster (left), slower (middle) and equal (right) than the complex turnover rate k eðRLÞ concentrations depend on the infusion rate k infus . They are given in Gabrielsson and Peletier [5].

Conclusion
We have seen how the complexity of the TMDD model is elegantly depicted by plotting the steady states of the compounds L, R and RL versus the infusion rate k infus (cf. Fig. 3). The graphs can be computed explicitly and depend critically on the initial amount of target R 0 and what one could call a generalised dissociation constant which combines dissociation and internalisation of ligand-target complex: If K m ( R 0 , a common situation (cf [10,15,18]), the graphs are almost piece-wise linear and the explicit expressions for the slopes are quite simple. Dependence of target-ligand complex formation RL and receptor suppression R on the ligand concentration, proves to be described by graphs of simple sigmoidal functions (cf. Fig. 4). This introduces a potency parameter L 50 which is related to K m by the quotient of the two target elimination rates k deg and k eðRLÞ : Thus, a close analysis of the steady state properties of the three compounds involved in target-mediated disposition reveals a great deal about the relative importance of the sub-processes involved, and offers the possibility to acquire quantitative information about them.

Using visual inspection to estimate model parameters
This case study is aimed at demonstrating how mathematical reasoning can be used to help the modeller in choosing appropriate models on the basis of pharmacokinetic and pharmacodynamics data sets. This becomes especially important when little is known about the underlying pharmacokinetics, e.g., when the drug is not supplied to the plasma. Thus, for such situations the only data available are for response over time for different doses and forms of administration. Study of such problems is often referred to as dose-response-time analysis. It goes back to early papers by Levy [19] and [20] and Smolen [21]. For further references we refer to [6,[22][23][24][25][26][27][28]. Specific questions such as (i) How does mathematical reasoning address what is observed in onset, intensity and duration of response, and (ii) How to choose an appropriate model are discussed.

Data, model and equations
The data set records the locomotor activity, measured by the number of times moving rats interrupted three light beams in a cage when they were supplied by a drug, (dexamphetamine). In the absence of dexamphetamine the number of interruptions was negligible, but it goes up when the drug is given. The exact mode of action of drug is not known, and therefore an empirical mathematical model is proposed. Data shown in Fig. 5 were obtained and digitized from Van Rossum and Van Koppen [29]. They recorded the locomotor activity score after administration of dexamphetamine to rats at two dose levels (3.12 and 5.62 lg kg -1 ). The data set was unusual because (i) the rise and drop of response were approximately linear with time and (ii) the slope of the increase and of the decline in the locomotor activity score was independent of dose. In addition, the transition from increase to decline was relatively rapid.
Because the exposure to testcompound (drug) was not known, a classic biophase model was fitted, one in which the drug is administered through an intravenous bolus administration. The amount of drug A b in the biophase (in lg) is described by the equation: where D denotes the dose (in lg), and F Ã the biophase availability, i.e., the fraction of the dose that reaches the biophase, and k the elimination rate of the drug out of the biophase. The pharmacodynamic response R i.e., the number of interruptions per minute, is assumed to be described by a nonlinear turnover equation in which FðA b Þ is the drug-mechanism function through which the drug in the biophase impacts the response, k out;max the maximal elimination rate and K m the response at which the elimination rate is half-maximal. The combined biophase-and pharmacodynamic model is depicted in Fig. 6. Prior to administration of the compound no activity is observed, i.e., Rð0Þ ¼ 0. Two observations inform the selection of the function FðA b Þ: (i) The data exhibit a zero baseline. This means that Fð0Þ ¼ 0. (ii) As the drug dose increases, the initial slope of the data curves appears to reach a maximum.
They suggest a saturable drug-mechanism function which vanishes when the drug vanishes. It is of the following form: in which F max (response units Át À1 ), FD 50 (dose units) and n H correspond to the maximum drug-induced efficacy, the potency and the Hill-exponent.
The particular form of the turnover eq. (14) was selected in light of the approximately linear elimination of response which suggest saturation.

Mathematical analysis
In order to proceed from qualitative observations to quantitative estimates about the model, we employ the following observations: • Decline of the response: After the time T max of maximal response the graph has three conspicuous properties: (i) it is straight, (ii) it does not change with drug dose, and (iii) it exhibits a sharp angle as it approaches the baseline.
These characteristics of the response curve offer us an unusual insight into the dynamics of the model.
(a) At the time of maximal response T max the drug has been eliminated and FðA b Þ % 0, so that for t [ T max , the turnover equation is approximately reduced to When R ) K m , then eq. (16) simplifies to which shows that the slope of the graph is Àk out;max . On the basis of the data we obtain the estimate: (b) The sharp angle of the graph of R(t) as it approaches the baseline, i.e., when R % 0, can be accounted for by a small value of K m . (c) The response-time course associated with the higher dose peaks at about T max ¼ 0:8 h. If one assumes that approximately four half-lives have elapsed before the drug has been cleared from the biophase, this means that 4 Â t 1=2 % 0:8 h, i.e., the half-life of the drug in the biophase can be approximated by t 1=2 % 0:2 h; ð19Þ • Rise of the response For the higher dose the initial segment of the graph of R(t) is approximately a straight line. This suggest that during this period the stimulatory function is saturated, so that FðA b Þ % F max , and,  13) and (14). Drug is supplied to the biophase, eliminated through a first order process (k), and the amount of drug A b has a stimulating effect in the turnover equation for the response through a nonlinear function FðA b Þ. Loss of response is modeled by a saturable function with maximal (zeroth order) loss rate (k out;max ), half of which is reached when R ¼ K m provided R ) K m , the turnover equation is well approximated by Since the data show an initial slope (dR / dt) of 168 interruptions/minute/h, we deduce that F max % 168 þ k out;max interruptions/minute/h. Using the estimate for k out:max from the first observation, we conclude that Remark It is evident from eq. (14) that at steady state, the amount A b;ss should be small enough so that the production FðA b;ss Þ is smaller than the maximal rate of loss in order to reach a steady state. Thus a basic assumption in this model is that • Time to maximal response T max . To obtain a ball park value for the time to maximal response T max , we approximate the function FðA b ðtÞÞ, defined by (15), by a step-function. This choice is based on the fact that, as we argued before, the up-swing is more or less linear, so that the function FðA b ðtÞÞ appears saturated, and the decline is also linear and in addition, dose-independent, which suggests that FðA b ðtÞÞ % 0 after the peak-time T max . Thus, if A b ðtÞ is decreasing and crosses the level FD 50 , say at time T, i.e., when A b ðTÞ ¼ FD 50 , then we postulate that the function FðA b ðtÞÞ can be approximated by the step function Here HeavðsÞ denotes the Heaviside function which equals ?1 if s ! 0 and 0 if s\0. Thus, As the Hill-coefficient becomes larger, this approximation improves, i.e., This property follows readily from the classical limit lim n!1 The turnover eq. (14) can now be approximated by except for when R is small, specifically, when R ¼ OðK m Þ. Starting at baseline, the solution R Ã ðtÞ of this equation is given by: Because the biophase is assumed to follow intravenous bolus dynamics, as described by eq. (13), it follows that We assumed that F Ã ¼ 1, had taken the larger dose D ¼ 5:62 lg kg -1 , and estimated the elimination rate out of the system to be k ¼ 3:5 h À1 . By visual inspection, the time of maximal response was estimated by T max ¼ 0:8. Thus, we conclude from eq. (30) that the potency can be estimated as follows: Observe that for fixed A b the steady-state response R ss is formally given by Recall that according to the assumption (23), we have k out;max [ FðA b;ss Þ.

Conclusion
Using mathematical methods we have been able to estimate many of the parameters, such as k in the biophase model, and F max and FD 50 , as well as k out;max and K m in the pharmacodynamic model. These estimates can serve as preliminary estimates for further refinement by statistical software. By fitting the model (cf. equation (14)) to the experimental response-time data in Fig. 5, the final parameter estimates of Table 2 were obtained using the nonlinear regression software WinNonlin 5.2 (Certara Inc.). The six model parameters generally had high precision and were close to the analytically and graphically derived initial parameter estimates. The mathematical reasoning in the pattern recognition process thus proves to be useful when unconventional data (as in this case) need to be analysed.

Model predictions beyond the experimental range
Drug intake over long periods of time, such as is common in the treatment of chronic diseases, may harbour risks which are less evident over the period in which experimental data is available. Simulations on the basis of mathematical models, although predicated by the limited availability of data, may then yield indications of what kind of long time behaviour can be expected and how it can be influenced.

Data, model and equations
We consider an example of such a situation discussed by Peletier, de Winter and Vermeulen [7], in which data for  the plasma concentration of a compound was available for a period of 480 hours. They are shown in Fig. 7. This case-study involves a compound that is administered into a pool compartment and absorbed in a central compartment, from where it binds to two receptors through Michaelis-Menten type reactions (cf. Michaelis and Menten [30]) and dissociates according to first order kinetics. The amounts (in mg) in pool-and central compartment are denoted by, respectively, A 1 and A 2 . Binding to one receptor, which is probably located in the red blood cells is fast (amount A 3 mg) and binding to the other receptor, which is located in what is called the ''remote'' compartment, is slow (amount A 4 mg). The distributional model is illustrated in Fig. 8.
The model used to reach an optimal fit to the data shown in Fig. 7 is based on a 2-receptor model due to Snoeck et al [31] and is shown in Fig. 8. For comparison, the corresponding 1-receptor model is obtained from the above model by putting k RMT ¼ 0.
The 2-receptor model for the amounts of compound in the four compartments (A 1 ; . . .; A 4 ) translates into the following set of differential equations: where k RBC and k RMT are the distributional rate constants to the receptors in the red blood cells and the remote receptors, B max;1 and B max;2 the maximal capacity of these receptors and K d;1 and K d;2 the associated dissociation constants multiplied by the corresponding volumes. The infusion rate is D Á q mg/h, where q is the unit infusion rate, i.e. q ¼ 1 mg/h and D is the amount of compound that is supplied per hour. It is assumed that initially, there is no compound in any of the compartments or bound to the receptors, i.e., From t ¼ 0 onwards the compound is administered to the pool compartment through a constant-rate infusion of D Á q mg/h, which in this study is taken to be 40 mg/h. Two models are used to fit the data, which are given in aqua. The black curves are the individual fits made with the 2-receptor model (33) and the magenta curves are the individual fits made with the 1-receptor model obtained from the system (33) by putting k RMT ¼ 0. The grey curves are the population fits made with the 2-receptor model. In Table 3 we give the parameter values obtained by fitting the 2-receptor model to the data.
Consistent with these parameter values we assume throughout this section that the receptors in the red blood cells have a high affinity to the drug, a small capacity and a short half-time, all relative to the remote receptor. Thus, throughout we assume about the parameter values that: As drug flows into the system at a rate of D mg per hour, the system eventually settles on a steady state: (A 1;ss ; . . .; A 4;ss ). Plainly, the steady state amounts in the different compartments are given by  Note that A 1;ss and A 2;ss are independent of the binding constants and the capacities of the receptors, and increase linearly with the infusion rate D.
Remark As a preliminary observation we note that the first equation in the system (33) can be solved explicitly, resulting in the following expression for A 1 ðtÞ: Thus, A 1 ðtÞ ! A 1;ss as t ! 1 with a half-life of t 1=2 ¼ lnð2Þ=k a ¼ 0:28 h, or 17 min. This is exceedingly short for the time scale of interest. So effectively, it is permissible to put A 1 ðtÞ A 1;ss . Therefore, we shall be mainly interested in the amount of drug in the central or plasma compartment and in the two types of receptors.

Simulations
In Fig. 9 we show how the amount of compound in the plasma compartment, A 2 ðtÞ, evolves over time after the infusion has been switched on. Evidently, two phases can be distinguished: in the first phase, shown in the left panel, A 2 climbs to what appears to be a stationary state, which we usually refer to as the Plateau value. Subsequently, during a second phase, shown in the right panel, which extends over a much longer time, A 2 continues to climbs towards its final steady state, albeit at a much slower pace. In Fig. 10 we show simulations of the amount of compound in the two receptors: A 3 in the red blood cells and A 4 in remote tissue. The left panel shows that A 3 quickly reaches a constant value, which is close to its final steady state

Mathematical analysis
In order to understand the observations made about the simulations and answer such questions as (i) Which receptor governs the dynamics in the initial phase? (ii) How high does A 2 rise in the first phase, i.e., what would be a good estimate of the plateau value A 2 ? (iii) What would be the rate of convergence towards A 2 ? and (iv) What is the rate of convergence towards the final steady state in plasma A 2;ss , and analogous questions about the amount of drug in the two types of receptors.
To answer these questions, we need to compare the relative impact of the different terms in the system (33). Because the amounts in the four compartments A 1 ; . . .; A 4 vary widely, as, do the rate constants given in Table 4, it is necessary to transform to dimensionless variables and normalise the amounts with well-chosen reference values. In light of the steady state values of A 2 , A 3 and A 4 shown in (36), using D, B max;1 and B max;2 as reference values seems appropriate. Thus, we introduce the variables where we used 1=k RBC as a reference time.
Because e ( 1 the right-hand side of the third equation of the system (38) is very small. This means that zðsÞ % 0 for times s up to order 10 À1 Â e À1 % 10 3 . Thus, up till such time we may approximate the full system (38) by putting z ¼ 0 and replace it by the simpler system which involves only two equations and two unknowns: x and y: Note that while e is very small, when it is multiplied by B max;2 , as it is in eq. (41), the product is no longer small (cf. eq. (40)). This system (41) has a unique steady state ðx; yÞ, where x can be computed as the root of the function f(x) defined below: and y can then be computed from the right-hand side of the second equation in (41). Equation (42) is essentially a second order equation in x which can be solved explicitly. However, because of the many parameters the expression for x is quite messy. However, an important property of x can easily be inferred by inspection: Note that (i) the function f(x) is monotonically decreasing, (ii) f ð0Þ ¼ / and (iii) f ðx ss Þ\0, where x ss ¼ A 2;ss =D ¼ /=l. Therefore, 0\x\x ss and hence 0\A\A 2;ss , in agreement with the simulations shown in Fig. 9. Note that this conclusion is independent of the values of the parameters. For the parameter values of Table 4, the computed values for x and y are found to be x ¼ 57:4 and y ¼ 0:965 For small perturbations ðx þ nðsÞ; y þ gðsÞÞ, where nðsÞ and gðsÞ are small compared to x and y, one can derive a linear system of differential equations for nðsÞ and gðtÞ and compute the rate with which orbits converge towards the steady state ðn; gÞ ¼ ð0; 0Þ as s ! 1, and thus compute the terminal slope or terminal elimination rate k z , defined by  Table 4 this amounts to t 1=2 ¼ k À1 RBC Á s 1=2 ¼ 50 min (cf. [7]). Hence after 4 Â t 1=2 ¼ 200 h the plateau value has approximately been reached, consistent with the simulations shown in Fig. 9.
Beyond the first phase, A 2 and A 3 are in quasi-equilibrium and we may assume that This enables one to reduce the system (38) to a different, smaller, system making it possible to estimate the half-life of the convergence towards the final steady state A 2;ss . In fact, it is found that the terminal slope of this phase is k ¼ k RMT , so that the half-life is given by This is also consistent with the findings in Fig. 9 (Note that t % s). For details of the derivation of these estimates we refer to [7].

Conclusions
The mathematical analysis of the multiple receptor binding system demonstrates that care should be taken when using the model for making long-term predictions since such predictions may involve extended periods which well exceed the duration of experimental data. The final steady state of both binding processes may then be significantly higher than what is reached within the experimental time span. Therefore, long term exposure data will be needed to validate the model if used for future risk assessment. The insights obtained from this mathematical analysis will support the development of alternative models that exhibit the same short to medium term kinetics, but different long term kinetics provided chronic exposure data are available for model validation. For example, they could quantify the impact of small leakage, which over extended periods, may well be large (cf. [32]).

Vetting a model that yields counter-intuitive concentration-versus-time graphs
In pharmacology, mechanistic mathematical models are commonly developed on the basis of a combination of what is known about the underlying physiology and statistical methods which attempt to estimate the parameters in the model using experimental data. The resulting model is then employed to make predictions about optimal drug dose and generally, about the temporal behaviour of the drug and its effect. In general, before using the model in a clinical environment, it is ''challenged'' against different drug doses for which data sets exist, or for completely different data sets. In some cases this yields unexpected results. To get to the bottom of such an apparent anomaly a mathematical study of the model is then advised. In this case study we study a recent mechanistic model for Prolactin (PRL) which yielded unexpected results about the dynamics of prolactin (cf. [8] and [33]) when fitted to pharmacokinetic data from Kozielska et al [34].

The prolactin model
This case study involves a model designed to investigate the response of PRL to antipsychotic drugs, such as remoxipride or risperidone, in rats. The model, developed by Movin-Osswald et al., [35], which is based on the classical precursor-pool model [36][37][38] which distinguishes between PRL in plasma, and PRL in lactotrophs that serves as a precursor pool for the PRL in plasma. If P denotes the PRL concentration in the lactotrophs and R the PRL concentration in plasma, then this model is given by the following system of equations: Here k s denotes the zeroth order synthesis rate of PRL, k r the first order rate of release of PRL from the lactotrophs into plasma, and k el the first order elimination rate of PRL from plasma. The drug, at concentration C(t) in the brain, stimulates the release rate from lactotrophs through a standard saturable drug-mechanism function S(C) defined by where S max is the maximal extent of stimulation, SC 50 the drug dose for which stimulation reaches 50% of its maximal effect, and n H the Hill exponent. Evidently, the baseline ðP 0 ; R 0 Þ is given by Stevens et. al., [33] incorporated the fact that release of prolactin by the lactotrophs into plasma has a stimulating effect on the production of prolactin resulting in a positive feedback loop (see Fig. 11). See also Friberg et. al., [39].
They incorporated that effect into the model (46) through a multiplicative function of the synthesis rate k s that depends on the prolactin concentration in plasma, generalising the classical model (46) to the following model In this model f(R) is a non-decreasing function of R endowed with the following properties: (i) For R ! R 0 , the function f(R) is strictly increasing which vanishes when prolactin is at baseline, i.e., when R ¼ R 0 . (ii) To guard against sudden collapse, the feedback is switched off when the PRL concentration drops below the baseline value, i.e., f ðRÞ ¼ 0 when R\R 0 . Stevens et al [33] chose a function of the following form: f ðRÞ ¼ f 0 ðRÞ Á HðR À R 0 Þ; in which E max denotes the maximal stimulatory effect on the synthesis rate k s and EC 50 the increase of the PRL concentration above the baseline value at which half the maximal effect is achieved. In addition, HðR À R 0 Þ is the Heaviside function which vanishes for R R 0 and so ensures that f ðRÞ ¼ 0 for R\R 0 , and equals to 1 for R ! R 0 so that f ðRÞ ¼ f 0 ðRÞ for R ! R 0 .

Counterintuitive behaviour
When Proost and Taneja (Private communication) used PK data for risperidone obtained by Kozielska et al [34] (cf. Appendix) to drive the model above, they found the following counterintuitive behaviour (cf. Fig. 12): 1) For drug doses D ¼ 0:1 and 1.0 mg/kg they observed that simulations return R to the baseline R 0 as t ! 1. 2) For drug doses D ¼ 0:05 and 2.0 mg/kg they observed that simulations do not return R to the baseline R 0 as t ! 1, but converge towards a higher constant level, which we denote by R 1 . 3) The PK is fast, in that the half-life of drug in the central compartment is around 1 h, whilst the time for R to reach R 0 is about 10 h and to reach R 1 is around 40 h.
This suggests that, (i) when CðtÞ 0, there exists besides the baseline R 0 an additional steady state R 1 [ R 0 , and (ii) the drug dependence of the dynamics is not monotone and quite sensitive to drug dose. For practical situations this is very critical so that it is important to find the reasons for this behaviour.
In order to understand these phenomena it is necessary to study the mathematical properties of the model (49) with positive feedback function (50) more closely. That will be done in the next subsection.

Mathematical analysis
To simplify the equations and reduce the number of parameters, we introduce dimensionless variables. To this end we scale P and R by their respective baseline values and time by 1=k el and put: Introducing these variables into the system (49) and the feedback function f(R) we obtain where wðsÞ ¼ 1 þ SðCðtÞÞ and

Stationary solutions
Suppose that the drug concentration is constant, i.e. CðtÞ C ! 0, and wðsÞ w ¼ SðCÞ. Then, by (52) a stationary solution ðx; yÞ satisfies the pair of equations Substituting the second equation into the first yields the following equation for y: This is a quadratic equation in y which has the roots: Plainly, y 0 corresponds to R 0 , the baseline in the absence of positive feedback. However y 1 corresponds to a new stationary solution which is introduced by the positive feedback, denoted by R 1 .
It is illustrative to follow the dynamics of the system in what one may refer to as the state space, the (x, y)-plane in which the state, defined by x and y, travels. It is often called the Phase plane and the trajectory, traced by the concentration pair (x(t), y(t)), is called the Orbit. Plainly, at each point (x, y) in this plane the velocity vector q ¼ ðdx=ds; dy=dsÞ can be computed from the system (52). States at which x increases with time ðdx=ds [ 0Þ and where x decreases with time ðdx=ds\0Þ are separated by curves C x where dx=ds ¼ 0. Similarly, C y separates states where y increases, respectively decreases. The curves C x and C y are called the Null clines. Clearly, the stationary states are located at the points where C x and C y intersect.
When C ¼ 0, then by the system (52), the null clines are given by Fig. 13 Null clines in the phase plane spanned by the dimensionless concentrations of PRL in lactographs (x) and in plasma (y) for three ranges of b with respect to c (Here c ¼ 1 and b ¼ 0:2; 1; 3). The arrows indicate the direction of the velocity field q: vertical on C x and horizontalon C y C x : y ¼ x and C y : y ¼ 1 þ cðx À 1Þ b À ðx À 1Þ ; These null clines are shown in Fig. 13 for c ¼ 1 and b ¼ 0:2; 1:0 and 3. Notice that C y is fixed and C x moves to the right and up as b increases. As predicted by eq. (55), we see that the corresponding steady states are, besides ðx 0 ; y 0 Þ ¼ ð1; 1Þ: ðx 1 ; y 1 Þ ¼ ð0:2; 0:2Þ when b ¼ 0:2, and ðx 1 ; y 1 Þ ¼ ð3; 3Þ when b ¼ 3.
The null clines are very helpful in determining various aspects of the dynamics of the system, such as (i) the stability of the steady states, (ii) the large-time behaviour of orbits: e.g. where they go to and how they approach the stable steady states and (iii) identification of invariant regions, i.e., regions in the plane which trap orbits. Thus, the arrows in Fig. 13 suggest that when the positive feedback is small, i.e., when b ¼ E max is small, then ðx 0 ; y 0 Þ ¼ ð1; 1Þ is stable because all the arrows point towards it. However, when the positive feedback becomes stronger, and specifically, when b becomes larger than c, then ðx 0 ; y 0 Þ loses its stability and arrows point to ðx 1 ; y 1 Þ.
When the parameters in the system explicitly depend on time, the situation is more complex. In the system (52) the parameters are all constants except the coefficient wðsÞ which depends on s. However, since CðtÞ ! 0 very quickly (cf. Fig. 12, right panel), for most of the orbit we may put C ¼ 0, and hence w ¼ 1, after a brief initial period.
The specific parameter values employed for the system (52) by Stevens et al [33] are given in Table 4.
They yield for the dimensionless constants: a ¼ 0:10, b ¼ 3:47 and c ¼ 1:99. Thus, for the data used in [33] we conclude that b [ c, so that the right-hand graph in Fig. 13 applies. For the baseline we obtain R 0 ¼ 6:24 ng mL -1 and, using eq. (55), we obtain for the upper steady state R 1 ¼ R 0 Â y 1 ¼ R 0 ð1 þ b À cÞ ¼ 15:48 ng mL -1 , in agreement with the simulation shown in Fig. 12.
Summarising and rephrasing the observations made in the simulations shown in Fig. 12 we can state that RðtÞ ! R 0 as t ! 1 when D ¼ 0:1&1:0 mg=kg; RðtÞ ! R 1 as t ! 1 when D ¼ 0:05&2:0 mg=kg ð57Þ and the question is, why the PRL concentration does not go back to the baseline R 0 for any initial dose. In the next subsection we attempt to shed light on this observation.

Behaviour explained
In order to understand why the behaviour of the PRLconcentration is so sensitive to the drug dose, we turn to the phase plane and follow the orbit traced by the concentration pair (x, y) from its starting point ðx; yÞ ¼ ð1; 1Þ all the way towards its limiting state ðx 1 ; y 1 Þ. Specifically, we wish to know when the orbit tends to ðx 0 ; y 0 Þ and when to ðx 1 ; y 1 Þ, and how the drug dose D enters into this selection. For simplicity we first study the dynamics of the system (52) in the absence of a cut-off of the positive feedback. For this case we show in Fig. 14 the orbits in the phase plane for the drug doses D ¼ 0:05, 0.1, 1.0 and 2.0 mg/kg.
We observe in Fig. 14 that after a rapid introduction of risperidone, all orbits leave the baseline ðx 0 ; y 0 Þ ¼ ð1; 1Þ move up and to the left along an orbit which is initially tangent to the line After describing a big loop the orbits all return to a neighbourhood of the baseline ðx 0 ; y 0 Þ ¼ ð1; 1Þ from where they started. However, because b [ c the baseline is unstable and orbits move away from it, with the exception of two orbits: one from above and one from below, which tend towards ðx 0 ; y 0 Þ. Orbits which pass above these ''stable orbits'' (cf. D ¼ 0:05 and 2 mg/kg) ultimately converge towards the second equilibrium solution ðx 1 ; y 1 Þ. Those which pass below (cf. D ¼ 0:1 and 1 mg/kg) leave the first quadrant and they do this through the y-axis, below the line y ¼ 1. This dichotomy is clearly demonstrated in Fig. 14. Evidently, orbits entering the neighbourhood of ðx 0 ; y 0 Þ are very sensitive to small changes in paramater values or to the drug dose: a small change may flip them to the other side of the two stable orbits and change their subsequent course dramatically. This is the phenomenon that Proost and Taneja observed.

Including cut-off
Of course, by definition, both x and y cannot be negative and therefore realistic orbits should lie in the first quadrant. In order to ensure that x [ 0 and y [ 0, Stevens et al. cut off the positive feedback as soon as R\R 0 or y\1. This affects the null cline C x below the line y ¼ 1, i.e. in the region fðx; yÞ : y\1g, as shown in Figure 15. The cut-off makes C x vertical below y ¼ 1. Since the velocity vector on C x is also vertical, this section of C x is invariant and orbits cannot cross it.
Thus, the orbits are unaffected as long as they stay above the line y ¼ 1. However, when they dip below this line they change direction, move to the right and generally are captured in the triangular region between the two null clines C x and C y and the x-axis. In this region q points up and to the right and orbits can readily be shown to converge to the baseline ðx 0 ; y 0 Þ ¼ ð1; 1Þ as time tends to infinity. Thus xðsÞ % 1 and yðsÞ % 1 This is what we see demonstrated in the simulations of Fig. 16.
As we have shown, thanks to the cut-off of the positive feedback, we now observe two types of large time behaviour: (i) Orbits converging to the baseline ðx 0 ; y 0 Þ and (ii) Orbits converging towards the steady state created by the positive feedback. In both cases, orbits approach their limits from below and from the left. The fact that orbits pass through a neighbourhood of (1,1) makes for very sensitive dependence on the drug dose.

Conclusion
The mathematical analysis of the prolactin model with positive feedback has explained the unexpected behaviour of orbits and in particular their sensitivity to the drug dose. However, in the process it has done much more and given us great insight in the properties of the model. Thus, it has revealed the existence of two constant steady states, even in the absence of any drug-driven stimulation, i.e., two baselines of which one corresponds to the old baseline ðP 0 ; R 0 Þ. Only if b\c, i.e., if E max \ðEC 50 Á R 0 Þ, is the old baseline stable.

Overall conclusions
We have highlighted the impact of mathematical pharmacology in four case studies derived from previously published data and analyses ( [5][6][7] and [8]). In case study 1, we derive relations between quantities involved in TMDD systems at steady-state, such as target versus ligand and ligand-target complex versus ligand. These relationships give rise to a new expression of Potency, L 50 . This potency parameter is a conglomerate of binding affinity (k on & k off ), target turnover (k deg ), and ligand-target complex removal (k eðRLÞ ). Its applications will range from descriptions of in vitro and in vivo correlations to assessment of determinants of pharmacologically efficacious concentrations. The second case study involves locomotor activity and demonstrates how mathematical analysis, when combined with pattern recognition, can serve as a refined instrument for extracting qualitative as well as quantitive properties out of data sets, such as model structure, rate constants and dose dependence. Thus, Visual Inspection yields valuable parameter estimates which can be used in statistical data analysis for further refinement.
The third case study serves as a warning against the dangers of using multiple rate binding models for extrapolating beyond the experimental time range. This study emphasises the need for long-term data for making small, but in the long term significant corrections to exposure or environmental conditions. Deeper mathematical study may be required to uncover mysteries in the dynamics of PK-PD systems, such as large computational times or unexpected behaviour. The fourth case study focusses on the latter: a semi-mechanistic poolprecursor model for the dynamics of prolactin is found to exhibit-what appears to be-random dependence on drug dose. Small changes of the dose are seen to cause step changes in terminal behaviour. Thus, mathematical analysis exposes and makes explicit weaknesses of the model.
The take-home messages from this Communication are that mathematical and analytical techniques serve their purpose when we need to chisel out new structures (Case study #1), quantify patterns (Case study #2), make what-if? predictions (Case study #3) and diagnose models with hidden pathologies (Case study #4).
Acknowledgements The authors would like to thank Willem de Winter and Suruchi Bakshi for their stimulating collaborations about, respectively, Case studies #3 and #4 and the referees for their careful reviewing.
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.   Fig. 13: D = 0.05, 0.1, 1 and 2 mg/kg. The drug PK is simulated using the PK model for iv dose of risperidone [34] (cf. Appendix A)