Comparison of the gamma-Pareto convolution with conventional methods of characterising metformin pharmacokinetics in dogs

A model was developed for long term metformin tissue retention based upon temporally inclusive models of serum/plasma concentration (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ C $$\end{document}C) having power function tails called the gamma-Pareto type I convolution (GPC) model and was contrasted with biexponential (E2) and noncompartmental (NC) metformin models. GPC models of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ C $$\end{document}C have a peripheral venous first arrival of drug-times parameter, early \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ C $$\end{document}C peaks and very slow washouts of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ C $$\end{document}C. The GPC, E2 and NC models were applied to a total of 148 serum samples drawn from 20 min to 72 h following bolus intravenous metformin in seven healthy mongrel dogs. The GPC model was used to calculate area under the curve (AUC), clearance (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ CL $$\end{document}CL), and functions of time, f(t), for drug mass remaining (M), apparent volume of distribution (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{d}$$\end{document}Vd), as well as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{1/2}\ f(t)$$\end{document}t1/2f(t) for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ C $$\end{document}C, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ M $$\end{document}M and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{d}$$\end{document}Vd. The GPC models of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ C $$\end{document}C yielded metformin \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ CL $$\end{document}CL-values that were 84.8% of total renal plasma flow (RPF) as estimated from meta-analysis. The GPC \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ CL $$\end{document}CL-values were significantly less than the corresponding NC and E2 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ CL $$\end{document}CL-values of 104.7% and 123.7% of RPF, respectively. The GPC plasma/serum only model predicted 78.9% drug \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ M $$\end{document}M average urinary recovery at 72 h; similar to prior human urine drug \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ M $$\end{document}M collection results. The GPC model \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{1/2}$$\end{document}t1/2 of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ M $$\end{document}M, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ C $$\end{document}C and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_d$$\end{document}Vd, were asymptotically proportional to elapsed time, with a constant limiting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{1/2}$$\end{document}t1/2 ratio of M/C averaging 7.0 times, a result in keeping with prior simultaneous \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ C $$\end{document}C and urine \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ M $$\end{document}M collection studies and exhibiting a rate of apparent volume growth of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_d$$\end{document}Vd that achieved limiting constant values. A simulated constant average drug mass multidosing protocol exhibited increased \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_d$$\end{document}Vd and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{1/2}$$\end{document}t1/2 with elapsing time, effects that have been observed experimentally during same-dose multidosing. The GPC heavy-tailed models explained multiple documented phenomena that were unexplained with lighter-tailed models. Electronic supplementary material The online version of this article (10.1007/s10928-019-09666-z) contains supplementary material, which is available to authorized users.

Terminal apparent volume of distribution

Introduction
Metformin (1,1-dimethylbiguanide), average molecular weight 129.164 g/mol, is a þ1 cation at physiological pH with apparent volumes of distribution (dogs) [1] so large that it may be problematic to use classical pharmacokinetic methods to calculate those volumes. Metformin is used among numerous indications as the default first-line treatment for type II diabetes and for cancer chemotherapy in humans. There is a developing interest in treating feline and canine cancer with metformin, which in dogs is more typically used as an anti-hyperglycaemic agent to treat obesity or insulin resistance [1,2]. In humans, the postulated cancer cell susceptibility to metformin is from defective mitochondrial oxidative phosphorylation and low glucose levels [3]. Metformin is non-metabolised and exclusively plasma-phase renal cleared at an estimated 90-100% of total renal plasma flow in humans [4]. Seen in association with impaired renal function is dose related lactic acidosis, a rare but serious condition [5][6][7] suggesting a maximum safe dosing. Also, antineoplastic chemotherapy with metformin potentiates temsirolimus side effects [8], such that secondary metformin effects during therapy are also concerning. For example, extramitochondrial secondary metformin effects are suppression of glycogenolysis, and enhanced cellular glucose uptake relevant to metformin's pharmacodynamics in type II diabetic therapy [9]. During oral dosing, washed jejunal tissue biopsies had about 30-300 Â higher metformin concentrations than that in plasma (human) [10]. Metformin accumulates in the cytosol of erythrocytes, small intestine, skeletal and cardiac muscle cells, hepatocytes, and brain [4,[11][12][13]. Metformin's primary effect is to suppress gluconeogenesis in mitochondria [14]. Metformin mitochondrial concentration from cultured rat hepatoma cells was 1000 Â that in plasma at 60 h [15]. Metformin's extra-mitochondrial drug effects appear to occur in other intracellular locations. For example, control of gluconeogenesis was time-delayed using a same-dose regimen [16], and unrelated to blood concentrations [17]. As mitochondria have a very high affinity for metformin, and most eukaryotic cells (except mature erythrocytes) contain many mitochondria, which occupy up to 25 percent of the cytoplasmic volume [18], and are the sites of metformin tumoricidal drug action, with elapsing time one would expect a predominance of body drug mass to be in close proximity to the mitochondrial effector sites. Intracellular localised metformin effects and plasma drug concentration are kinetically distinct as supported by the observations in rats, dogs and humans that half-lives of drug mass in urine or in erythrocytes are approximately 4.2 to 11.1 times longer than in plasma [4,13,[19][20][21][22]. The difference of half-lives relates to persistent redistribution. To cause persistent redistribution, it would not be enough to add a solitary extra compartment to our model, but with elapsing time even more compartments would be needed. Variable-volume, variable half-life pharmacokinetic modelling is a relatively recent development that allows volume of distribution to be added to a model with elapsing time [23]. The additional concept of a half-life that is not static may seem novel, but was shown to arise in nonequilibrium states (i.e., when the mass half-life is longer than the plasma half-life), and occurs transiently before a compartmental model achieves dynamic equilibrium (i.e., when mass half-life and plasma half-life are equal). Thus, rather than just accept that metformin plasma concentration halflives are poorly characterised [11], we note that Xie et al. [4] as well as Sambol et al. [24], the latter for oral multidosing, have commented that metformin half-life increased with elapsing time. Indeed, bolus intravenous metformin plasma half-life was 1.5 h as last sampled at 8 h [20], 4.5 h as last sampled at 12 h [22], and 20.4 h, as last sampled at 72 h (dogs) [1], such that there is good evidence that metformin plasma half-lives increase with elapsing time.
To duplicate metformin half-lives of drug mass in urine and erythrocytes that are persistently multiple times longer than in plasma, classical drug models like the sum of exponential terms (SET) and noncompartmental (NC) exponential models are perhaps not the best choices as only a terminal drug mass to plasma concentration half-life ratio of one is possible for those models. SET functions yield transient nonequilibrium models whose time to dynamic equilibrium can be prolonged 7.5 times but not eliminated using gamma distributions with greater than exponential log-convexity enforced using adaptive Tikhonov regularization [23,25,26]. It is not surprising then that for many bolus drug experiments, a late-time exponential tail underestimates the actual amount of drug remaining in the body [5,27]. Many intravenously administered drugs has been observed to follow the very heavy-tailed, power-law at late times [28][29][30][31][32][33][34]. The statistical form for a power law is the Pareto distribution (PD), which like the Cauchy distribution, has tails so heavy they confer unusual statistical properties and have been given the name fat-tailed distributions. Power laws are scale independent and intrinsically fractal. Fractal drug models are consistent with carrier transport [35], and metformin is principally hOCT transported [7]. However, a power-law model does not accurately model the first few hours of concentration [36,37]. The incorporation of a power function tail into a semi-infinite support model, 0 t\1, was explored with fractional calculus [38], but without consideration of mass conservation and nonequilibrium dynamics. Explored is whether the inclusion of power function tail into a temporally more inclusive model of serum or plasma concentration will yield the permanent nonequilibrium that explains and duplicates the observed persistent mass to concentration ratio of half-lives and whether that forms an accurate enough concentration model to suggest clinical usefulness.
Convolution of a washout model (a monotonic decreasing function) with a gamma distribution (GD) has been used before for correcting the contribution to total residence time of early concentrations [37]. In the current paper, a gamma-Pareto type I convolution (GPC) concentration density function of time was used to model the pharmacokinetics of single bolus intravenous metformin injections is each of seven dogs, most having 22 timesamples collected over 20 min to 72 h. In a GPC model, the PD models the terminal tail of the drug's plasma concentration after the lighter-tailed gamma density (GD) has decayed. Hyperglycaemia therapy with metformin multidosing has shown a progressively increasing effect of decreasing fasting plasma glucose for at least the 8 weeks of testing [39]. Thus, the longer metformin drug mass in tissue may be better linked to drug effect than serum concentration [17], and the modelling introduced should generate accurate drug mass curves from serum data alone, and to be useful should allow for a multidosing regimen that rapidly establishes a constant drug mass in the body, as opposed to dosing for many weeks for that to occur. This may be relevant for cancer chemotherapy as the alternatives are either to wait for months for tissue build up of metformin, which may not provide optimal therapy for rapidly growing tumours, or to give large dosages without knowing when unnecessarily high metformin tissue concentration will occur.
The residual drug masses imply tissue retentions and excretions that were contrasted with published urinary excretion and erythrocyte retention results. The modifications of pharmacokinetic parameters resulting from the application of the new models are presented and discussed for bolus intravenous metformin disposition curves and body drug mass as functions of time and their half-life functions together with some distinct pharmacodynamic implications arising from them. As some of the results are compared to exponential tailed model results, some statistical characterisation and comparative data analysis of exponential based models is also provided.

