Predicting topical drug clearance from the skin

For topical drug products that target sites of action in the viable epidermal and/or upper dermal compartment of the skin, the local concentration profiles have proven difficult to quantify because drug clearance from the viable cutaneous tissue is not well characterised. Without such knowledge, of course, it is difficult—if not impossible—to predict a priori whether and over what time frame a topical formulation will permit an effective concentration of drug within the skin ‘compartment’ to be achieved. Here, we test the hypothesis that valuable information about drug disposition, and specifically its clearance, in this experimentally difficult-to-access compartment (at least, in vivo) can be derived from available systemic pharmacokinetic data for drugs administered via transdermal delivery systems. A multiple regression analysis was undertaken to determine the best-fit empirical correlation relating clearance from the skin to known or easily calculable drug properties. It was possible, in this way, to demonstrate a clear relationship between drug clearance from the skin and key physical chemical properties of the drug (molecular weight, log P and topological polar surface area). It was further demonstrated that values predicted by the model correlated well with those derived from in vitro skin experiments. Supplementary Information The online version contains supplementary material available at 10.1007/s13346-020-00864-8.


Introduction
Skin disease affects millions of people worldwide [1]. The effective delivery of a locally acting, dermatological drug demands knowledge of its 'skin pharmacokinetics' to determine the rate and extent with which it reaches its site of action in the epidermis/dermis. Of necessity, this requires understanding of not only the input rate of the drug into the skin but also its clearance from the 'skin compartment' into the systemic circulation. In this work, the phrase 'clearance from the skin' is used to mean the volume of skin from which a drug is completely removed per unit time.
For topical drug products that target sites of action in the viable epidermis and/or upper dermal compartment of the skin, the local concentration profiles have proven difficult to quantify because both drug input into the viable cutaneous tissue and its clearance therefrom are not well characterised [2]. Without such knowledge, of course, it is difficult-if not impossible-to predict a priori whether and over what time frame a topical formulation will permit an effective concentration of drug to be achieved within the skin compartment.
Mathematical and pharmacokinetic modelling has made a substantial contribution to the interpretation of drug movement and disposition in the skin [3,4]. However, given the multistep nature of the dermal absorption process, the many formulation types and the complex nature of the physiological barrier, all models suffer from one or more limitations. For example, the need to 'guesstimate' several parameters to permit simulations to be performed, or the incorrect relation of 'rate constants' to drug physicochemical parameters, or the inability to obtain a prediction of drug concentration in the viable epidermis, means that many models (for most dermal products) are of little practical use.
There is a need, therefore, to develop simple, yet realistic and mechanistically meaningful, models to estimate the key dermato-pharmacokinetic parameters. Among the approaches currently under investigation are physiologically based pharmacokinetic (PBPK) models, which consist of physiologically realistic compartmental structures into which input parameters from different sources (e.g. in silico predictions, in vitro or in vivo experiments) can be combined to predict plasma and/or tissue concentration-time profiles [5,6].
PBPK models take into account physiological properties of the tissue in question, which are not dependent on the drug and can, therefore, be applied to any compound, as well as characteristics intrinsic to the drug. PBPK models employ a 'bottom-up' approach, as opposed to the 'topdown' approach of classical pharmacokinetic models (e.g. one-or two-compartment approaches) [7]. That is, rather than estimating model parameters based on in vivo data (commonly derived from plasma/blood concentration versus time profiles), PBPK model parameters are determined a priori from in vitro experiments, in silico predictions or, if required, in vivo data.
Nonetheless, most PBPK models require a level of parameter calibration and/or optimisation. In general, a drug's concentration in plasma is determined by the systemic volume of distribution at steady-state (V SS , L), defined as the total amount of drug in the body divided by the drug concentration in the plasma [8], and clearance (Cl, L/h), which is the volume of fluid (plasma or blood) cleared of drug per unit time. Assuming a simple one-compartment model with 1st-order elimination from the tissue compartment, the ratio of these independent physiological parameters provides the systemic elimination rate constant k e (Eq. 1): V SS is an apparent volume that describes the extent of drug distribution and binding to the tissues and plasma (Eq. 2): where V plasma is the volume of the plasma and V tissue,i is the volume of the ith tissue; K tissue/plasma is the tissue-toplasma partition coefficient; E i is the tissue extraction ratio and, for non-eliminating tissues, as is generally true of the skin [9]; and E i equals 0.
In this paper, we test the hypothesis that valuable information about drug disposition, and specifically its clearance from the skin, can be derived from available systemic pharmacokinetic data for drugs administered via transdermal delivery systems. When a transdermal patch is applied, the drug delivery rate to the skin has been determined (the 'input rate') and the resulting systemic plasma concentration versus time profile has been measured both during patch wear and after its removal. The decline in the systemic plasma concentration post-patch removal is characterised by a terminal systemic rate constant (k terminal ) that typically is much smaller than the elimination rate constant determined following intravenous administration, demonstrating clearly that 'flip-flop' kinetics are operative [10]. In other words, in transdermal drug delivery, the skin desorption rate is normally much slower than the systemic elimination. As a result, the disposition of a drug following transdermal application is usually rate-controlled by skin desorption. We have therefore assumed that the terminal rate constant post-patch removal reflects the elimination rate constant from the skin (i.e. k e,skin = k terminal ). In addition, we hypothesise as illustrated in Fig. 1 that, in flip-flop conditions, k terminal is related to the drug's clearance from the skin (Cl skin ), via the corresponding 'local' volume of distribution (V SS,skin ), and therefore, k e,skin can be related to key physicochemical parameters of the drug.
To test this hypothesis, the transdermal delivery literature has been searched and values of k terminal identified for the 18 drugs present in 25 FDA-approved products (FDA Orange Book database of the end of 2017) for this route of administration (see Table 1). Information was also included (Table 1) on a lidocaine patch, for which systemic pharmacokinetic data are available, even though this product is administered to elicit a local rather than systemic pharmacological effect. Importantly, the physicochemical properties of these transdermal drugs are quite broad: for example, molecular weights (MW) between 160 and 470 Da and logarithm of the octanol-water partition coefficient (log P) values from ~ 1.0 to 5.0 ( Table 2).  Table 1 Reported values of k terminal , V SS /BW, F i,pH=7.4 , f u,p and predicted values for K (skin/p) and V SS,skin /BW of 19 transdermal drugs a Calculated from literature values of t 1/2 assuming first-order terminal phase kinetics; sources are provided in Table S2 of the Supplementary Information b Deduced from pharmacokinetic data obtained after administration of the drug; sources are provided in Table S1 of the Supplementary Information c Calculated using Eq. 4 for the pK a values listed in Table 2 d Extracted from [11] e Calculated using correlations described in Table 3 f Calculated using Eq. 3 g Also reported as estradiol h Also reported as norethindrone acetate i Also reported as norelgestromin j Also known as nitroglycerin(e)

