Exact solutions and equi-dosing regimen regions for multi-dose pharmacokinetics models with transit compartments

Compartmental models which yield linear ordinary differential equations (ODEs) provide common tools for pharmacokinetics (PK) analysis, with exact solutions for drug levels or concentrations readily obtainable for low-dimensional compartment models. Exact solutions enable valuable insights and further analysis of these systems. Transit compartment models are a popular semi-mechanistic approach for generalising simple PK models to allow for delayed kinetics, but computing exact solutions for multi-dosing inputs to transit compartment systems leading to different final compartments is nontrivial. Here, we find exact solutions for drug levels as functions of time throughout a linear transit compartment cascade followed by an absorption compartment and a central blood compartment, for the general case of n transit compartments and M equi-bolus doses to the first compartment. We further show the utility of exact solutions to PK ODE models in finding constraints on equi-dosing regimen parameters imposed by a prescribed therapeutic range. This leads to the construction of equi-dosing regimen regions (EDRRs), providing new, novel visualisations which summarise the safe and effective dosing parameter space. EDRRs are computed for classical and transit compartment models with two- and three-dimensional parameter spaces, and are proposed as useful graphical tools for informing drug dosing regimen design.


Introduction
Mathematical models for the absorption, distribution and elimination of drugs are common in the pharmacokinetics (PK) literature. Typically, a drug's route through the body to its pharmacological effect site is modelled as a number of compartments, with transfer between compartments being governed by pharmacokinetic rate laws. It is common to consider only one or two compartments, with linear pharmacokinetics, resulting in low-dimensional linear ordinary differential equation (ODE) systems. However, such models are not sufficient to capture delay-type effects, whereby some time passes before the drug appears at measurable levels in the systemic circulation [44]. If a significant ''drug absorption delay'' [44,46] is observed, then a lag-time is sometimes introduced into solutions to the simple models to account for the delay, while avoiding any mechanistic considerations of the underlying delay roaaataaat esses. This simple approach may be used to paramaterise a system delay, but it is known that absorption delay is a complex process that is not switch-like. As such, lag-time models can give a poor characterisation of absorption-phase PK.
Transit compartment models have been proposed to capture delay effects in PK time courses, through a semimechanistic approach of increasing the number of compartments through which the drug is transferred en route to the central compartaament (blood) [26,27,32,33,44,46,47]. The development of ''full'' or accurate mechanistic physiologically-based PK models requires much experimental data and knowledge which may be unavailable. For systems exhibiting delays, transit compartment models therefore represent a physiologically plausible and mathematically practical alternative to the change-point approach of lag-time models.
While transit compartment models add additional complexity beyond one-or two-compartment models, their outputs in response to a single bolus dose to the first compartment may be derived analytically in certain cases. An analytical solution permits relatively straightforward approaches to both sensitivity analysis (particularly when varying the number of compartments) and model fitting. In [44], the response to a bolus dose is considered for an ntransit compartment model with an additional ''absorption compartment'' between transit compartment n and the central circulation. An analytical solution for the drug concentration in the nth transit compartment is presented, and used as an input to the absorption compartment ODE, together with the Stirling approximation, to transform a discrete optimisation problem to a continuous one for the purpose of data fitting. The analysis and parameter estimation is limited to the case of single bolus dose, and exact solutions are not found for the absorption and central compartments. Further mathematical properties for transit compartment models in pharmacodynamics are presented in [55].
Drug dosing regimens often use multi-dosing treatments, whereby a regular dose is given at regular specified dosing intervals. For intravenous (IV) or oral administration of a drug, the analysis of a one-or two-compartment model with periodic bolus input yields analytical solutions for the drug level in the central compartment (e.g. [11,24,43,45]). Time courses of drug level (or concentration) simulations show transient and steady-state (periodic) profiles which may then be compared with minimum effective and maximum safe levels which define a therapeutic window. We will introduce the idea of an equidosing regimen region of (dose,interval)-space, which gives a guide for selecting therapeutic equi-dosing regimens.
Given that it is accepted that drug absorption delay may be a significant pharmacokinetic effect, analysis of transit compartment models, incorporating multi-dose inputs, appears to be a valuable pursuit. Some attention has been paid to this problem in the PK literature. Shen et al. [46] extend the work of [44] to derive a solution to the multidose problem at the nth transit compartment using the method of superposition. However, the challenge remains to solve for the drug level in the central compartment exactly, and to use this result as a platform for further analysis, including design of safe dosing regimens.
In this paper, we present new mathematical and graphical results which both generalise transit compartment PK models and summarise dosing regimen constraints given by therapeutic ranges imposed on these models. In 'Methods: multi-dosing models with and without transit compartments-formulation', we formulate linear ordinary differential equation (ODE) models for one-compartment, twocompartment and transit compartment pharmacokinetics, extending the work of [44] to consider drug level in the central compartment. In 'Analytical solutions for equidosing regimens', we review standard analytical results for multi-bolus and multi-infusion dosing, and derive a new exact solution for the transit compartment model with general number of compartments and doses. This new generalised solution and its improvement over existing models comprise our first main contribution. In 'Results: equi-dosing time courses', we present simulations and data fitting using the new analytical solutions, illustrating their predictive capability, and highlighting the error between the new exact solution and Stirling approximation solution of [44] for a single bolus dose. In 'Results: equi-dosing regimen regions', we present our second main contribution, namely the new concept of equi-dosing regimen regions (EDRRs), which provide a novel visualisation to summarise constraints on dosing regimen parameters. We conclude in 'Discussion' with a discussion of our main results, highlighting our contributions to the PK and mathematical modelling literature.