Background
To obtain ratios of half-lives of mass and concentration that do not tend to unity from persistently non-equilibrated models required defining half-life as a function of time. Said ratios naturally arise in the context of variable apparent volume of distribution models, V d ðtÞ. The first mathematically correct V d ðtÞ models were proposed by Niazi for SET functions [40]. Those models were simultaneously compartmental models and variable volume models, but the interpretations of those different models are distinct. Niazi's work was generalised by some of us for any density function supported on the time is ½0; 1Þ interval having a varying apparent volume of drug distribution in time, V d ðtÞ, with half-life expressions that vary in time [23]. For variable volume modelling V d ðtÞ is a drug-ocentric volume of distribution in time. For metformin, the drug is either in the body, or has been eliminated to the urine. We do not know where the M(t) drug mass is in the body. We only know that mass is concentration times volume. That apparent volume of distribution is the imaginary volume that the drug would occupy if it were everywhere at the same concentration as it was observed to be in plasma at that time, and CL eliminates drug mass from that single volume no matter how big that volume is.
The general basic equations of variable volume modelling with variable half-life and constant clearance are listed as Eqs. (1)- (10). Note well that these equations are for bolus conditions; dose delivered rapidly under intravenous conditions, and do not include dose that remains in any syringe, tubing or vial, nor any depo caused by inadvertent subcutaneous extravasation, i.e., a common, but largely unrecognised problem [41]. Drug mass, M(t), means drug mass in the body, not drug mass excreted. Equations (1,2) require clearance CL to be proportional to plasma/serum concentration and constant in time. Strictly speaking, the observed concentrations, C obs ðtÞ, for these two equations should be measured where elimination is occurring, which for metformin is where the renal artery contents are cleared in the kidney, and the C(t) model fit to those observations. Note that Eq. (1) reduces to D ¼ CL AUC for t ¼ 0, where D ¼ Mðt ¼ 0Þ is the dose, CL is plasma/serum clearance and AUC is the area under the curve of C(t) from t ¼ 0 to 1. Equation (2) is a mass ratio; the quotient of Eq. (1) evaluated at t divided by its evaluation at t ¼ 0, is the complementary cumulative density function and defines a density function's right-hand tail area (pdf often have two-tails). As a matter of convenience, the complementary cumulative density function is referred to as a survival function SðtÞ ¼ def 1 À CDFðtÞ, where CDF(t) is the cumulative density function and CDFðtÞ ¼ def R t 0 pdfðuÞdu. This is common practice despite the confusion it causes. That is, a survival function in the strictest sense is discrete, which is not how the term is being used here. Although the name survival arose from actuarial usage as the fraction of an initial population surviving in time, here it is the fraction surviving in the body of a unit dose mass at time t. Sometimes, survival functions are also used to calculate ''Relative tail heaviness'' as per the Appendix Section. Survival functions and CDF are dimensionless (units cancel).
A list of equations for bolus-dose, constant-clearance, intravenous conditions follows: Mass as an f(t), Survival function definition, Conservation of mass, Concentration parsing, Two time-sample exponential half-life, Instantaneous half-life, Half-life of drug mass, Half-life of concentration, Ratio of M(t) to C(t) half-lives, Harmonic sum of half-lives, Equations (3)-(10) are variable apparent volume of distribution equations and/or variable half-life equations. Of these equations, Eq. (3) and the derivative of its logarithm, parsing of concentration Eq. (4), are the direct results of conservation of mass. For these equations, V d ðtÞ at t ¼ 0 is either undefined for a short time when C(t) models a peripheral venous sample or a C(t) is used that models more proximal concentrations at the injection site itself. Equation (4) relates that the negative relative change in concentration (i.e., a positive number for decreasing concentration) is the negative relative change of mass (also typically positive) plus the relative drug volume of distribution apparent growth, which latter dilutes/rarefies the relative concentration. For a constant volume of distribution, V d-RT ¼ CL RT can be used, where RT is whichever type of residence time is appropriate to the model under consideration-see the ''Residence Time'' Appendix Subsection.
Equation (5) is common knowledge for any two timesamples, ft i ; C i g and ft iþ1 ; C iþ1 g. Use of this equation for non-compartmental conditions when t iþ1 À t i is small is done with median values from all of the dogs' samples to obtain tractable results-see the ''Noncompartmental timesample group, median half-life'' Appendix Subsection for details. However, Eq. (5) only applies for monoexponential conditions or in the limit if we substitute f ðtÞ % C i , and t iþi ¼ t i þ Dt and let Dt ! 0 þ , where f(t) is any continuously differentiable function. That yields Eq. (6), which is of classical mechanical type, and applies equally well to the evolution of half-life of satellite orbital decay [42] as to variable half-life of concentration of any drug. Half-lives can be both positive and negative. For an f(t) that increases in time, the t 1=2 is negative. For example, this occurs when concentration is increasing in time in a peripheral vein shortly after intravenous injection. It easier to understand a negative half-life of an increasing function of time when some maximum total value of a measurement is being approached in time by subtracting the measured values from that maximum and viewing the result as a positive half-life of a decreasing function in time, as is done, for example, using total dosage administered for a urinary drug recovery t 1=2 calculation. However, growth of any function in time, including drug mass accumulating in urine, actually has a negative half-life.
For the half-life of drug mass from the Eq. (6) definition, we substitute Eq. (3) for f(t), and the derivative of Eq. (1) for f 0 ðtÞ, which yields Eq. (7). To find the plasma/ serum clearance half-life, Eq. (8), we eliminate the mass terms in concentration parsing Eq. (4) using MðtÞ ¼ V d ðtÞ CðtÞ and M 0 ðtÞ ¼ ÀCL CðtÞ, then rearrange terms to put them into the form of Eq. (6). V 0 d ðtÞ is concentration dilution having the same units as CL, e.g., ml min À1 kg À1 . Next, we form the ratio of the mass and concentration halflives above to yield Eq. (9). We hypothesise that this can have different concentration parsing limits, [23] and Eq. (33) below, respectively. In the first case, V 0 d ðtÞ eventually shrinks to zero, volume dilution of concentration stops, and dynamic equilibrium is achieved. Thereafter, mass half-life is the same as the halflife in plasma/serum. In the second case, V 0 d ð1Þ is a positive constant and t 1=2 ; MðtÞ remains longer, e.g., as witnessed as Àt 1=2 of metformin mass in urine, than t 1=2 ; CðtÞ in plasma/serum. The half-life of V d ðtÞ adds harmonically, i.e., reciprocally, with mass and concentration half-lives, as follows. Substituting the mass and concentration half-lives, Eqs. (8,9), into concentration parsing Eq. (4) yields Eq. (10).

Constant infusion
Let us assume that elimination rate is proportional to concentration, and that superposition is valid. Using as impulse-response the infinitesimal of an intravenous bolus C(t) model and as constant input of that model, hðtÞ, the unit step function in the role of its transfer or control function, their convolution is a constant infusion concentration model. Then as hðxÞ Ã pdfðxÞ ðtÞ ¼ CDFðtÞ, where CDF is the cumulative density function, where the 'in' of C in ðtÞ is for 'infusion', and where C SS ¼ D R =CL is the steady state concentration [23]. D R , an acronym of dosing rate, was used and elsewhere R 0 or k 0 are often used. A CDFðtÞ has a maximum amplitude of 1 terminally; the occurrence time of a steady state concentration, C SS . Thus, to scale CDFðtÞ, we multiply it by C SS the magnitude of which is R 1 0 CðxÞdx from the above convolution, such that C SS AUC bolus . Note that there is no bolus physically, just a result that can predict the magnitude of C in ðtÞ in terms of a previously observed bolus function.

Methods
This section includes new material such as constant average drug mass multidosing and the explanatory statistical treatment of models as scaled density functions.

Data collection
The metformin serum concentration data was published elsewhere [1] was reused to reexamine basic modelling assumptions. In that study, seven healthy adult mixedbreed dogs were used including 2 sexually intact males, 2 neutered males, 1 sexually intact female, and 2 spayed females. The dogs were between 2 and 3 years old and weighed between 25.7 and 29.2 kg. The animal study was approved by the University of Saskatchewan's Animal Research Ethics Board, and adhered to the Canadian Council on Animal Care guidelines for humane animal use.