Drug
Trade name k e,terminal a (h −1 )

Materials and methods
Evaluation of dermal drug clearance (Cl skin ) involved the following steps: (i) identifying the terminal half-life (t 1/2 ) and the corresponding terminal rate constant (k terminal ) from the drug's systemic plasma concentration versus time profile after removal of a transdermal patch, (ii) estimating the drug's volume of distribution in the skin (V SS,skin ) and, finally, (iii) calculating drug clearance from the skin assuming k e,skin = k terminal .

Identifying the terminal half-life (t 1/2 ) and rate constant (k terminal )
A literature search of pharmacokinetic studies performed on 19 transdermally delivered drugs provided information from which the terminal half-lives (t 1/2 ) were derived. The literature search used PubMed and different combinations of the keywords: 'cutaneous', 'skin', 'transdermal patch', 'pharmacokinetic' and 'clearance'. Pharmacokinetic information was also obtained from relevant Drug@FDA public repositories (https:// www. acces sdata. fda. gov/ scrip ts/ cder/ daf/, FDA-Clinical Pharmacology and Biopharmaceutics Review(s),). The data used were from healthy adults (18-71 years; n ≥ 5), and the C p versus time profiles analysed had at least 3 measurements after patch removal. There were, however, two exceptions: (i) as almost all studies involving testosterone were performed on patients with hypogonadism, these data were accepted for analysis; and (ii) similarly, for methylphenidate, the only data available were from children (6-12 years).
When not specifically reported, t 1/2 values were extracted from C p versus time profiles assuming that the terminal phase kinetics were first-order (i.e. t 1/2 = ln(2)/k terminal ). The half-lives were derived from graphs using WebPlotDigitizer software (version 3.10, Ankit Rohatgi; Austin, TX, USA), and k terminal values were deduced (Table 1 and Table S2) (together with the corresponding references) in the Supplementary Information). If data from more than one source were available, an arithmetic mean of the k terminal values was calculated. Equation 2 indicates that each tissue/organ contributes to the overall total volume of distribution (V SS ) of the drug and suggests that V SS,skin can therefore be approximated by the following: where V skin and A skin are respectively the volume and area of the skin compartment, and K (skin/p) is the drug's partition coefficient between the skin and the plasma. It is relevant to point out that protein binding can occur within the skin (in the stratum corneum (SC) and/or the viable skin) and that K (skin/p) is likely to be greater than 1, therefore. Also, as most transdermal drugs are either neutral compounds or weak bases, the degree of ionisation of the latter is an important parameter to consider as well.