General compartmental model schematic
We use a compartmental approach to model a drug's route from administration to the systemic circulation. Ultimately, the systemic circulation is treated as the final compartment in a cascade, hereafter referred to as the central compartment. The central compartment drug concentration (the drug amount per volume of distribution) is responsible for responses at drug effect sites [45], and we consider the drug level a c as the output in each of our models.
For intranvenous (IV) dosing, the drug immediately appears in the central compartment upon administration. We will refer to the corresponding model as a singlecompartment or one-compartment model (Fig. 1, model (M1)). The amount of drug in the central compartment (the ''drug level'') in this model is governed by an ordinary differential equation (ODE) which describes linear pharmacokinetics, whereby the drug is eliminated from the compartment as a first order process with elimination rate constant k e . A two-compartment model (Fig. 1, model (M2)) in which drug appears in the central compartment via an absorption compartment is often used to model oral dosing [43,45], where the absorption compartment is representative of the gastrointestinal (GI) tract. In response to a bolus dose input, the central compartment drug level a c in (M2) is both delayed and smoothed in comparison with the absorption compartment level a b . In order to model a more pronounced delay by way of a semi-mechanistic compartmental schematic, we consider a transit compartment cascade feeding the absorption compartment, as in [44] (Fig. 1, model (Mt)). We note that such a modelling approach corresponds to the so-called linear chain trick [23].

Dosing regimen inputs
Dosing patterns in therapeutics often consist of multidosing treatments, whereby doses are administered periodically [11,43,45,46]. Models for multi-dosing are typically analysed under the assumption of equi-dosing, whereby both the dose D 0 and dosing interval (time between doses) T are constants. In this case, the couple ðT; D 0 Þ constitutes the dosing regimen. Here we principally investigate equi-dosing regimens in which the input is a fixed bolus dose administered periodically to the central, absorption or first transit compartment (see Fig. 2, regimen (Beq)). We also consider a simple perturbation to this regimen, where a loading dose D L (greater than D 0 ) is administered at t ¼ 0, followed by equi-dosing (Fig. 2, regimen (BeqL)). This regimen is common in practice, such that a loading dose helps to achieve therapeutic drug levels rapidly, while the subsequent equi-dosing maintains therapeutic levels [45].
We further consider the case of equi-infusion dosing, whereby for model (M1), the input is periodic constant infusions of drug to the central compartment, over fixed ''on'' time intervals, separated by fixed ''off'' intervals ( Fig. 2, regimen (Ieq)).