Fitting metformin concentrations
Area under the curve (AUC) inherits its units from those of the curve being used to calculate that area. For example, density functions have a total AUC pdf of 1, i.e., R 1 0 pdfðtÞ dt ¼ def 1, with no units (dimensionless) whereas for concentration AUC CðtÞ has units like mgÁh L . The acronym pdf was used for density functions even though the p for probability was irrelevant. That is, concentration is not a probability, and a fraction of total AUC per unit time at time t density function, i.e., pdfðtÞ ¼ CðtÞ=AUC CðtÞ (e.g., units h À1 ), obeys similar rules to those that would apply to any other density per x-axis unit function (pdf) in whatever context the x-axis units and (x, y) coordinated values provide. In our case, fitting data with CðtÞ ¼ j pdf ðtÞ yields a concentration curve model, C(t), where AUC CðtÞ ¼ j.
Peripheral venous drug concentration is observed, C obs ðtÞ, to be smoothed and delayed by its passage through the circulatory system causing zero initial concentration, C obs ð0Þ ¼ 0, to be observed peripherally. In this paper, a virtual washout signal with a t ¼ 0 start time, i.e., a gamma distribution, is convolved with a Pareto distribution having a 25-30 s start time to yield a peripheral C(t) model. Ideally, one should minimise bolus error; to not leave drug in the needle track within the dermis by using a large flush bolus injectate to wash the needle and to deliver a highconcentration, low drug volume quickly enough to improve and standardise circulatory mixing [41].
To select disposition curve models one needs to be mindful of which statistical distributions have the best chance of being useful. These selection criteria include inspection of each distribution's parameter types. A distribution's parameter types influence the distribution's derivative fidelity for following the curve shape, and influences the well-posedness of the AUC determination, i.e., whether a unique (global) solution that is not sensitive to noise, called a stable solution exists. Stable solutions were obtained in high precision using the Nelder-Mead global optimisation method [43], The practical choice for the modeller is to either find a robust concentration model whose fit to the curve is accurate or to relax the goodnessof-fit criteria, for example by treating AUC as an ill-posed integral of the first kind with use of an inverse method for fitting [25]. Distributions generally take one or more of three types of parameters, those of shape, location and scale. For example, the exponential distribution has a scale parameter, but no shape or location parameters. A location parameter indicates where a distribution is positioned on the x-axis. Depending on the distribution, the location may be a best measurement of centrality, but more generally just indicates an x-axis position. Compared to a distribution without a shape parameter, if we choose a distribution with one or more shape parameters, not only will our goodness of fit generally improve, but the curve shape will be better fitted. Some functions, like the gamma distribution (GD), which has a shape and a scale parameter, may still have an ill-posed AUC, such that fitting with Tikhonov regularization rather than ordinary regression may be needed to find a Tikhonov well-posed clearance at the price of a small penalty in goodness-of-fit [25].
To fit concentrations we used proportional error minimisation, i.e., 1=C 2 obs ðtÞ weighted ordinary least squares fitting of models to data [37,[44][45][46][47]-see the ''Concentration data and fitting it'' Subsection of the Appendix for details on this and on proportional assay error, how to quantify root mean square proportional error (rrms), and noncompartmental methods. The GPC function was the experimental model fit to data. It is calculated from hundreds of sums using 65 decimal place arithmetic to mitigate roundoff error. The GPC algorithm with Nelder-Mead regression was designed for precision and accuracy to converge to 30 decimal places. The GPC algorithm used, Eq. (25) below, purpose written in Mathematica 12.0.0.0 code, was not optimised for speed such that the run time for each case was 34 AE 12 min on a 2.7 GHz Intel Core i7 CPU using four parallel processors. In post hoc testing this was reduced to 28 min using series acceleration methods, i.e., not much faster and that despite reducing the call to the GPC algorithm from 45 to 19 ms average. Thus, much of the elapsed time for regression is due to factors other than GPC algorithm run time, for example performing constrained rather than unconstrained regression and the use of extended precision. Individual subject parameter statistics (CV%) are not a feature of the global optimisation search method used, and see the ''Parameter errors from modelbased bootstrap'' Appendix Section for a computationally intensive analysis of algorithmic robustness. For E2, parameter CV% are typically calculated from local optimisation gradients, e.g., using the Levenberg-Marquardt method and similar methods, which may not be the CV% from least error, i.e., global, solutions, and only indirectly examine for robustness. Sample-time groups are used in the text, especially in ''Noncompartmental time-sample group, median half-life''. These are merely groups of samples having in common the same times following bolus injection that were used for drawing samples from all of the subjects.

Dose regimen models
The simplest dose regimen is repeat constant dosing. Another dose regimen of possible interest for vasoactive drugs, e.g., dipyridamole, would be to keep the mean plasma drug concentration constant for each dose interval. A dose regimen for drugs like metformin that have major effects in cytosol and mitochondria and that accumulate in cytosol, may have as a goal to dose for the same mean drug mass retained in the body during each dose interval.
Let us assume first order elimination and superposition. Then, a protocol for keeping a constant mean unit dose drug mass in the body starts with a unit intravenous dose from the identity function for a density, R 1 0 pdfðtÞ dt ¼ 1, break it into dose periods of s, and divide that series by s to calculate the average fraction of a unit dose mass per unit time during each dosing interval, which is the average mass over the dose interval, which written from the survival function equals DSðtÞ=s, i.e., SsðiÞ ¼ 1 During the first dose interval of duration s, to augment the dosage to yield an average drug mass during s to be a unit dose, the dose administered is augmented by multiplying the unit dose by the reciprocals of the first terms of right sides Eqs. (12a) and (12b) to yield the dose factor, DF To obtain a dose factor for a constant average unit dose mass during the second s interval one includes the contribution from the first dose second dose interval, and the second dose first dose interval, which simplifies to yield For n ! 2, this equation generalises to be the recursion, Note that as n ! 1, DFðnÞ goes to that fractional dose multiplier, [ 0, that allows an average of a unit dose to be maintained in the body during the n th dose interval. A word of caution, for accurate results from Eq. (14), the terminal tail of the model used must accurately match plasma/serum drug concentrations.

Distributions
This subsection presents the distributions (pdf) used. Some properties of these distributions are: The convolution of two pdf is a pdf. The AUC of a pdf is 1. CðtÞ ¼ j pdf where j is the AUC of C(t). A survival function, S(t), is the t to 1 integral of the pdf. For intravenous bolus conditions, regardless of distribution type, when t ¼ 0, Sð0Þ ¼ 1, and the entire dose fraction is within the body.

Gamma distribution
The gamma distribution (GD) is given by where the gamma function satisfies CðaÞ ¼ R 1 0 e Àt t aÀ1 dt, and hðtÞ is the unit step function, i.e., 0 for t\ 0 and 1 for t ! 0. The GD is an exponential density (ED, below) when a ¼ 1. The GD also has a (þ 1) discontinuity at t ¼ 0 when 0 \ a \ 1. However, that discontinuity is integrable with zero area in the limit as t ! 0. Thus, the survival function (S GD ) is defined at t ¼ 0, and, where the upper incomplete gamma function satisfies Cða; zÞ ¼ R 1 z u aÀ1 e Àu du, and Q(a, z) is its regularised form. The GD rate parameter is b, whereas 1 b is the scale parameter. The shape parameter for the GD is a. The shape parameter aids in fitting disposition curves and their shapes. There is no location parameter.

Pareto distribution
Negative power functions of time have an undefined (þ 1) discontinuity at t ¼ 0 that is not integrable (þ 1). Thus, a power function as a density function could have a minimum positive start time, b [ 0, thereby allowing for integration that avoids the discontinuity. Just such a function is the type I Pareto distribution (PD), Note the shift to the right of the unit step function; hðt À bÞ, i.e., the PD function is zero for t\ b. The PD survival function (S PD ) is, b is the scale parameter. The shape parameter is a. There is no location parameter.

Sums of exponential term distribution
Sums of exponential terms (SET) models are equivalently C SET ðtÞ ¼ hðtÞ P n i¼1 c i e Àk i t , where n is any positive integer including 1, and without loss of generality k 1 ! k 2 ! k 3 ::: ! k n [ 0. SET models have been incorrectly criticised as having meaningless c i coefficients. In a statistical context note that ke Àkt is an exponential density, then let P n i¼1 p i ¼ 1, where the p i (not probabilities here) are the fractions of the model's total density attributed to each exponential density term. Then, where the k i decay coefficients are rate (i.e., 1/scale) parameters. The c i of SET functions relate to the ED n parameters as Note that a scaled monoexponential, an E1, is not a monoexponential density, ED, as E1 ¼ AUC E1 ED. So when we refer to a biexponential we write E2, where a biexponential density is an ED 2 . The survival function is Note that ED n , S EDn and SET only have scale parameters (or their reciprocals). There are no shape parameters to aid for fitting curve shapes or for extrapolation, and the solutions for SET n ! 2 lack robustness. For example, biexponentials (E2) solved for 4 time-samples sometimes have solutions as 2 C, not 2 R [36]. There are no location parameters.

The gamma-Pareto convolution distribution, basic formulas
Derived from series expansion of the GD exponential, Eq. (15), followed by convolution with a PD, Eq. (17) and summation of the expanded parts, 1 the GPC model inherits two shape parameters, a and a, and two scale parameters, b and b from its convolved pdf's, where the incomplete beta function is B z ða; bÞ ¼ R z 0 u aÀ1 ð1 À uÞ bÀ1 du. The GPC function has a concentration start time, b [ 0, and a circulatory mixing start time of 0. In practice, GPC models chain, i.e., convolve, two independent washout functions (monotonic decreasing; GD iff a 1 and a PD) that only when combined produce a circulatory peak concentration. The GD washout models relate to the circulatory mixing that starts at t ¼ 0, whereas the PD models relate to the injection to sampling venous-venous delay of a few dozen seconds (¼ b) and the terminal rate of the increase in apparent volume of distribution in time [CL=a, Eq. (33), below].
For use with M(t), and other functions, we integrate the GPC function from 0 to t to construct its cumulative density function (proof by differentiation), where for z ¼ 1, B z ða; bÞ becomes the beta function Bða; bÞ ¼ CðaÞ CðbÞ CðaþbÞ . As it is needed for the half-life of concentration, we take the derivative of the GPC, Asymptotes of the gamma-Pareto convolution Two asymptotes to the GPC function were identified for t is a sufficiently long time. The first is the Pareto Type I distribution itself, and is the desired result; to write a power function tail into a more general PK model, and is a straight line on a log-log plot, where ( $ ) means asymptotic to. This occurs because the GD tail decays quickly and the PD tail decays very slowly, see ''Relative tail heaviness''. The second GPC asymptote was obtained by Taylor series expansion of Eq. (21), with the proof appearing as the .pdf file labelled Supplementary Materials 2, ðaÞ k ðbÞ k z k k! is the regularised form of Kummer's confluent hypergeometric function of the first kind wherein ðaÞ k and ðbÞ k are Pochhammer sym bols, i.e., generalised ascending factorials, The closed form term (beginning þh) of Eq. (25) is positive when 0\a\1 and b t are sufficiently large, while the limit of the sum converges to zero. This forms a second GPC asymptote that approaches the GPC more rapidly than the PD tail, and finds use for computing the GPC when t is long, e.g., approximately [ 100 h. That is, for 0\a\1 and t sufficiently long, that term and the GCP are asymptotic to each other, What the concentration models imply

Drug mass
The drug mass remaining in the body, is from the GPC CDF, Eq. (22) applied to Eq. (2), where Mð0Þ ¼ D is the dose.

Half-life functions
The half-life function of metformin model concentration, t 1=2 ; CðtÞ, as observed in a peripheral venous sampling site, comes from substituting the GPCðtÞ, Eq. (21), and its derivative Eq. (23) for the f(t) and f 0 ðtÞ of Eq. (6), Next, the half-life of drug mass remaining in the body from Eqs. (2), (6), (21), and (22) is The asymptotic half-life of the GPC concentration model was obtained by substitution of the first GPC asymptote, Eq. (24), and its derivative into half-life Eq. (6), and is linear Similarly, the half-life of drug mass is asymptotic to a linear function of time. Equations (30) and (31) which as a [ 0, is [ 1. That this limit is not equal to 1, as it would be for lighter-tailed models, exemplifies why Pareto distribution tails have been called fat-tailed historically; because of their provocatively different properties from those of lighter-tailed distributions. 2 There are now two expressions for the mass to concentration ratio, one of them for any variable volume model, Eq. (9), and the other from the PD asymptote to the GPC, Eq. (24). Combining these expressions yields the terminal rate for volume of distribution apparent growth for our fat-tailed model of As a is dimensionless, V 0 d ; the rate of growth of V d , has the same units as CL. Rewriting of the half-life harmonic sum Eq. (10) and solving for the asymptote of t 1=2 ; V d ðtÞ yields its negative half-life (growth)

Results
The gamma-Pareto convolution (GPC) distributions' parameters for metformin in seven dogs appear in Table 1.
As a global search routine, the Nelder-Mead method does not generate case-wise CV% results. For that information, please refer to the ''Parameter errors from model-based bootstrap'' Appendix section. Figure 1 shows GPC models fit to the data over four orders of magnitude of serum concentration of metformin and 72 h of data sampling for seven subjects. Also shown as straight lines on the same log-log plots are the GPC's Pareto power-function asymptotes, Eq. (24). A back extrapolated GPC concentration model is shown in Fig. 2, with zero initial concentration that cannot be shown in the Fig. 1 log-log plots.
For semi-log plots please see the .pdf file labelled Supplementary Materials 3. Goodness of fit was shown as relative error plots of residuals of fitting in Fig. 3 for the GPC models and Fig. 4 for biexponential (E2) models. Recall (Appendix ''Concentration data and fitting it'') that relative errors are identically 1=C 2 weighted residuals, and note that the sample-times were approximately equally logarithm of time distributed. Regression in y does not see the x-axis spacing, thus the equal spacing between grouped sample-times on these plots. From Table 1 the GPC models mimicked experimental concentrations and achieved an early concentration peak ranging from 0.66 to 6.17 min total elapsed time. GPC shape parameters from the GD were a\1 for which there were infinite concentrations at t ¼ 0. These GD curves are washout or monotonic decreasing curves with a rate (1/ scale) parameter, b % 1 h À1 . The GPC PD parameters are a and b. The PD shape parameters, a, in the power functions' exponents, were less than one, which obviates calculation of MRT as per Appendix Subsection ''Residence Time''.
The GPC earliest, or start, time is b. The start times, b, of the PD (and thus of GPC distributions), known as the PD scale parameter could not be measured directly (no data at that time) and the unconstrained parameter introduced surfeit variability into the concentration and mass results. By testing b every 5 s from 10 to 35 s and plotting the errors of fit, a good compromise turned out to be an estimated circulatory first arrival time of drug of 25 to 30 s. The GPC harmonic mean residence time (H-MRT) was calculated as in Appendix Subsection ''Residence Time''. Table 1.
Equation (19) gave the statistical context of constant multipliers of the exponential terms of SET functions. For E2 models Table 1 Shown are parameters from gamma densities (GD), Pareto densities (PD) and both from Gamma-Pareto convolution (GPC) fitting of concentrations data for seven dogs In this, any f ðtÞ ¼ ke Àkt term is a pdf, where AUC pdf ¼ 1.
Biexponential fit parameters in this format appear in Table 2.
Fitting for proportional error was performed as per the ''Concentration data and fitting it'' Appendix Subsection. The results of the selection of the GPC and E2 fit functions and their parameters are next characterised by goodness of fit. As per Table 3, the GPC mean relative root mean square (rrms) error was 8.6%, where 10% or less was taken to be good a priori as elsewhere [37,51], and where The blue curves are the GPC concentration models and the green lines are the Pareto (power function) asymptotes of the GPC models. Note the maximum distance between the blue and green curves at 1-2 h. That is the peak effect time of GD convolution. At those times, the convolution of GD function's slow mixing increased the PD function magnitude several fold Fig. 2 This shows the early time pdf, i.e., each having an AUC of 1, of the GPC model (blue) and the gamma density (GD, in red) for dog 5. The GD is shown truncated at its top, and is the faster decaying of two GPC constituent functions. Note that the GPC starts a small time (25 s here) later than the GD, and the GD is less than the GPC for late time Journal of Pharmacokinetics and Pharmacodynamics (2020) 47:19- 45 29 the residuals appear in Fig. 3. The GPC residuals were narrow in range, and approximately normally distributed. Biexponential model (E2) residual errors-see Fig. 4averaged 16.1% rrms or nearly twice that of the GPC fit error, using the same fitting methods. 3 The GPC intravenous bolus arrival time, b, was highly constrained, which adds only a fraction of a degree of freedom, df, to the 4 other parameters. To put it another way, using a 25 to 30 s window for b gives an 8.61% residual error of fitting compared to 8.65% for a fixed 25 s time, which latter has exactly 4 df comparable to the 4 df of an E2, but which had a significantly smaller rrms than E2's 16.1%. The E2 residuals are much more structured than the GPC residuals, see the ''Parameter errors from model-based bootstrap'' Appendix Section for details. The E2 model residuals had an overall 'M' shape, which was also identified for each individual E2 curve. The two 'M' humps are one each from the two local exponential density distributions of a biexponential (a triexponential has three humps). E2 model residuals were non-normally distributed (approximately a three parameter Weibull distribution) and had much wider 95% confidence interval than GPC residuals with an early and late time elongated tail in the underestimating direction (down in the figure). AUC was thus systematically underestimated using E2 models. As implied by the structured residuals, SET models do not systematically approximate the derivatives of concentration. That is, biexponentials fit to this data had slopes roughly analogous to how a broken 24 h clock shows the correct time three times within 72 hsee Fig. 5. This is the effect of not having explicit shape parameters in SET models; the shape of the data, i.e., the derivatives, are only transiently approximated.
Unexplained fraction is 1 À R 2 , a measure of goodness of fit, herein from the multiple correlations of lnðmodelÞ and lnðdataÞ-see Eq. (39). Table 3 shows these unexplained fractions from the E2 models (0.84% mean) to be multiple times greater than from GPC models (0.22% mean) for which the two-tailed Wilcoxon signed-rank of p ¼ 0:0156 suggested significant difference. 4 Clearance In addition to the fit errors, Table 3 shows clearance (CL) values for three models: the GPC, E2 and non-compartmental (NC) models, the latter as per the last paragraph of the ''Concentration data and fitting it'' Appendix Subsection. The NC CL was significantly less than the biexponential CL, but still significantly greater than the GPC CL from Eq. (27) (both p ¼ 0:0078, 1-tailed paired Wilcoxon tests). The NC and GPC correlation was good, r 2 ¼ 0:9840, and the NC CL was 3:6 ml min À1 kg À1 greater than GPC CL. Note that for the GPC model V 0 d ðtÞ, i.e., rate of apparent volume growth with the same units as CL, achieved a limiting value of 6:0 AE 0:6 (mean AE SEM) times the CL.   3 A display in sample groups may seem unusual, and the reader may be more familiar with residuals plotted over time. However, regression in y alone generally injects correlation with x when the x-values are not exactly equidistant, and to be seen, this has to be examined following display of residuals in the same way that regression sees them; as a sequence or list of values (that are then sum-squared and minimised). 4 Statistics: (1) A paired nonparametric test of a correlated sequence of R 2 -values has greater power than comparing two r-values from independent random variables. (2) Without maximum likelihood fitting, and with different residual distributions an AIC (information criterion) comparison of E2 and GPC is intractable.

Half-lives
Half-life defined at the instant in time it is measured was presented as Eq. (6) and in [23]. The fat-tailed GPC distribution, unlike lighter-tailed functions, yielded half-lives of mass and concentration that did not converge to the same value. Instead, they converged to different constants times the elapsed time-Eqs. (31) and (30)-where the ratio of those half-lives converged to a constant as shown in Fig. 7. The time following peak ratio for 95% of the ratio of half-lives of M(t) and C(t) to equal terminal occurred at a mean time of 9.8 h (range 7.4 to 13.1 h), well within the 72 h data collection time. The mean terminal ratio-obtained without urine collection-was 7.0 (range 4.8-9.3)-see Fig. 7 and Eq. (9). The mean terminal ratio of mass to concentration half-lives is 1 for SET functions. Figure 8 shows this effect for our E2 models for which by 21 h 8 min the average t 1=2 ; MðtÞ to t 1=2 ; CðtÞ ratio collapsed to \1:001 and rate of apparent volume growth, i.e., V 0 d ðtÞ, became vanishingly small as the models entered dynamic equilibrium. From Eq. (34) as illustrated in Fig. 7d, the GPC models Àt 1=2 ; V d ðtÞ eventually converged to lnð2Þ t. For an E2 model's t 1=2 ; V d ðtÞ the equivalent asymptotes are where p 2 ¼ 1 À p 1 and where k 1 [ k 2 , which functions reached distinctly different asymptotes that grew much faster (% 10 19 by 72 h) than the single-valued-line GPC asymptote. Thus, the result for the GPC model that Àt 1=2 ; V d ðtÞ $ lnð2Þ t is simpler than the comparable results for ED n ! 2 models.    Figure 9 plots the stable noncompartmental (NC) median values at the geometric mean time for each pair of sample times of the each of 7 dogs' t 1=2 ; CðtÞ for the earliest sample-time pair, then the next earliest and so forth for all 21 pairs of sample-time groups using Eq. (5). This plot is overlaid with the (stable) median values of t 1=2 ; C GPC ðtÞ model curves evaluated at those same geometric mean times. In that figure, the NC model median exhibited a significant trend for increasing values with elapsing time. There is independence between NC half-lives that do not share common time samples, i.e., the odd sample groups, f1; 3; 5; 7; . . .21g, and the even sample groups f2; 4; 6; 8; . . .22g are each independent for trending. Thus, the overall trending of NC half-lives was from the data and not from any assumption made. That is, the NC model median t 1=2 ; CðtÞ trending was due to metformin concentration kinetics. The GPC concentration halflives smoothly followed the trending of the NC median t 1=2 ; CðtÞ results, due to the two shape parameters in the GPC models that allowed reduced-error, derivative matching to serum concentration samples, which in turn allowed for accurate half-life functions to be calculated using Eq. (6). No convincing trend toward a terminal t 1=2value was seen within the data. When expressed 5 parametrically ln t 1=2 ; C NC ðtÞ Â Ã ¼ 1:014 ln t 1=2 ; C GPC ðtÞ Â Ã with 95% CI of 0.966 to 1.062 for slope, and a discarded, not significantly different from 0 intercept, had multiple correlation R 2 ¼ 0:9820 and Pearson r 2 ¼ 0:9823. Lin's concordance correlation [52] was q c ¼ 0:9906, i.e., in the almost perfect range of the McBride scale [53] of q c [ 0:99, with most of the \ 2% error being NC noise. From this, these two very different measurements are almost perfectly measuring the same thing. To be clear, each adjacent sample group pair contained up to 14 samples for extraction of the NC median t 1=2 . This has the same meaning as a plot of the terminal exponentials of NC models fit to data ending at the last time of that sample group pair, and, the more temporal data was included, the longer the t 1=2 was and Fig. 9 shows that t 1=2 ; CðtÞ was a function of time for our data that was also almost perfectly fit by the median GPC model. When this same NC data was parametrically log-log plotted with the median E2 t 1=2 ; CðtÞ, rather than median GPC models, the plot with NC t 1=2 ; CðtÞ was a highly nonlinear sigmoidal function and the concordance dropped to q c ¼ 0:9463, or in the 0.90 to 0.95 moderate concordance range, and with no 95% confidence interval (CI) overlap with the GPC 95% q c CI's. That is, E2 t 1=2 ; CðtÞ did not measure the same thing as t 1=2 ; CðtÞ from GPC and/or sample group median NC t 1=2 .

Continuous infusion simulations
Constant infusion curves with limiting concentration, C SS , were calculated indirectly from their IV bolus models' CDF as implied by Eq. (11). These are shown for the E2 and GPC models for dog 6 as Fig. 10. For dog 6, the GCP infusion simulation of the right panel of Fig. 10 suggests that metformin's persistent volume of distribution growth-see Eq. (4)-appeared to cause a major delay in the time (t ) 72 h) it takes concentration to approach a steady state, C SS , suggesting that measuring metformin CL based on saturation during a constant infusion experiment would be intractable. However, constant infusion is considered by many to be a reference standard for measuring CL. To see why this might be the case, let us assume that the GPC model is the more accurate measurement and that we observe a curve height of C SS ¼ D R CL , Eq. (11). Thus, if one drew a reduced C SS height thinking that the data is from a quickly convergent E2 infusion function, a CLvalue that is 20% too large (0:991=0:826 À 1) would be produced. From Table 3 the NC and E2 CL-values were respectively 25% and 41% greater than the GPC CL-values for dog 6, such that a constant infusion measurement from an underestimated C SS would have appeared to reflect a better standard than other methods that also use exponential approximations to data. However, assigning an underestimated terminal asymptote can be more misleading in another context, i.e., the completeness of collection of drug mass in urine to establish the percent drug recovery. For that problem, estimating percent recovery from the administered dose is numerically much more reliable than from asymptotic completeness based on curve shape.

Simulated dosing regimen for constant drug mass in the body
To illustrate a constant mean body drug burden during any dose interval, s, Eqs. (13) and (14) were applied to dog 6 data, which dog had moderate 72 h drug retention (17.4%- Table 4). The intravenous dose loading needed to maintain a constant average body drug mass in every s ¼ 12 h dose interval is shown in Fig. 11. As the plot is semilog, a constant half-life would be a straight line, whereas in each dose interval, the plot is predominantly log-convex suggesting that the half-life increases for elapsed time during each dosing interval. As time progressed, the peak to trough ratios decreased and half-lives increased as per Eq. (6). Overall, the drug mass peak to trough variation varied less than one-tenth as much as it did for concentration for the same dosing intervals. The serial dose factor multipliers for constructing Fig. 11

Discussion
This work gathers together treatments of models of drug concentration as AUC scaled density functions, variable apparent volumes of distribution equations with half-life defined as a real number variable of time. Also, the density function approach to convolution, a new definition for sums of exponential terms coefficients, an alternative expression for constant (inhomogeneous) volume of distribution from residence time that arises from a new density function model for concentration for power function tailed drugs; the gamma-Pareto type I convolution (GPC) were presented. The gamma-Pareto convolution (GPC) model chains a faster gamma distribution (GD) with a slower type I Pareto distribution (PD), the latter a normalised left truncated power function. As data during the first venous drug arrival times were not available, the b-parametervalues of the PD models were constrained to be 25-30 s and mimic reasonable vascular-drug arrival times. As per Table 1 both a & a\1, and both the GD and PD were logconvex. However, their convolution; the GPC distributions, like the concentration curves they mimic, only became logconvex a short time after the model's peak concentration circa 2.7 min was achieved-see Table 1. The GD models were lighter-tailed than the first compartment E2 k 1 but nonetheless augmented the GPC amplitudes out to relatively late times, circa 20 h-see Fig. 1. GPC models are asymptotic with elapsing time to their second constituent distributions, the Pareto distribution (PD) power functions, Eq. (24). As a\1 for both the PD, and the GPC function asymptotic to it, their mean residence times (MRT) were undefined, and different mean times were found using harmonic MRT (H-MRT). Unlike V d ðtÞ, which grew with time, V H-MRT ¼ CL GPC H-MRT were constant with values \1 (L/kg) that were significantly Spearman rank correlated to E2 V SS values-see Table 1 and the Appendix ''Residence Time'' Subsection. Neither non-equilibrated states nor variable volume modelling prohibit multicompartmental analysis. For example, a two-compartment GPC model would result from assigning a central compartment and a peripheral volume assigned that concentration as a function of time that conserves systemic mass. One might think that without classical multicompartmental modelling of cell membrane transport that pharmacogenomic effects would be undetectable. Dog 1 had the lowest CL and highest 72 h mass eliminations while its terminal rate of volume growth, V 0 d ð1Þ, a parameter only meaningful for fat-tailed models, also had the least value in this short series of seven dogs. The reduced V 0 d ð1Þ-value indicates a relatively reduced rate of tissue absorption of the drug and may be due to a genetically deficient transporter protein. In effect, the requirement for conservation of mass in a single compartment variable volume model allowed physiological parsing into mass loss only to renal clearance and redistribution to be quantified as concentration dilution from apparent volume growth. GPC distributions imply a single fractal recirculatory structure that could be elaborated mathematically. For an introduction to variable volume models and fractals the reader is urged read reference [54]. Although compartments and fractals are of theoretical interest, a model's usefulness arises from the scope of experimental results that it clarifies. Accordingly, this discussion focuses next on what GPC imply models, and how they might be used.

Implications
The usual definition of renal clearance for a non-metabolised renal-only cleared marker is where U is drug concentration in urine, V is urine volume, and P is drug concentration in blood plasma. The assumption is that such a system is in dynamic equilibrium. U V=P for intravenous bolus conditions can be generalised to nonequilibrium states by multiplying it by either side of Eq. (9), i.e., To see this is so, one multiplies Eq. (10) by t 1=2 ; CðtÞ to obtain the fractions of concentration that term-wise parse to mass elimination and volumetric dilution, as follows Then, to correct plasma concentration, P, in the denominator of Eq. (35) to mass elimination conditions, we multiply its numerator by the reciprocals of the mass terms of Eqs. (37a) and (37b), because these reciprocal mass terms are both sides of Eq. (9), and are equal to each other. In other words, for a non-metabolised, renal-excreted drug, the molecules of drug that have shown up in urine, are only those molecules within the plasma that are parsed at that instant of time into urine as opposed to being 'redistributed' by ongoing apparent-volumetric concentration dilution. How long the nonequilibrium state appears to persist depends on the model used. Monoexponential models lack nonequilibrium; they are instantly mixed. As all higher ED n models have an instantly mixed volume, a premature achievement of a terminal volume carries over to biexponential and higher models. In transient nonequilibrium V 0 d ðtÞ eventually disappears-see Fig. 8. To delay the time to dynamic equilibrium one can use NC models or heaviertailed GD models. For example, a population biexponential model of a GFR marker was observed to have a 7.5 times shorter time to 95% of terminal mass half-life than GD models with enforced log-convexity [23]. 6 For a fat-tailed model, the nonequilibrium state is permanent and the ratio of half-lives achieved converged to aþ1 a ¼ 7:0 AE 1:5 (mean AE standard deviation)-Eq. (32) and Table 1. As a (dimensionless) is the shape parameter of the Pareto distribution, and as the Pareto distribution is the GPC terminal tail asymptote, only the plasma concentration tail information would need to be extracted from data to quantify a late-time nonequilibrium U V P experiment. On average CL ¼ 7:0 U V P to within 5% for t [ 9:8 h-see Fig. 7c.
The other shape parameter from the gamma distribution (GD) part of the GPC model, a, aids in shape fitting of the early data. The two shape parameters of the GPC model allow for good derivative fitting between the data and the model; C obs ðtÞ and GPC a b a b t . Even when GD models were used as a stand-alone; without convolution, provided that the 0\a\1 shape had been selected by adaptive regularization, that shape parameter variation allowed matching to the different shapes of the disposition curves for fluid overloaded patients [55]. That E2 models were inefficient for metformin curve fitting was suggested by their poor goodness-of-fit-see Figs 3 and 4-and by their derivatives being unrelated to the data shape-see Fig. 5. From Table 3, the E2 models averaged 16.1% rrms error with 0.84% unexplained fraction, the GPC models' errors were significantly better at 8.6% rrms with 0.22% unexplained fraction. Thus far, there have been interesting implications from GPC modelling, for example, an apparent volume of distribution growth within the body that was terminally 6.0 times the renal clearance, that U V=P needed correction, that terminal half-life of mass appearance in urine was 7.0 times the serum drug concentration half-life value, that the drug mass retained in the body was substantial and that V d and t 1=2 increased with elapsed time during multidosing. Although direct testing of these findings was not within the scope of the current study, some of these results can be contrasted with prior study results, as follows.

Clearance
If metformin is only cleared by renal plasma flow (RPF) extraction, then its clearance should be RPF. We estimated RPF from p-aminohippurate ( The above results suggested that NC curve methods using exponential extrapolation largely prevented underestimation of slope magnitude under the time-samples themselves, but did not prevent a significant underestimation of AUC in the extrapolated tails. This supported the proposition that heavier-tailed GPC models captured more of the otherwise under-extrapolated AUC-see Fig. 3 and Table 4. From a prior work with a GFR marker, backextrapolation using 120-to 240-min data validated with withheld early time-samples was more accurate using a A À B lnðtÞ fit function than using exponential back extrapolation, where the exponential back-extrapolation led to an underestimation of AUC [36]. This agrees with the 6 Adult human expected physical volume of drug distribution from mean residence time; EðVÞ ¼ 15-litres, clearance; CL ¼ 100-ml=min, using a glomerular filtration rate (GFR) marker; 169 Yb-DTPA, from Table 1 of [23]. 7 Geometric mean (GM) is a rarely-reported, better estimator of CLvalues, and GM CL E À Á ¼ GMðCLÞ GMðEÞ . Using GM E PAH of 35 controls (of 65 above) yielded 84.7% (84.8% above) of RPF for GPC metformin. NC overestimation of CL-values seen elsewhere from exponential under-extrapolation of early AUC [60] and of terminal tail areas [5,25,27]. Thus for multiple reasons, the CL-values from the GPC models seemed more plausible than the alternative models applied to the same data.

Metformin drug mass
The survival of drug mass in the dogs' bodies at 72 h in Table 4 of 21.1% remaining represents a 5 fold decrease from the initial dose. By way of comparison, the serum concentration decreased 1000-fold from estimated peak concentration over 72 h-see Fig. 1. As metformin extraction in the dog appeared to be less than estimated RPF, Clearance above, we assumed that drug mass decrease was from renal excretion. Thus, at 72 h on average 78.9% was GPC model renal excreted. Sheen and Sirtori et al. list human urinary metformin recoveries at 60 to 72 h, i.e., 86% and 78.9%-see Table 1 of [19] and Table II of [20]. A 99.9% 48 h urine recovery from IV metformin was calculated by Pentikäinen et al. [21], but may have been inflated from a graphic estimation of the location of maximum cumulative urine mass, thus subject to error as suggested in Fig. 10 and surrounding text. Overall, the drug excretion implied here was similar to published human values, such that we cannot reject the hypothesis that the 72 h GPC model excreted metformin estimates were accurate.

Metformin half-life
As summarised in Table 5  in erythrocytes, a mitochondrial depleted tissue marker, than in plasma [13]. In Table 5, the tissue (or urine) drug mass to plasma (or serum) ratios of half-lives averaged 7:3 AE 1:6 SEM. The GPC model ratios averaged 7:3 AE 0:1 SEM. The average E2 model half-life function ratios in this series at the same times (8; 10; 12; 72 h) would be 3.2, 1.9, 1.3, 1.0, which ratios are not as large as any of the published or GPC half-life ratios. Note that light-tailed modelling does not accommodate mass to concentration half-life ratios persistently greater than 1 until at least 168 h as RBC versus plasma t 1=2 in rats [4].
The GPC model half-life of metformin serum concentration, Eq. (28), is asymptotically proportional to elapsed time, and predicts this trend with approximately twice the half-life observed for the exponential models in this series and in Table 5. This latter is not dissuasive as Fig. 4 clearly shows biexponential underestimation of late time-samples, which is consistent with underestimation of terminal halflives for E2. Examination of the plot of median t 1=2 of grouped time-samples of metformin concentration halflives in Fig. 9 provided strong evidence that longer data collection times yielded longer half-lives.
i.e., old observations only now explained using new modelling concepts.

A dosing strategy for possible use for therapy
In human 850 mg three times daily metformin for 5 days, in comparison to a single dose, Sambol et al. noted drug effect versus no (or lesser) effect for a single dose [24]. Also noted were the effects with elapsed time noted above, increased C max , higher V d consistent with drug mass accumulation, and a prolongation of t 1=2 from 7.2 to 19.8 h. The tentative explanation for this was 'The V d ...is probably larger after multiple doses...' Moreover, the lethargic pace adopted clinically of two weeks before changes are made to metformin dosing, the several months allotted to measuring outcomes in clinical trials [61], the relative sequestration of metformin into erythrocyte tissue versus plasma with elapsing time [13], the decreasing fasting plasma glucose for at least 8 weeks of multidosing [39] are consistent with the expanding apparent volume of drug distribution, increasing half-lives in time, and increasing sequestration of drug in tissue that form the basis for the ''Simulated dosing regimen for constant drug mass in the body''. Although diabetes, as present in Sambol et al.'s patients may reduce glucose transport and retard the time course of drug effects, another hypothesis would be that the majority drug effect is retarded due to slow tissue accumulation in all subjects, including our normal dogs. In any population, a faster initiation of therapeutic cytosol drug content may be available by using the dose loading method hypothesised in this paper, which may find eventual relevance for metformin dosing for cancer chemotherapy and/or other indications.

Synopsis
Sum of exponential terms models have non-zero initial volumes problematic for calculating t 1=2 ; MðtÞ. For example, a monoexponential model has the same concentration half-life as its half-life of drug mass, lnð2Þ k , and as V 0 d ðtÞ ¼ 0, there is no volume growth (or redistribution). Using monoexponentials in 1959, Walser et al. [62] found a multiple day stubbornly persistent delay between plasma disappearance of radiolabeled urea and its appearance in urine of approximately one hour. This implicates ongoing redistribution. It is only hindsight 8 that resolved the paradox in assuming no redistribution and finding irreducible redistribution. This adynamical constraint on apparent volume growth is only transiently mitigated by the inclusion of additional terms in a SET function, e.g., to 21 h for our E2, see Fig. 8.
Herein, fat-tailed metformin models offered for the first time, results consistent with the observed half-life of mass appearance in urine as multiple times the serum drug concentration half-life values with drug mass versus serum drug concentration half-lives which agreed by their ratios and amounts to otherwise irreconcilable results. There were plausible drug masses at 72 h, a plausible explanation for increased V d and t 1=2 with multidosing with a plausible strategy for better multidosing for drug effects, and more plausible clearance-values than from the other methods tested. All of these things suggest potential usefulness for this type of modelling for extracting information and prediction that is problematic using current models for metformin kinetics.
We sought to determine why metformin half-lives did not appear to agree between publications, to develop a hypothesis concerning why metformin control of plasma glucose levels had delayed onset, e.g., following 4-weeks of oral dosing [16], and to speculate concerning the lack of a direct correlation between drug effect and blood metformin concentrations [17]. A single explanation for these phenomena may very well be that the mass half-life is persistently longer than the serum half-life and that mitochondria and cytosol are absorbing much of this mass, although reversibly. All these findings are consistent with high affinity of the drug in or in proximity to the sites of drug action.

Limitations
The Introduction Section listed a number of references to drugs with heavy tails [5,[27][28][29][30][31][32][33][34][63][64][65][66]. However, the potential clinical penetrance for fat-tailed data-spanning models is not known. For example, a power function tail for metformin does not seem to have been previously identified. There was no evidence for nonlinearity for the intravenous route of administration at the drug levels in this study. However, nonlinearity might occur at higher drug levels. A PBPK-type study where blood, urine, and tissues are all collected in an animal study may be appropriate for further investigation. Only a prospective study that measures both drug mass and drug effect can be used to validate how useful the drug mass with drug effect relationship hypothesised for the GPC model is for metformin.
The assay used allowed for 72 h of serum data collection [67]. Due to a technical error, the earliest time-samples were at 20 min [1]. Approximately 20.4% of the AUC occurred before 20-min-see Table 4. This was verified indirectly by comparing the CL to renal plasma flow as estimated from metadata, implying that a lot of this early AUC would be missed using NC methods. Due to this early time back-extrapolation some parameters may be somewhat imprecise but not others, e.g., AUC 20Àmin 0 but not t 1=2 ; Cðt [ 20-min). The last two time-samples were erratic especially for dogs 3, 4 and 5 but not for dogs 1, 2, 6 and 7 as seen in Fig. 1. The apparent increased magnitude and variability of some of the last two time-samples had no clear cause and was not prior identified [1].
The GPC model with simplex method global search optimisation has a long run time, and further refinement of the algorithm for speed would be useful. As a new PK model convergent to a power function tail, the GPC model has only been applied to metformin. How well the GPC model would estimate the plasma/serum concentrations of other intravenously injected substances is unknown at this time. In common with other first applications of a fit function for intravenous drug dosing, there was no modelling of oral dose data, tissue drug concentrations, or correlation with simultaneously measured drug excreted in urine, the applicability is only for linear systems, population kinetics were not explored in the usual sense, and the fitting algorithm used does not allow for case-wise error statistics of parameters, per se, necessitating calculation by other means, for example, as appears in the ''Parameter errors from model-based bootstrap'' Appendix Section.

Conclusion
Fat-tailed models of metformin concentration were fit to data from 7 dogs. The fat tails altered the usual context of some pharmacokinetic parameters. Metformin was cleared at 84.8% of estimated total renal plasma flow. At 72 h there was 21.1% unexcreted metformin. Half-lives eventually became linear functions of the sample-times. The mean ratio of mass to concentration t 1=2 95% converged to 7.0 by 9.8 h. Multidosing for constant drug mass rather than serum drug levels was simulated and may allow for an earlier onset and better maintenance of drug effects but is in need of validation.
Acknowledgements The authors thank Charlotte A. Johnston, DVM; Valerie S. MacDonald Dickinson, DVM and M. Casey Gaunt, DVM as well as the American Journal of Veterinary Research for permission to use data collected, analysed and published in reference [1].
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 https://creativecommons. org/licenses/by/4.0/.

Concentration data and fitting it
The samples ranged from 9.78 to 19; 100 ng=ml metformin. The 20 min sample for dog three was discarded due to quality assurance failure in the original study. In the current study, probable quality assurance failure was additionally identified with squared rank testing in five of 14 of the last two sample-time groups (at 48 and 72 h) leading to discard. Thus, a total of 148 samples of the original 154 (96%) were used for analysis here. Serum samples were analysed using a simple flow injection analysis-tandem mass spectrometry (FIA-MS/MS) by Michel et al. [67], also developed at the University of Saskatchewan. This method has linear response from lowest limit of quantification (LLOQ) of 5 ng=ml to a maximum of 2340 ng=ml metformin base with a relative root mean square of proportional error (rrms) of 5.2% calibration proportional error. The dog doses were weighed and diluted as metformin hydrochloride (165.625 DA average molecular weight) and dose mass corrected to metformin base (129.164 DA average molecular weight) for the computations in this publication.
Noncompartmental analysis of AUC and CL was performed with Prism 5.0 (GraphPad Prism, San Diego, CA, USA) for Johnston et al. [1] and reused here. The AUC values were from the log-trapezoidal rule, with 4-10 early and late time-sample (a generous number) used for extrapolation. This NC technique has the advantage that it is commonly used, but should not be considered optimally comparable to proportional modelling used with better fitting and extrapolating functions than monoexponentials.
This expansion ln f ðxÞ y % ðu À 1Þ ¼ f ðxÞ y À 1, has as its first term the target for the proportional norm minimisation of relative error. Proportional error values closer to 0 are asymptotically ln f ðxÞ y ¼ ln f ðxÞ À lnðyÞ. Thus, log transformations of f ðxÞ and y were used to calculate 1 À R 2 values to make unexplained fractions of total variance indexed to goodness-of-fit.

Relative tail heaviness
The limiting ratio as elapsed time increases of two competing models' survival functions is the statistic of choice for testing relative tail heaviness applicable to density functions, mass functions, functions with discontinuities, functions without derivatives, piecewise defined distributions, and empirical distributions. The use of classification schemes, i.e., indirect comparisons using derivatives and relying upon hazard function theorems, is an unreliable substitute for the use of smoother integral methods for random variable data using survival functions [69]. For continuous functions, pdf binary comparison through survival function ratios avoids false attribution, for example, classifying the GD as having an ED terminal tail, whereas in fact, the exponential has tail heaviness that is within the GD tail heaviness spectrum. 10 This is shown theoretically as follows.
The relative terminal tail heaviness of the S GD ¼ Qða; b tÞ, Eq. (16), and the monoexponential survival function, S ED ¼ e Àk t , is obtained from their ratio which is a 0 0 type indeterminate form, for t ! 1. Applying L'Hôpital's rule yields b k CðaÞ For a; b; k [ 0, the limit as t ! 1 of the right hand fraction of this is dominated by the exponential, which goes to infinity for k [ b and to zero when k\b, which means that the GD tail is only heavier than the ED tail when k [ b. k b does not occur following regression because as b ! k, then a ! 1, and GD becomes statistically indistinguishable from ED. Thus, ED has a GD spectrum tail heaviness, not the obverse. Experimental demonstration of this theoretical result includes the following. k [ b always occurred only when special methods were used for GD fitting. The slower E2 rate constant, k 2 [ b, with b from adaptively regularised GD fitting that minimised AUC error; mean [23]. Herein, the GD tails were lighter even than the E2 fast decays, i.e., all k 1 \b.
Classifications of tail heaviness consist of broad categories of very different functions each having a different spectrum of tail heaviness. Even so, the exponential distribution tail heaviness category did not overlap with the category for Pareto distributions [70]. Even when two functions are in the same category of tail heaviness their range of heavinesses might not overlap, which may seem counterintuitive, but implies once again that only binary tail heaviness comparisons make sense. With regards to PD tail relative tail heaviness, taking the logarithm of the ratio of survival functions simplifies the analysis [69]. The relative terminal tail heaviness of PD versus ED follows from ln S PD ðtÞ S ED ðtÞ This result is unequivocal; the terminal tail of a PD is always heavier than that of an ED.

Residence time
Geometric mean residence time, when robust, can be found using Monte Carlo simulation of a pdf, a tedious process. The median and geometric mean residence times are robust for the GD, PD, exponential and GPC distributions. The harmonic mean residence time (H-MRT) for PD and GPC 9 OLS regression failed Veng Pedersen's sign run criterion for the terminal tail [68] confirming Toutain and Bousquet-Mélou's ( Fig. 17) [46] recommendation to fit data for relative error. 10 Some authors, for lim SðtÞ , where f(t) is a density function and S(t) is its survival function, mistakenly apply L'Hôpital's rule to that limit without having an indeterminate form. The confusion likely arises from the ED's unique property; f ðtÞ  Table 6. Dog 1's V H-MRT is an extreme proportional outlier; its reciprocal being [ 3 IQR (interquartile range) larger than the third quartile V H-MRT reciprocal. However, the V H-MRT -values sort best with E2's V SS (Spearman rank correlation coefficient, rs ¼ 0:96). Indeed, this agreement is nominally better than the agreement between the two E2 apparent volumes, V SS and V d-area (rs ¼ 0:89, apparent steady-state volume of distribution and terminal apparent volume of distribution, respectively). Although the significance of this is unknown, it would appear that if E2 V SS is measuring something, then GPC's V H-MRT is as well.

Noncompartmental time-sample group, median half-life
Given time-samples, ðt i ; C i Þ and ðt iþ1 ; C iþ1 Þ, the half-life of an exponential connecting them occurs once or more during that time interval, and as per Eq. (5) is For subjects having groups of samples drawn at the same times, t iþ1 À t i ¼ Dt i , the numerator of Eq. (5) is a constant. The D lnðCÞ denominator as sampling noise goes to zero is asymptotically a proportional data model-see Eq. (39)-thus is approximately normally distributed. In the limit as Dt ! 0, the signal reduces in the denominator to an asymptotic mean of 0 leaving constant noise with standard deviation, k. Thus, Eq. (5) is a constant times the reciprocal of a normally distributed random variable, which is not the same as the reciprocal of the normal distribution pdf. Rather we substitute 1 x ¼ z in N z ð0; kÞ to yield, which is a special case of the reciprocal normal distribution. For Eq. (5), an attempted mean calculation becomes more erratic when more subjects are included because Eq. (44)'s mean and variance do not exist. This is well known and its tails using ''Relative tail heaviness'' methods are Cauchy distribution similar, i.e., AEx ! 1; S N 1=x ð0; kÞ ðxÞ=S Cauchyð0;kÞ ðxÞ ! 1 . Although the constant multiplier, lnð2ÞDt, of Eq. (5) gets smaller as Dt ! 0, the behaviour of Eq. (44) does not improve (its variance does not exist). N 1=x ð0; kÞ is a bimodal distribution with left and right sided fat-tails for which the median measures centrality as CDF N 1=x ð0; kÞ x ¼ 0 ð Þ¼ 1 2 . Thus, one can use the median values of the subjects' t 1=2 (include negative ones) from Eq. (5) to plot t 1=2 versus sample times. However, the mean is pathologic, i.e., 1 n R n i¼1 X i diverges as n ! 1, and should never be used.
A final consideration is what mean time on the interval t i to t iþ1 of Eq. (5) to assign a half-life as having occurred. As there are 18 to 21 time intervals for the data here, there is only a small numerical range of ½t i ; t iþ1 for assigning the time of occurrence of a half-life. As the logarithm of time approximately linearised the half-life functions examined, the geometric mean time, GMðtÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi t i Â t iþ1 p , was used as the i th time of t 1=2 .

Parameter errors from model-based bootstrap
Nelder-Mead, also called the simplex method, is a fast global optimisation search routine, and is a practical choice for regression of the GPC function to data. The simplex method is not easily fooled by multiple local minima and if searching a space build like an ice cube tray, will find the least error point in that space. However, the simplex method does not calculate the gradients, Hessians or Jacobians generated during Levenberg-Marquardt or other local optimisation routines that rely on gradients to direct the regression to a lesser, but not necessarily least-error, fit location. Those matrices are frequently used for generating and reporting estimated propagated error of case-wise pharmacokinetic model parameters. Another possibility for obtaining case-wise parameters errors is to use model-based bootstrap of regression residuals [71]. The model-based bootstrap procedure randomly resamples residuals with replacement and adds those randomised residuals to the model amplitudes at all the original sample times to create synthetic data and that procedure is then repeated multiple times for each case. The synthetic data are then regressed to yield multiple parameter values for each case from which errors and distribution information not otherwise available are extracted. Because of what is being done, we required normality of residual magnitude. The GPC residual values were approximately normally distributed (Results and Fig. 3). The E2 residuals were non-normal (Results and Fig. 4). The other requirement for model-based bootstrap is that the residuals are relatively unstructured as a list sorted by earliest time of acquisition. We measured the residuals' structure from the data used for Figs. 3 and 4 using analysis of variance (ANOVA). 11 ANOVA of residuals with their sample group average values as the test model yielded an adjusted R 2 adj ¼ 0:1407 for GPC models and R 2 adj ¼ 0:5826 for E2 models. That is, the E2 models had 4.1 times more explained fraction, i.e., structured variance, than GPC models. Thus, model-based bootstrap E2 parameter errors would likely be inaccurate and were not calculated, and for that problem a method that conserves list-wise, i.e., sequence, dependent non-normal loss functions during nonlinear regression would be more appropriate.
Thirty-five GPC simulations were performed and were consistency checked. The bootstrap fit errors ranged from 4.35 to 15.8% rrms% and were not significantly differently distributed from the fit errors from the original data (p ¼ 0:51, Cramér-von Mises, two-sample test), and were approximately log-normally distributed (p ¼ 0:81; n ¼ 35, Cramér-von Mises, one-sample test) with a bootstrap GM 12 of 7.5% and original solution GM of 8.2%, implying that the bootstrap distribution was not significantly inaccurate. The largest value, 15.8%, was a solitary near outlier. 13 This outlier was not implausible for its log-normal distribution (binomial p ¼ 0:47 of at least one value ! 15:8%).
The trial with second largest fit error is illustrated in Fig. 12. This photogenic simulation was selected to illustrate the results of the process of simulated data creation. Visually, the original and simulated data appear to have similarly large noise, and the large noise permits good visual separation of the original data and simulated data.
The 35 model-based GPC simulations, i.e., only 5 per case, were performed as currently limited by the long runtime for each model solution. That is insufficient to establish precise CV% for each case, but as the ratio between gradient based error propagation results and bootstrap is not unusually a factor of two larger or smaller [72], our results were sufficiently precise to roughly gauge CV% magnitude. At the best of times, individual parameters errors do not represent enough information to fully characterise those parameters and for constrained parameters do not contribute useful information. Additional and more precise quality assurance information arose from the 35 bootstrap parameter distribution statistics for all seven cases, which is more relevant for the examination of outliers and robustness of parameters. Table 7 summarises the GPC regression model parameter results from the simulations with all 35 sets of parameter results from bootstrap as well as screening parametric and nonparametric analysis including 95% confidence intervals for mean and median listed in worksheet bootstrap results .xlsx file labelled Supplementary Materials 1.
In this table, the minimum and maximum values correspond to the nonparametric 94.4% confidence intervals for values with n ¼ 35 (Weibull method). The rrms% of case-wise CV% values were used to calculate the global case-wise processing errors of parameters, and the CV% values of all 35 simulations were the population deviations for parameters. Note that the rrms% of individual case  11 ANOVA applied to sample groups provided a population summary of residual variance, rather than case-wise variances as would be a more typical ANOVA application. 12 The GM (geometric mean) of a lognormal distribution parameterised as N lnðxÞ ðl; rÞ is its location, which is also its median value because logarithm is a monotone transform of quantiles. 13 Outliers are usually defined in terms of the first and third quartiles (Q1, Q3) and the interquartile range (IQR ¼ Q3 À Q1). A near outlier is \Q1 À 1:5 IQR or [ Q3 þ 1:5 IQR, but is not so extreme as to be a far outlier. Outliers for a normal distribution occur for erfc 4 erfc À1 1 2 À Á Â Ã of the realisations. Far outliers are \Q1 À 3 IQR or [ Q3 þ 3IQR and occur for erfc 7 erfc À1 1 2 À Á Â Ã , where erfc is the complimentary error function.
parameters CV% were only a fraction of the total CV% for all simulations and all parameters, i.e., the population of dogs had more variability of parameters than the expected case-wise parameter errors. b was a constrained parameter. Its values were located in a 16, 5, and 14 split between, respectively, being at the lower constraint boundary of 25 s, between the 25 and 30 s boundaries, and at the 30 s upper boundary. Even if that fairly even 16/14 split at the constraint boundaries were due to noise, it would be suspiciously unbiased noise. For example, if all the values were grouped at the lower constraint boundary of 25 s, one would have reason to suspect that a better regression answer would occur at some constraint value less than 25 s. The fact that 5 of the 35 b-values occurred between 25 and 30 is not inconsistent with having set plausible constraints, and that despite not having primary study b-values in that time window. This is the physical equivalent of hitting a one-inch wide, elongated strip of masking tape with 5 darts of 35 thrown from 20 feet away. The most variable parameter in the table is the (approximately lognormal distributed) time of peak concentration, which varied from 0.56 to 9.1 min but had only a single (i.e., not implausible) near outlier. That variability did not adversely affect AUCvalues as those latter were of narrow range with no outliers, had low CV% values and insignificantly (p ¼ 0:036; n ¼ 35) negative correlation with peak time.
Both times of peak concentration and b first peripheral venous arrival times are unobserved events predicted by the GPC model whose variability would be better assessed using earliest observations available sooner than the 20 min time-samples in the data set used.
As it goes to robustness, the examination of outliers is important. For a normal distribution near outliers occur for 0.698% of the realisations, whereas far outliers occur for only 0:000234%. 13 Four of the seven (upper) near outliers were from CL, which inherited a heavier than normal right tail from reciprocation because CL ¼ D=AUC, and AUC as well as 1=CL ¼ AUC=D ð Þwere more normally distributed with no near outliers. The other three near outliers were similarly insufficient to suggest any robustness problem, and infrequent near outliers often suggest non-normality, which do not require investigation as erratic values. However, far outliers are so infrequent in the normal case that consideration should be given to pathologic or extreme non-normality of a parameter, to effects arising from some inappropriate quality of the data itself, and/or to some systemic methodological problem. The outliers listed in the table were near outliers, were few in number and were all upper range suggesting heaviness of right tails for some parameters, and no far outliers to suggest a lack of robustness of the methodology were seen.