Estimating the drug's volume of distribution in the skin
Correlations listed in Table 3 have been developed by Yun and Edginton [8] to predict K (skin/p) in the rat from physicochemical descriptors (log P, degree of ionization (F i ) and plasma protein binding (f u,p )) together with organism-specific information (specifically, the ratio of the systemic volume to body weight (V SS /BW, L kg −1 ) in the rat) for moderate to strong bases (pK a ≥ 7.4, equation A) and for acids and neutral compounds, zwitterions and weak bases (pK a ≤ 7.4, equation B) [8].
To adapt this approach to estimate V SS,skin in humans using Eq. 3, V skin and A skin were assumed to be 2.6 L and 1.8 m 2 (based on a standard body weight of 70 kg, [12]) and K (skin/p) was determined from either equation A or B (Table 3) using the appropriate V SS /BW. We note that the dermis makes up 95-96% of the total skin weight in humans and the epidermis makes up the remainder [12]. As Yun and Edginton [8] did not distinguish between the different skin layers, it was assumed here that K (skin/p) is the average value for all skin layers.
Ideally, in the selection of an appropriate model for determination of K (skin/p) , a determination must be made as to how accurately the selected model reflects the properties of the skin and the plasma. One needs to bear in mind that the correlations developed by Yun and Edginton [8] were built using V SS /BW from rats. In contrast, in this work, K (skin/p) predictions were made using V SS /BW from humans. This approach was taken due to the absence, at least to the authors' knowledge, of a model/correlation for K (skin/p) prediction using human data. Although rat skin is commonly used in in vitro and in vivo percutaneous studies [13,14], its properties do not perfectly mimic human skin [14][15][16][17]. However, in comparison with human skin, rat skin does have a similar stratum corneum (SC) thickness, although a slightly thinner epidermis and total skin thickness [14].
(3) V SS,skin ∕A skin = (V skin ∕A skin ) × K (skin∕p) Values of the steady-state systemic drug volume of distribution per BW in humans (V SS /BW) and the plasma protein binding (f u,p ) were obtained from the literature ( Table 1). The relevant drug physicochemical properties are listed in Table 2. The degree of ionization (F i ) at physiological pH 7.4 (F i,pH=7.4 ) for a chemical was calculated in the normal way using Eq. 4: where g = + 1 for a monoprotic acid and − 1 for monoprotic base. Table 1 lists V SS,skin /A skin as well as the parameter values used in its calculation (i.e. V SS /BW, f u,p , F i,pH = 7.4 and K (skin/p) ).

Calculating drug clearance from the skin
Finally, assuming that skin pharmacokinetics can be described using a one-compartmental model with first-order elimination kinetics, Cl skin /A skin (which has units of cm h −1 , i.e. the same as those for the skin permeability coefficient) was calculated using Eq. 5.