Model assumptions and considerations
Transit compartment models (TCMs) for PK typically take the drug amounts in each compartment (a i for i ¼ 1; . . .; n, a b and a c in Fig. 1) as state variables (see, e.g., [26,44,46]). The bioavailability factor F, which is the fraction of drug dose ultimately absorbed into the systemic circulation, is an important consideration. Existing TCM models and simpler models introduce this ''correction'' factor at different points in the cascade [23,26,34,43,44,46]. Here we follow [23,43] in introducing the factor immediately, such that the first compartment in the cascade is fed by the effective dose (F Â dose).
We consider equi-dosing regimens for both IV infusion and n-transit compartment cascades, to allow analysis of drug level time course features in general. Together with prescribed therapeutic ranges, our analysis will indicate regions of equi-dosing regimen parameter space which give safe and effective treatments at steady-state. A therapeutic range is typically defined by minimum effective and maximum safe drug concentrations in the central compartment. For a given drug amount a c , the corresponding drug concentration C c is given by a c =V, where V is the calculated volume of the central compartment [25,43,45]. Therefore, for a fixed, known volume V, we can state a corresponding therapeutic range in terms of drug amounts, requiring D me \a c \D MS ; ð2:1Þ where D me and D MS are the minimum drug level for therapeutic effect and maximum safe drug level respectively. Transfer from absorption to central compartment is a first order process, with rate constant k a (we consider k a [ k e [45])-and k 6 ¼ k e ; k a [44]. (Mt) Transit-compartment model-the input dose first appears in the first of n transit compartments, which contain the amounts a 1 , a 2 ; . . .; a n . Drug is transferred through the cascade of n transit compartments, as first order process with rate constant k for each compartment. From transit compartment n, drug is transferred to the absorption compartment. For all three models, the (first order) elimination rate constant is k e For the transit compartment model, the transit cascade prior to the absoprtion compartment consists of n compartments, each with an elimination rate k. The mean transit time (MTT) for this cascade is given by (see [44]) ð2:2Þ

Ordinary differential equation formulation
For a multi-dosing regimen, the problem can be stated as an initial value problem (IVP) for the state variables (e.g. a i for i ¼ 1; . . .; n, a b and a c for model (Mt)) with impulsive drug inputs to the first compartment described by a dosingrate forcing function comprising Dirac delta functions for discrete impulses (as in [26,34]). Firstly, the IVPs we consider for single-compartment IV multi-dosing and two-compartment oral multi-dosing are summarised in Table 1. For equi-bolus dosing with M doses D 0 at time intervals T starting at t ¼ 0 (regimen (Beq)), the input rate to the first compartment is FD 0 dðt À ðj À 1ÞTÞ: ð2:3Þ If the first bolus dose is replaced with a larger loading dose D L (regimen (BeqL)), then the forcing rate function is FD 0 dðt À ðj À 1ÞTÞ: ð2:4Þ For equi-infusion dosing (regimen Ieq) with infusion rate k in , infusion ''on'' duration T, infusion ''off'' duration t f , and M infusions, the forcing rate function is where H is the Heaviside function. We will consider solutions to the IVPs given in Table 1 in constructing the associated equi-dosing regimen regions. Beyond these relatively simple models, our analysis extends to the n-transit compartment model with input into the first transit compartment (model (Mt)). The governing with drug infused into central compartment, starting at time t ¼ 0, at a rate k in (mg h À1 ). (Ieq) Equi-infusion dosing, with drug infused into central compartment at a rate k in , periodically with period T. Each dosing interval infusion ''on'' for duration t f , then infusion ''off'' for duration ðT À t f Þ equations consist of ðn þ 2Þ ODEs, which may be written in matrix form as

Analytical solutions for equi-dosing regimens
Here we present exact solutions for central compartment drug levels under a variety of equi-dosing regimens. Single-compartment and two-compartment model solutions are included for comparison with transit compartment model (TCM) solutions, and to aid the construction of equidosing regimen regions in 'Results: equi-dosing regimen regions'. The exact TCM solutions represent an improvement on existing approximate solutions [44].
Exact solutions for one-compartment and twocompartment models with equi-dosing In Table 2, we list exact solutions for drug level in the central compartment a c ðtÞ for one-compartment and two-compartment model IVPs formulated in Table 1, under equi-dosing inputs given by Fig. 2 and Eqs. (2.3)-(2.5). These include well-known solutions (e.g., [6,19,45]); for comparison with the TCM solution, their derivations may be found in detail in Appendix 1. Throughout, H is the Heaviside function, and ð3:14Þ Solutions may be written compactly without summation notation by considering, for example, the central compartment drug level after the Mth dose, a M c ðt M Þ. The steadystate (T-periodic) drug level function is denoted a 1 c ðt 1 Þ, where t 1 is the time since the start of the dosing interval. a c ð0Þ ¼ 0:

M2 Beq
Exact solution for transit compartment model with equi-bolus dosing More generally, the transit compartment model (Mt) comprises a multi-compartment oral absorption process, and the ODEs may be written in matrix form as in (2.6). We consider equi-bolus dosing, i.e. forcing input (Beq), so that g 1 ðtÞ ¼ g B ðtÞ. The exact solution may be written using the matrix exponential [9,27], or by using Laplace Transforms (Appendix 1.4). For the transit compartments, we find that Hðt j Þt iÀ1 j e Àkt j ; i ¼ 1; . . .; n; ð3:15Þ where t j ¼ t À ðj À 1ÞT. We note that this result is equivalent to the multi-dose analytical result of Shen et al [46] which was derived via superposition arguments, to be used as an input to their central compartment module, and also to the Savic single-dose solution [44] if M ¼ 1. We now extend our calculations to establish solutions for the absorption and central compartment drug levels, which in effect gives a general, analytical multi-dose solution to the Savic problem [44]. The solution for the absorption compartment (see Appendix 1.4) is Hðt j Þe Àk a t j c À n; ðk À k a Þt j Á ;

Steady-state behaviour
The derivation of the steady-state solution is more involved than for the earlier models (see Appendix 1.4.1). For the transit compartments, we find that The coefficients a i ð0Þ (the steady-state dosing interval initial values) may be found using ð3:20Þ together with the recurrence relation (for i ¼ 2; . . .; n) where S is the Stirling number of the second kind [40], given by ðÀ1Þ p q p ðq À pÞ n ; where q p ¼ q! p!ðq À pÞ! is the binomial coefficient; ð3:24Þ and taking Sð0; 0Þ ¼ 1. For the absorption compartment, we find cðn À p þ 1; ðk À k a ÞTÞ; ð3:26Þ and b a ¼ e Àk a T 1 À e Àk a T : ð3:27Þ Finally, and ultimately, the steady-state profile in the central compartment is given by e Àk e t 1 " cðn À p þ 1; ðk À k e Þt 1 Þ À e Àk a t 1 cðn À p þ 1; ðk À k a Þt 1 Þ

#)
; ð3:29Þ and b c ¼ k a ðk a À k e Þð1 À e Àk e T Þ : ð3:30Þ We now have in place new analytical solutions for a general M-equi-dose input to an n-transit compartment model with absorption and central compartments. These solutions may be used to predict drug level dynamics exactly, rather than approximately (see (4.5)). Further, the steady-state solutions may be used to guide safe and effective dosing regimen design, given a specified therapeutic range.

Computational evaluation of the lower incomplete gamma function
In order to use the analytical results of the previous subsection, a computational method for evaluating the lower incomplete gamma function is required. The definition itself (3.17) immediately suggests several approaches for a given n, including numerical evaluation of the integral for given t, computation of the truncated exponential sum, and using symbolic computation to derive an exact expression for the integral, then evaluating at given t. In fact, more efficient methods for evaluating this function have received attention in the mathematical literature, with many involving series and continued fraction expansions [1,2,17,41,50]. In MATLAB, the built-in function gammainc may be used [35], and is our preferred evaluation method due to its accuracy and run time (see Appendix 1.4.2). In software where such a function is not avaialble, the following relationship between the lower incomplete gamma function c and the cumulative distribution function F C for the gamma distribution may be used ( [35]): cðn; tÞ ¼ CðnÞF C ðt; n; 1Þ: ð3:31Þ Here, n is taken as the shape parameter of the distribution, and the scale parameter is unity. Furthermore, the loggamma function is a built-in function in many software packages, and the exponentiated log-gamma function is often used in situations where numerical difficulties may arise in evaluating the gamma function directly [35].

ð3:32Þ
The log-gamma function is available in PK analysis packages and languages including NONMEM [4], MLXTRAN/MONOLIX [31], PharmML [48], and also in MATLAB [35] and Excel [37]. Each of these packages also has the exponential function and cumulative gamma distribution F available.

Results: equi-dosing time courses
Here, we present simulated time courses of drug levels, using the analytical solutions given in 'Analytical solutions for equi-dosing regimens'. In all cases, we have computed using MATLAB [35].

IV equi-bolus dosing-one-compartment model
In Fig. 3a, we plot a typical drug level time course for the IV equi-bolus dosing problem (M1,Beq), which has solution given by (3.2), and steady-state profile given by (3.3). The characteristic features of the time courses include jump discontinuities at t ¼ jT, exponential decay over each dosing interval, and approach to a periodic steady-state [45]. Clearly, acceptable dosing regimens would only exist for certain ðT; D 0 Þ choices. Also shown is the continuous infusion profile given by a c ðtÞ ¼ Fk in k e ð1 À e Àk e t Þ; ð4:1Þ 16)). The dosing interval average drug level approaches the corresponding infusion steadystate level, which is apparent from the graph, and from .3) and comparing with the steady-state of (4.1). Finally, we note that the administration of a loading dose at time t ¼ 0 may give a treatment that is immediately and always therapeutic, and drug levels close to the steady-state.

IV equi-infusion dosing
In Fig. 3b, we plot a typical drug level time course for the IV equi-infusion dosing problem (M1,Ieq), which has solution given by (3.6), and steady-state profile given by (3.8). The characteristic features of the time courses include derivative discontinuities (but continuous drug levels) at t ¼ jT and t ¼ jT þ t f , exponentially decaying rise followed by exponential decay over each dosing interval, and approach to a periodic steady-state [45]. Again, acceptable dosing regimens would only exist for certain ðT; D 0 Þ choices.