Drug physicochemical parameters
The molecular descriptors are as follows: log P, molecular weight (MW) and melting point (MP) are from EPA's CompTox Chemistry Dashboard (https:// compt ox. epa. gov/ dashb oard), which includes experimental values that were used when available. When several experimental values were available, the mean value was used. In addition, the logarithm of the octanol-water distribution coefficient at pH = 7.4 (log D 7.4 ), the number of rotatable bonds (RotB), the numbers of hydrogen-bond acceptors (HBA) and donors (HBD) and their sum (HBT), molecular volume (MV), topological polar surface area (TPSA) and pK a were calculated using ACD/Labs (Toronto, Canada, version 5.0). The physicochemical properties of the 19 drugs are listed in Table 2.

Multiple linear regression model development
Multiple linear regression (MLR), using the ordinary least squares (OLS) method, was used to develop an empirical Table 3 Correlations for predicting drug skin-to-plasma partition coefficients (K (skin/p) ) in rats. relationship that best described the dependence of the log transformed area normalised skin clearance (log Cl skin /A skin ) on the key physicochemical properties of the drug. Stepwise MLR was performed using The Unscrambler® X software (Version 10.5, Camo A/S, Oslo, Norway). In each regression analysis, a variable was either added or removed until the fit obtained had the highest adjusted and predicted coefficient of determination (R 2 ) and when all the predictors were statistically significant (p value ≤ 0.05). In addition, collinearity between predictor values was assessed by screening the variance inflation factor (VIF) for each equation; a VIF = 5 was the cut-off criterion [18]. Subsequently, internal validations were undertaken to estimate the predictive value of the final model, defined by the determination coefficients of leave-25%-out (Q 2 25% ) and of leave-one-out cross-validation (Q 2 LOO ) [19]. In leave-25%-out cross-validation, the data set for the 19 compounds were randomly divided into training (n = 14) and test (n = 5) groups. Then, coefficients of determination (R 2 ) and prediction (Q 2 25% ) were obtained by regressing the parameters of the model equation to random combinations of the 14 training observations. The calculation of Q 2 LOO involved the omission of the data for one drug and re-determining the regression model using the remaining 18 data. The resulting equation is then used to predict the dermal clearance of the omitted chemical. The correlation between the predicted and observed values in the newly generated dermal clearance data set is used to judge the fit. Q 2 LOO is therefore able to validate the model without additional compounds or splitting the data.

Model verification using in vitro skin permeation
In vitro permeation experiments were performed to measure the rate at which three drugs (buprenorphine (BUP), nicotine (NIC) and diclofenac (DF)) are cleared from the skin following application of examples of commercially available transdermal patches. The transdermal patches tested were the following: Transtec® (35 μg h −1 , 20 mg of buprenorphine over 25 cm 2 ) from NAPP (Cambridge, UK), Nicotinell® (7 mg/24 h, 17.5 mg of nicotine over 10 cm 2 ) from Novartis (Camberley, UK) and Voltaren® (medicated plaster, 180 mg of diclofenac epolamine over 140 cm 2 ) was from GlaxoS-mithKline (Munich, Germany). Experiments were carried out in static, Franz diffusion cells (Permegear, Hellertown, PA, USA) with a receptor volume of ∼ 7.4 mL. Porcine skin (thickness ~ 750 µm) from a single pig, sourced, stored and prepared as previously described [20] was used. A 1.54cm 2 disk of the patch was applied to the skin surface before mounting in the Franz cell; 10 passes of a custom-made, weighted roller ensured complete adhesion between patch and skin. After assembly of the diffusion cell, the lower compartment was filled with a drug-specific receptor solution ( Table 4). The patch was applied for a specific 'uptake time' and then removed; subsequently, the skin remained mounted in the diffusion cell for a further 'clearance time' (or times) ( Table 4), at the end of which the experiment was terminated. Six replicates of all measurements were made. 'Uptake times' were chosen to be sufficiently long that steady-state diffusion had been achieved; 'clearance times' were selected such that an obviously significant reduction in the quantity of drug taken up into the skin had occurred without compromising the ability to detect that remaining.
After dismantling the diffusion cell, the SC was removed by repeated tape-stripping. Templates (Scotch® Book Tape, 3M, St. Paul, MN, USA), with a circular internal area that matched the 1.54-cm 2 patch area, were positioned on the skin, and then a total of 20 adhesive tape strips (2.0 cm × 2.5 cm, Scotch® Book Tape) were sequentially applied, pressed down firmly and quickly removed. Drug was extracted from the individual tapes (Table 4) and the total amount therein was quantified by HPLC (Dionex, UK) ( Table 5). The quantity of drug in the remaining skin tissue was also determined with the same analytical procedure after extraction. The concentrations of the three drugs in the receptor solution (measured at the end of each experiment by HPLC) were always less than one-tenth of the corresponding aqueous solubilities confirming that sink conditions  Table 5.

Results and discussion
The analysis of the literature for 19 different drugs identified more than 70 specific studies that yielded a total of 160 terminal half-life values. In some cases (such as scopolamine and norethisterone acetate), only a single half-life was found while, for other drugs (such as nicotine and ethinyl estradiol), 20 individual half-life values were discovered. The distribution of the half-lives for each drug is presented as a box-and-whisker plot in Fig. 2. From the half-life value, k terminal (= (ln 2)/t 1/2 ) for each drug was then calculated (see 'Supplementary Information', Table S2 together with the corresponding references).
To estimate the area-normalised volumes of distribution in the skin (V SS,skin /A skin ) of the 19 drugs, the skin-to-plasma partition coefficients (K (skin/p) ) were first determined using the equations in Table 3; values ranged from 0.8 for scopolamine to 14.5 for rotigotine ( Table 1). The relatively high K (skin/p) of the latter drug was not unexpected as its high lipophilicity (log P = 4.7) and substantial systemic volume of distribution (V SS /BW = 53.8 L kg −1 ) already suggest that rotigotine is likely to accumulate in tissues. Likewise, scopolamine's smaller lipophilicity (log P ~ 1) and relatively small volume of distribution (V SS /BW = 1.0 L kg −1 ) are completely consistent with its much smaller K (skin/p) . It is worth mentioning that there is only an 18-fold difference in skin partitioning between rotigotine and scopolamine despite the more than 4000-fold difference in their octanol-water partition coefficients. This is consistent with a skin compartment that is dominated by the 'viable' skin layers rather than the SC. The skin-plasma partitioning represents, therefore, the equilibrium between the relatively hydrophilic layers of the skin and hydrophilic plasma.
From the calculated K (skin/p) for each drug and the V skin /A skin , V SS,skin /A skin was then calculated from Eq. 3; the results are listed in Table 1. Finally, drug clearance from the skin (Cl skin /A skin in cm h −1 ) was assessed using Eq. 5 and the resulting values are in Table 6.
An empirical model was then derived, using multiple linear regression (MLR), to predict dermal drug clearance. MLR is a statistical technique that can use a number of molecular descriptors to identify predictive, albeit empirical, relationships in data sets. The advantage of MLR is its simplicity and the easily interpretable mathematical results. The sign of the coefficient derived for each molecular descriptor indicates whether it contributes positively or negatively to the predicted parameter and its magnitude is a measure of the relative importance. However, MLR works best when (i) the structure-activity relationship is linear, (ii) the set of molecular descriptors is independent (i.e. descriptors do not show collinearity) and (iii) the number of compounds in the training set exceeds the number of molecular descriptors by at least a factor of five [22].
At the outset, MLR was performed using ten molecular descriptors (MW, log P, MP, logD 7.4 , RotB, HBA, HBD, HBT, MV and TPSA) as potential predictors of drug clearance from skin; these values for the 19 drugs are listed in Table 2. In the end, a model based only on MW, log P and TPSA (in units of Å 2 ) best explained the calculated dermal clearance (Eq. 6; Table 7; Fig. 3): with an adjusted R 2 = 0.67, Q 2 LOO = 0.61 and p < 0.05 for all three variables; adjusted R 2 is the square of the determination coefficient adjusted for degrees of freedom; Q 2 LOO is the cross-validated (leave-one-out) square of the  0.14 0.12 0.10 determination coefficient; and the p value is related to the significance of the parameters ( Table 7). The general principle of cross-validation is to split data into training and test sets. The former is used to fit the model while the latter serves to evaluate the fitted model's predictive adequacy. Leave-one-out (LOO) cross-validation repeatedly partitions the data set into a training set which consists of all data points except one and then evaluates the predictive density for the held-out data point where predictions are generated based on the leave-one-out posterior distribution. The LOO estimator is nearly unbiased. Moreover, there was no evidence of collinearity in the predictors, with all VIF values being less than 5.
On the other hand, due to the small size of the data set here, the leave-25%-out cross-validation, which is normally employed when the data set is large, is highly variable and depends on which observations are in the training and test sets. To illustrate this point, the leave-25%-out cross-validation was performed 10 times and yielded the results in Table 8. Clearly, regressions using the same data set (but different subsets) produced quite different results (e.g. compare regressions II and X). Therefore, for the relatively small data set involved in this study, the LOO is a preferable and more appropriate method of validation.
In general, drug permeability across biomembranes, including the skin, increases with increasing permeant lipophilicity and decreases with increasing molecular size [23][24][25][26]. It is, therefore perhaps, not surprising that both log P and MW appear in the empirical relation describing drug clearance from the skin. TPSA was also found to be a significant predictor of drug clearance from the skin. TPSA is a molecular property related to the polarity, hydrogen-bonding potential and water solubility of organic molecules [27,28]. TPSA has been shown to be inversely correlated with drug transport across the brain-blood barrier [29,30] and the intestinal membrane [31][32][33], and the positive correlation found here is therefore somewhat surprising. While an increase in hydrogen bonding activity (both acceptor and donor) has been shown to result in a decrease in the partitioning into the organic phase due to the free energy cost associated with the disruption of the hydrogen bonds in the aqueous phase [34], the positive correlation observed with skin clearance may indicate the eventual importance of aqueous solubility in the sequential processes involved in drug elimination from the skin [35,36]. Further speculation on this issue, in the absence of additional data, however, is not warranted at this time.
In an attempt to further validate the predictions of the model, a series of in vitro experiments were performed with three transdermal drug products. The amounts (A, normalised by patch area) of NIC, BUP and DF in the SC (removed by tape stripping), and in the remaining skin, were measured (a) immediately after patch removal (uptake) and, separately, (b) after further periods of time subsequent to patch removal (clearance) and used to calculate the elimination rate constant of the drugs from the skin (k e,skin ). Assuming first-order kinetics, k e,skin was estimated from the slope of the log-transformed mass (in the SC plus in the epidermis/dermis) versus time data  Fig. 3 Relationship between log Cl skin /A skin calculated from experimental data and that predicted by the MLR-derived Eq. 6. The solid line is the line of identity. The number against each point corresponds to that assigned to each of the 19 drugs in Table 6. The blue triangles show the predicted log Cl skin /A skin values compared with those calculated from in vitro experiments (see Table 9 below) for buprenorphine (B), nicotine (N) and diclofenac (D) Table 8 Coefficients of determination (R 2 ) and determination coefficient of leave-25%-out (Q 2 25%) for 10 regressions on the same data set