Single continuous infusion as limit of IV equibolus dosing
In Fig. 3c, we plot a drug level time course for IV continuous infusion (4.1) for infusion rate k in , together with an equi-bolus dosing (M1,Beq) solution (3.2) for which the dosing rate is D 0 T ¼ k in with very short dosing interval T. It is apparent, and intuitively known, that continuous infusion represents a limit of a corresponding equi-bolus regimen for high dosing frequency. In Appendix 2.1 we offer a proof of this result using l'Hopital's rule.

Oral bolus equi-dosing-two-compartment model
In Fig. 3d, we plot a typical drug level time course for the oral equi-bolus dosing problem (M2,Beq), which has solution given by (3.10), and steady-state profile given by (3.11). The characteristic features of the time courses include derivative discontinuities at t ¼ jT, a two-phase profile (absorption then elimination) over each dosing interval, and approach to a periodic steady-state [45]. We note that the administration of a loading dose at time t ¼ 0 may give a treatment that reaches therapeutic level earlier, but for which there is still some nonzero waiting time before the therapeutic level is reached. It is clearly possible to give a loading dose which ensures that therapeutic drug level is both reached within the first dosing interval and is maintained thereafter.
The a c ðtÞ time course approaches a T-periodic steadystate profile, which will be non-monotonic for all regimens (even those for which the time course is monotonic for the Fig. 3 Drug level time courses. Throughout, we take F ¼ 1. Where shown, D me and D MS represent hypothetical minimum effective and maximum safe levels respectively, giving the therapeutic range ½D me ; D MS . Where shown, the steady-state profile overlays the final dosing interval for comparison. a IV equi-bolus dosing (M1, Beq), with and without loading dose. D 0 =500 mg, T = 12 h, k e ¼ 0:0692 h À1 (taken from [10]). Loading dose (M1,BeqL) has D L ¼ 800 mg. Continuous infusion at a rate k in ¼ D0 Oral equi-bolus dosing (M2,Beq) with D 0 = 500 mg, T=12h, k e ¼ 0:0692 h À1 , k a ¼ 0:7 h À1 . Loading dose (M2,BeqL) has D L ¼ 800 mg early dosing intervals), exhibiting both absorption and elimination phases. Expressions for the peak drug level and the peak timing are given in Appendix 2.2. Maximum and minimum drug levels for a number of models will be used in 'Results: equi-dosing regimen regions' to construct equidosing regimen regions, which give a summary guide for regimen design.
Transit compartments-smoothed delays, lag time and data fitting (single-dose) In Fig. 4, we demonstrate the delaying effect of a train of transit compartments for a single dose regimen. It is clear (Fig. 4a) that a delta-function-like bolus dose input to the first transit compartment effects a ''spread-out bolus'' dose in later transit compartments. Eventually a spread-out bolus profile is seen for the final transit compartment, which becomes the input to the absorption compartment, centered around t ¼ MTT ¼ n k . As in [44], we consider the transit compartment cascade delaying the appearance of bolus dose in the absorption compartment of a standard two-compartment oral dosing model. We see that (Fig. 4b,  c), for a fixed time lag t lag ¼ MTT taking k ¼ MTT n , the pure delay (time-lag) profiles (see [45]) ð4:2Þ and a lag c ðtÞ ¼ 0 0 t t lag k a k e À k a FD 0 e Àk a ðtÀt lag Þ À e Àk e ðtÀt lag Þ t [t lag 8 < : ; ð4:3Þ are approached by the equivalent transit compartment approximations for increasing n. Such ''smoothed delay'' profiles may well capture experimental data better than nodelay or pure-delay models [44]. Indeed, we find a better fit to published data for a single dose of the drug glibenclamide [44] using a transit compartment model than using a pure time-lag model (Fig. 5). For the least-squares data fitting shown in Fig. 5, we use the optimisation function fminsearch in MATLAB [35], with the objective function being the sum of squares between data and simulation at the data points. For each fixed n in turn, and for the time-lag model (for which t lag is one of the fitted parameters), the optimisation routine is run for 1000 iterations.

Parameter identifiability
We note that applying optimisation routines to estimate PK parameters for n ¼ 1 (a single transit compartment) should be with caution, since in this case, k e is identifiable but k and k a are unidentifiable. This is readily seen by considering a single bolus dose to the transit compartment for cases (i) k ¼ k 1 and k a ¼ k 2 , and (ii) k ¼ k 2 and k a ¼ k 1 , for rate constants k 1 and k 2 . In both cases, the inflow rate to  ,B1). Throughout, we take F ¼ 1, D 0 = 500 mg, k e =0.0692 h À1 k a =0.7 h À1 (for hypothetical drug described in [10]). Here, t lag ¼ MTT ¼ 3h, and k ¼ MTT n . a Drug level a i ðtÞ for compartments i ¼ 3; 20; 60; 100 of a cascade with n ¼ 100 transit compartments. b Absorption compartment level a b ðtÞ for cascades with n ¼ 3; 20; 60; 100 transit compartments, with equivalent time-lag profile for pure delay to absorption compartment. c Central compartment level a c ðtÞ for cascades with n ¼ 3; 20; 60; 100 transit compartments, with equivalent time-lag profile for pure delay to absorption compartment the central compartment k a a b ðtÞ is found from (3.15) to (3.16) to be The central compartment drug level will be identical for both cases, hence k a and k are not uniquely identifiable from the single output a c ðtÞ. This is a manifestation of the so-called flip-flop phenomenon for two-compartment kinetics [45]. The parameter k e is identifiable. An example computation is shown in Appendix 2.3. For n [ 1 this phenomenon is avoided, and all parameters are uniquely identifiable (see Appendix 2.3).

Transit compartments-exact versus approximate solutions
The original transit compartment schematic presented in [44] has been key to our analysis. A significant advance in our work is the development of an analytical solution which solves the problem exactly. The original single bolus dose analysis [44] uses exact solutions for each transit compartment, but employs the Stirling approximation for the factorial, to solve the following system for the final two compartments.
t nÀ1 e Àkt ; ð4:5aÞ ð4:5bÞ A comparison between a typical numerical solution of the approximate model (4.5) and corresponding exact solutions given by (3.16)-(3.18) for a single bolus dose is given in Fig. 6. It is clear that the new exact solutions provide a significant improvement in accuracy over the approximate solutions, particularly for n 4, for which the relative error in the peak a c value is between 3% and 8.4%.

Transit compartments-equi-dosing
While the published data and modelling in [44] focus on a single dose, we naturally wish to use such models for simulating time courses under multi-dosing regimens. The local maxima and minima in a multi-dosing timecourse prediction as the system approaches a periodic steady-state must be considered in regimen design [12]. In Fig. 7, we show predicted equi-dosing drug levels together with data used for fitting to a single dose, as in [12]. The approach to a periodic steady-state is clear in each case. For each timecourse, we note that the transit compartment profile is , with fitting to experimental data for drug glibenclamide, taken from [44], using WebPlotDigitizer [42], in response to a 3.5 mg dose. Original data converted from concentration to drug level using volume of distribution of 3.79l [44].

Results: equi-dosing regimen regions
It is common to consider potential therapeutic protocols which will give steady-state drug levels within a prescribed therapeutic range by simulation with multiple dosing regimens and observing whether the steady-state falls within that therapeutic range [12,25,43,45]. One can approach dosing regimen design iteratively, simulating in this manner and adjusting base parameters until a theoretically safe and therapeutic regimen is found [25]. Here we present a novel and alternative analysis which will capture the constraints imposed on the dosing regimen parameters by the therapeutic range, which is given by D me ¼ minimum effective level and D MS ¼ maximum safe level: ð5:1Þ We propose that equi-dosing regimen regions (EDRRs), which are regions of the parameter space giving acceptable regimens, may be used to summarise steady-state constraints and guide regimen design from the outset of any investigation.

Two-parameter dosing regimens
Equi-dosing regimen region for IV equi-bolus dosing Now a 1 c ð0Þ ¼ FD 0 1Àe ÀkeT and a 1 c ðT À Þ ¼ FD 0 e ÀkeT 1Àe ÀkeT , with a 1 c ðt 1 Þ decreasing, so we require that D me \ FD 0 e Àk e T 1 À e Àk e T and FD 0 1 À e Àk e T \D MS ; so that the region of ðT; D 0 Þ-space for acceptable dosing regimens is that corresponding to the inequality constraints The upper bounds on T and D 0 for the safe and effective dosing region are given by We see in Fig. 8a that the acceptable EDRR for equi-bolus dosing is given by a petal-shaped region. The two curves divide the ðT; D 0 Þ parameter space into four regions; the three regions other than the EDRR correspond to steadystate drug levels which are unsafe, sub-therapeutic, or both unsafe and sub-therapeutic over subintervals of their periodic timecourses. Illustrative drug level time courses are shown in Fig. 9.
Clearly the EDRR could be used at the outset of any investigation to guide regimen design, prior to any time course simulation or experiment.

Equi-dosing regimen region for oral equi-bolus dosing
For model (M2) with forcing (Beq), we consider the steady-state solution (3.11). It is straightforward to show (Appendix 2.2) that the minimum and maximum levels a 1 c;min , a 1 c;max , and peak time t Ã 1 are given by ; ð5:5aÞ 1 À e Àk e T 1 À e Àk a T ; ð5:5bÞ ka Àke : ð5:5cÞ Again requiring that D me \a 1 c ðt 1 Þ\D MS , we find that the region of ðT; D 0 Þ-space for acceptable dosing regimens is that corresponding to the inequality constraints    Fig. 11 Equi-dosing regimen region (EDRR) for transit compartment model with F ¼ 1, k e =0.0692 h À1 , k a =0.7 h À1 , and hypothetical minimum effective and maximum safe drug levels D me ¼ 300 mg and D MS ¼ 1000 mg. In both plots, the solid light blue EDRR is that for standard oral dosing with no transit compartments, and the four EDRR boundary curves are for transit compartment models with varying number of transit compartments n. a Varying n with fixed mean transit time MTT ¼ 4:4 h, so that k ¼ n MTT . b Varying n with fixed transit rate constant k ¼ 0:45 h À1 , so that MTT ¼ n safety and effectiveness combinations. A computed EDRR with illustrative drug level time courses is shown in Fig. 10.

Equi-dosing regimen region for oral equi-bolus dosing with n transit compartments
For the transit compartment model (Mt) with forcing (Beq), we consider the steady-state solution (3.28). Locating extrema in the time course is no longer viable analytically; simple inequalities such as in (5.3) and (5.6) cannot be found. Instead, to construct the EDRR, we discretise the ðT; D 0 Þ-space, and compare numerically-found maxima and minima of a c ðtÞ with D MS and D me respectively. The numerical method for constructing transit compartment EDRRs is summarised algorithmically and compared with the analytical approach for the simpler models in Table 3.
In Fig. 11, we show EDRRs for transit compartment models with varying number of transit compartments n, while fixing either the transit rate constant or the mean transit time (which effects the ''smoothed delay''). In each case, the EDRR is petal-shaped. Since the smoothed delay gives a narrower band for the steady-state timecourse than pure time-lag, the transit compartment EDRRs always contain the standard oral-dosing EDRR as a subset. Increasing either the timing of the delay (by decreasing n with k fixed) or its ''spread'' (by increasing n with MTT fixed) results in extension of the EDRR. For safe, conservative dosing regimen design, any regimen within the standard oral-dosing EDRR can, of course, be chosen.

Three-parameter dosing regimens
While equi-bolus dosing (input (Beq), with parameters T and D 0 ) is often of interest, we note that the cases of equibolus dosing with loading dose (input (BeqL)) and equiinfusion dosing (input (Ieq)) are also important clinically [43]. Appropriate three-dimensional EDRRs for these three-parameter regimens for IV administration may also be constructed using the analytical results.

Equi-dosing regimen region for IV equi-bolus dosing with loading dose
For model (M1) with forcing (BeqL), we seek constraints on the three equi-dosing regimen parameters D 0 , T and D L such that the time course given by (3.4) gives drug level that is safe and therapeutic immediately and always. The peak and trough levels are monotonic in time, so the necessary and sufficient constraints are that the steady-state time course given by For any given loading dose, the steady-state constraints require T and D L to be within the petal-shaped two-parameter EDRR as before. The added constraints (5.8) limit the dosing interval duration such that T\ 1 k e log FD L D me : Thus, for fixed D L , the two-parameter (T; D 0 ) dosing regimen region is a ''chopped petal'' shape, as shown in Fig. 12. The three-dimensional EDRR is given by the union of all such chopped petals, with D L ¼ D me F e k e T as a boundary surface.
A computed three-parameter EDRR with illustrative drug level time courses is shown in Fig. 13. The upper boundary surface projected onto the ðT; D 0 Þ plane gives the two-parameter EDRR for IV dosing (see further visualisation in Appendix 3). Only regimens (b) and (c) are within the three-parameter EDRR. Regimen (d) gives time course which is not therapeutic immediately, but it is at steadystate; it is clear that addition of a larger loading dose as in (c) gives an acceptable regimen. Practical dosing considerations such as frequency of administration, and therapeutic considerations such as drug level fluctuation, dosing interval averages and safety margins, can lead to specific dosing protocols [43]. Regimens following such protocols can easily be found by exploring the EDRR intuitively, for example, taking into account proximity to EDRR boundaries.

IV infusion equi-dosing
For model (M1) with forcing (Ieq), we have the three parameter dosing regimen ðT; k in ; t f Þ, and seek constraints for safe and therapeutic drug levels at steady-state. From (3.8), since the maximum and minimum drug levels at steady-state occur at t 1 ¼ t f and t 1 ¼ 0 respectively, we readily find that : ð5:10Þ The acceptable three-dimensional equi-dosing regimen region (for t f \T) is therefore given by D me 1 À e k e T 1 À e k e t f \ Fk in k e \D MS 1 À e Àk e T 1 À e Àk e t f : ð5:11Þ It is instructive first to consider a two-parameter ðT; k in Þ EDRR, given fixed infusion off-time t f . We find another chopped petal two-dimensional region, now as illustrated in Fig. 14. It is straightforward to show that the cross-over value of infusion rate is The petal does not extend to the origin in the ðT; k in Þ-plane; it is chopped at at T ¼ t f due to the simple constraint that t f \T. The three-dimensional EDRR is given by the union of all such chopped petals as t f varies, bounded by the planes T ¼ t f and T ¼ t f þ 1 k e log D MS D me . In Fig. 15, we show a computed three-parameter EDRR for a range of t f values, with illustrative drug level time courses. It is clear that dosing regimens can easily be chosen from within the EDRR. Only regimens (a) and (b) are within the threeparameter EDRR. Furthermore, systematic, direction-based regimen perturbations within the EDRR can be made to affect time course properties such as fluctuation level. The EDRR is proposed as a useful summary towards dosing regimen design.

Discussion
We have derived new analytical solutions for a generalised multi-dose transit compartment model (TCM), extending the analysis of the popular model in [44]. These solutions provide a means for analysing and parameterising delayed drug-level time courses without the need for nonsmooth time-lag models. The smoothed delay profile may be seen to give a better fit to experimental data than pure time delay models. The generalised model allows for simulation of realistic repeated dosing regimens that have traditionally been analysed in detail for simpler one-and two-compartment models [10,43,45]. In this sense, we importantly bridge the gap between traditional multi-dose analysis and the time-delay and TCM literature [23, 26-28, 38, 44, 55], providing a powerful method for capturing delays.
The exact solutions for the multi-dosing TCM will serve primarily as a tool for PK analysis. Further complexity may also be added as a future modelling extension by considering the ''body'' as a two-compartment schematic including central and peripheral compartments. Since firstorder transfer is typically considered between these compartments [20,43], the resulting ODE system including the full transit compartment cascade will be linear. As such, we expect analytical solutions to be available again, via the Laplace Transform method. Our new TCM analysis may also have application beyond PK. For example, signal transduction dynamics can often be modelled with linear transit compartment cascades [15,30,47,54]. Recently, solutions comprising incomplete gamma functions have been found for linear signal transduction cascades under different input conditions [5]. We expect that further analysis of periodic impulsive inputs to such systems will be valuable.
Drug dosing regimen design is an important consideration in therapeutics, from the stage of drug development [39] through to personalised regimens [14,36]. Given PK parameters, prediction of drug levels based on regimen parameters is common. With analytical expressions for drug levels as functions of time, we have shown that constraints on dosing parameters are readily available at the outset of any simulation-based study. Furthermore, the corresponding equi-dosing regimen regions (EDRRs) provide a novel, clear and succinct summarising visualisation of the acceptable dosing regimen parameter space, which may be explored intuitively to design effective and nontoxic treatments.
Predictive modelling using ODE models for PK is common, with end-user pharmacologists widely using exact solutions for low-dimensional compartmental models via a range of computational tools. Rapid computation of predicted time courses for multiple dosing regimens has further been facilitated by the development of usedfriendly simulation packages (e.g. [16,18,22]). EDRR visualisation can easily be achieved through a variety of computational tools, and we suggest that EDRRs could easily be incorporated into a number of packages to aid regimen design studies.
Finally, we remark that the mathematical detail of our work is also interesting in its own right, under the banner of mathematical pharmacology, which is now a recognised and growing field [51]. In [29], simple compartmental PK models are proposed as a starting point for biomathematics study and research. We propose a number of model perturbations and related mathematical directions beyond the scope of the current work. Constructing exact solutions for the TCM model relies on evaluation of the lower incomplete gamma function cðn; tÞ, which we have explored in some detail. An assessment of the practicality of using the analytical TCM results here versus numerical ODE solutions would be useful. Incorporation of more efficient and accurate approximation for cðn; tÞ, especially for parameter estimation purposes, may be a valuable pursuit, as discussed in [1,2,7,50]. We have proposed a thorough practical and sturctural identifiability analysis of our TCM ('Transit compartments-smoothed delays, lag time and data fitting (single-dose)' and Appendix 2.3). A theoretical comparison of the new TCM results and delay-differential equation (DDE) modelling approaches (see [26,38]) under impulse train inputs is also warranted. Linear pharmacokinetics is studied in our work, and much of the solution method relies on time-invariance of the PK parameters. However, chronopharmacokinetics is an important phenomenon that should be considered in PK ODE modelling [8,21,52]. Extension of our models and methodology to incorporate time-dependent parameters will be explored in future, but may be limited to numerical computation. Further, wider applicability of the TCM and EDRR methods will be achieved by consideration of nonlinear Michaelis-Menten elimination, which is discussed mathematically in [49,53]. Also, importantly, the PK models here may be linked to pharmacodynamics (PD) models to explore predicted drug responses; PD models are described in detail in [20,25,43,45].
Acknowledgements FH was supported by an Erasmus traineeship agreement between Universität Ulm and Swansea University. LB was supported in part by a Vice Chancellor's Early Career Researcher Award at UWE Bristol. We thank the expert reviewers for their valuable comments, which have resulted in an improved manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
Author Contributions Both authors devised the project, developed models, performed analysis and computation, and wrote the manuscript.

Appendix 1: multi-dosing solutions
Equi-dosing solutions are derived here. While the results are well known for the one-compartment and two-compartment models, it is instructive to see their solutions derived by Laplace Transform methods prior to deriving the full transit compartment model solution. Further, the steady-state solutions are vital for construction of EDRRs.

Single compartment-equi-bolus dosing
The solution to the (M1,Beq) problem (see Table 1 Taking the inverse Laplace Transform, we find the solution for the drug amount in the central compartment as where H is the Heaviside function. Since The loading dose effect is transient, so the steady-state behaviour is unaffected by D L .

Single compartment-equi-infusion dosing
The solution to the (M1,Ieq) problem (see Table 1 e ÀðjÀ1ÞTs À e ÀððjÀ1ÞTþt f Þs sðs þ k e Þ : Using the inverse Laplace Transform result [56] L À1 e Àas sðs þ k e Þ & ' ¼ Hðt À aÞð1 À e Àk e ðtÀaÞ Þ; we find the solution for the drug amount in the central compartment as Hðt j Þð1 À e Àketj Þ À Hðt j À t f Þð1 À e ÀkeðtjÀtf Þ Þ; ðA:12Þ where t j ¼ t À ðj À 1ÞT. Again, a geometric series results, and for 0 t M \T, we may express the drug level as ðA:13Þ Letting M ! 1, we see that the steady-state (T-periodic) level is given by e k e t f À e k e T 1 À e k e T e Àk e t 1 & ÀHðt 1 À t f Þð1 À e Àk e ðt 1 Àt f Þ Þ o ; for 0 t 1 \T:

ðA:14Þ
Alternatively, we may express this as e k e T À e k e t f 1 À e k e T e Àk e t 1 þ e Àk e ðt 1 Àt f ÞHðt 1 Àt f Þ & ' :

ðA:15Þ
We note here that the standard result (see, e.g., [45]) for a single continuous infusion is found by taking M ¼ 1 in (3.6) to give a c ðtÞ ¼ Fk in k e ð1 À e Àk e t Þ: ðA:16Þ

Two compartments-equi-bolus dosing
With a view to generalising the model to multiple transit compartment absorption, it is convenient to write the system (M2) (from Table 1) in matrix form as :

ðA:17bÞ
We solve the matrix ODE problem again using the method of Laplace Transforms. Firstly, we have We readily find ðsI À BÞ À1 ¼ 1 ðs þ k a Þðs þ k e Þ s þ k e 0 k a s þ k a ; so that the Laplace Transform of the solution is Taking the inverse Laplace Transform, we find the solution Hðt j Þ e Àk a t j À e Àk e t j Â Ã 0 B B @ 1 C C A : ðA:20Þ

Transit compartments-equi-bolus dosing
The solution to the (Mt,Beq) problem (see Table 1) can be found using Laplace Transforms. Taking the Laplace Transform of the IVP (2.6) (with g 1 ðtÞ ¼ g B ðtÞ) gives where the only nonzero element of GðsÞ is G 1 ðsÞ ¼ G B ðsÞ. Now, ðsI À BÞ is lower bidiagnoal, given by: ðA:22Þ Hence (A.21) may be solved easily by a process of forward substitution, or by noting that ðsI À BÞ À1 is lower triangular, given by: The Laplace Transform of the drug level a i ðtÞ in the i th transit compartment is then given by The Laplace Transform of the drug level a b ðtÞ in the absorption compartment is where, as before, For the transit compartments, we find from (A.24) and (A.27) that Now, given that L À1 e ÀðjÀ1ÞTs ðs þ kÞ i ( ) ðA:28Þ we find that Setting s j ¼ s À ðj À 1ÞT, we readily compute the following inverse transform as a convolution: k nÀ1 ðs þ kÞ n k nÀ2 ðs þ kÞ nÀ1 k nÀ3 ðs þ kÞ nÀ2 Á Á Á 1 s þ k k n ðs þ kÞ n ðs þ k a Þ k nÀ1 ðs þ kÞ nÀ1 ðs þ k a Þ k nÀ2 ðs þ kÞ nÀ2 ðs þ k a Þ Á Á Á k ðs þ kÞðs þ k a Þ 1 s þ k a k n k a ðs þ kÞ n ðs þ k a Þðs þ k e Þ k nÀ1 k a ðs þ kÞ nÀ1 ðs þ k a Þðs þ k e Þ k nÀ2 k a ðs þ kÞ nÀ2 ðs þ k a Þðs þ k e Þ Á Á Á kk a ðs þ kÞðs þ k a Þðs þ k e Þ k a ðs þ k a Þðs þ k e Þ Hðs j Þs nÀ1 j e Àks j e Àk a ðtÀsÞ À e Àk e ðtÀsÞ ds; Hðs j Þs nÀ1 j e Àks j e Àk a t e k a s j þðjÀ1ÞT Àe Àk e t e k e s j þðjÀ1ÞT ds j ; e Àk a t e k a s j þðjÀ1ÞT À e Àk e t e k e s j þðjÀ1ÞT ds j ; ðA:39Þ

Steady-state behaviour
To complete the analysis (and with a view to constructing EDRR's), we require solutions for the steady-state time course profiles. For each of the transit compartments (i ¼ 1; . . .; n), one approach is to solve the ODE for 0 t 1 \T, to give a solution parameterised by the initial value a 1 i ð0Þ (where t 1 ¼ 0 corresponds to the timing of the bolus input to a 1 ), then use periodicity to determine the appropriate value of a 1 i ð0Þ. In the following, we use lower case for the state variables and upper case for Laplace Transforms. For clarity and illustration, we derive the solutions for the first five transit compartments before writing the solution for transit compartment n in general.

ðA:49aÞ
The coefficients a i ð0Þ (the steady-state dosing interval initial values) may be found using together with the recurrence relation (for i ¼ 2; . . .; n) p! / p :

#)
; ðA:51gÞ ðk a Àk e Þð1Àe ÀkeT Þ : ðA:51hÞ Computational evaluation of the lower incomplete gamma function Our preferred method for evaluating the lower incomplete gamma function is using the built-in function gammainc in MATLAB, which computes the normalised function Pðn; tÞ ¼ cðn; tÞ CðnÞ . In Fig. 16, we show a typical result comparing run times and computed values of cðn; tÞ for a given t and range of n, for five different methods. For computations in MATLAB [35], it is clear that the quickest method is using the built-in command gammainc. All methods give values in agreement, except for the symbolic computation of the sum in (3.17) for large n. This is due to rounding and truncation errors introduced when evaluating computed symbolic expressions. Hard-coding the expressions in this sum results in the second shortest run time, but this method is not practical for n [ 10. Use of the log-gamma function is an order of magnitude slower, but is still much faster than employing a numerical integration method. We conclude that built-in functions, either evaluating cðn; tÞ directly or via the log-gamma and distribution functions, represent a practical approach to evaluating cðn; tÞ. More efficient methods involving series and continued fraction expansions may be investigated [1,2,17,41,50], but are beyond our scope here.
time since the first dose is t end (fixed constant). Then t end ¼ MT, and so a c ðt end Þ ¼ a M c ðTÞ ¼ FR 1 À e Àk e t end 1 À e Àk e T Te Àk e T : Then taking infinitesimally small dosing interval (and so an infinitesimal dose), our result becomes a c ðt end Þ ! lim FRð1 À e Àk e t end Þ Te Àk e T 1 À e Àk e T : Now, using l'Hopital's Rule, we see that a c ðt end Þ ! FRð1 À e Àk e t end Þ lim T!0 e Àk e T ð1 À k e TÞ k e e Àk e T : Thus, as T ! 0, we have a c ðt end Þ ! FRð1 À k e t end Þ Â 1 k e ¼ FR k e ð1 À e Àk e t end Þ:

Two-compartment oral equi-bolus dosing
Here, we note some properties of the two-compartment equi-bolus dosing solutions (3.10) and (3.11). Firstly, the minimum drug level on the nth dosing interval is given by a n c ð0Þ, while for a local maximum (a peak level) in the time course for a n c ðt n Þ at t n ¼ t Ã n , we need a n c 0 ðt Ã n Þ ¼ 0. We find, by using (3.10), that: 1 À e Ànk a T 1 À e Àk a T 1 À e Àk e T 1 À e Ànk e T :

ðB:2Þ
If 0\t Ã n \T, then there is a peak blood drug level in the nth dosing interval, given by   h. For validation, we show results from analytical solution and a numerically computed solution. The dosing interval average at steady-state is indicated. Also, the steady-state profile is shown for each compartment a n c ;max ¼ a n c t Ã n À Á ¼ FD 0 k e k a ke 1 À e Àk a T 1 À e ÀnkaT k e 1 À e Ànk e T 1 À e ÀkeT k a ( ) 1 ka Àke : ðB:3Þ So we may write the maximum and minimum values of a n c on [0, T] as: a n c ;min ¼ a n c ð0Þ ¼ FD 0 k a k a À k e 1 À e Ànk e T 1 À e Àk e T À 1 À e Ànk a T 1 À e Àk a T

& '
; ðB:4Þ a n c ;max ¼ a n c minðt Ã n ; TÞ À Á ; (using (B.2)) : ðB:5Þ At steady-state, a local maximum in the time course is ensured, occurring at t Ã 1 ¼ 1 k a À k e log k a k e 1 À e Àk e T 1 À e Àk a T : ðB:6Þ The steady-state minimum and maximum drug levels are then given by a 1 c ;min ¼ a 1 c ð0Þ ¼ FD 0 k a k a À k e 1 1 À e Àk e T À Transit compartment parameter identifiability In Fig. 17, we show time courses for a transit compartment model with n ¼ 1 with a single bolus dose. Non-identifiability of k and k a is clear, due to the flip-flop phenomenon discussed in 'Transit compartments-smoothed delays, lag time and data fitting (single-dose)'. For our transit compartment model, this problem is unique to n ¼ 1, and no such issues are seen with n [ 1. In Fig. 18, we plot a typical parameter estimation result for n ¼ 10 compartments, using the method described in 'Transit compartments-smoothed delays, lag time and data fitting (singledose)'. With ''perfect'' pseudo-data, we see the base parameter set being recovered exactly, and with noisy pseudo-data, the parameter estimates are in good agreement with the base set. A formal structural and practical identifiability analysis is proposed as future work, and will use transfer functions and other methods described in [13].

Transit compartment equi-bolus dosing
In Figs. 19-20, we show time course results for a transit compartment model under multi-dosing input (Mt,Beq), using parameters taken from the fitted data in Fig. 5. Each figure shows the analytical results for drug levels in a number of transit compartments, with n ¼ 10 transit compartments in total. The absorption and central compartments are also shown. For validation, results from a numerical ODE solver are also shown. For tracking the approach to steady-state, we mark the steady-state profile and dosing interval average for each compartment plotted. For each of the transit compartments, the steady-state dosing interval average may be computed, for example, from (2.6) by integrating each of the compartment ODEs in turn. We find the transit compartment averages to be In both figures, we observe the delay effect, as the peak drug levels appear later further through the transit cascade. In Fig. 19, the transit rate constant k is taken from the earlier data fitting, and the system almost reaches steadystate by the second dosing interval. In Fig. 20, we take k four-fold lower, which increases the smoothed delay, and hence delays the approach to steady-state. In this case, the delay effect is more evident throughout the transit cascade.