Regression number
Drugs omitted (see Table 6 for code)   Fig. 4 Amount of drug in the SC + epidermis/dermis (mean ± standard deviation; n = 6) as a function of the time of clearance. The slopes of the linear regressions (dashed lines) provide the drug elimination rate constants from the skin (k e,skin ) Table 9 Measured and deduced dermato-pharmacokinetic parameters for three transdermal drugs, including the skin clearances assessed experimentally in pig skin and empirically predicted (Eq. 6) a Physicochemical properties of diclofenac: MW = 411.32 Da (diclofenac epolamine), log P = 4.16 (diclofenac acid), TPSA = 49.33 Å 2 , F i = 1, f u,p = 0.99, pK a = 4 b Drug amounts (A, mean ± standard deviation, n = 6) recovered after uptake or clearance from the SC, epidermis/dermis, or the sum of both as indicated c k e,skin is the first-order elimination rate constant describing drug clearance from the 'skin compartment' in the in vitro experiments (see Fig. 4) d Values derived from in vivo studies after IV administration of the drug (the values were collected from drug approval packages for the product name listed on the FDA website; see Table S1 in the Supplementary information) e Value derived from in vivo studies after oral administration of diclofenac epolamine [37] f Calculated using Eq. 3 g Calculated using V SS,skin /A skin (from Eq. 3) and the experimentally (in vitro) determined k e,skin h Calculated using V SS,skin /A skin (from Eq. 3) and the reported in vivo k e,terminal (