Application of RP-18 TLC Retention Data to the Prediction of the Transdermal Absorption of Drugs

Several chromatographic parameters (RM0 and S obtained from RP-18 TLC with methanol—pH 7.4 phosphate buffer mobile phases by extrapolation to zero concentration of methanol; Rf and RM obtained from RP-18 TLC with acetonitrile—pH 7.4 phosphate buffer 70:30 v/v as a mobile phase) and calculated molecular descriptors (molecular weight—MW; molar volume—VM; polar surface area—PSA; total count of nitrogen and oxygen atoms—(N+O); H-bond donor count—HD; H-bond acceptor count—HA; distribution coefficient—log D; total energy—ET; binding energy—Eb; hydration energy—Eh; energy of the highest occupied molecular orbital—EHOMO; energy of the lowest unoccupied orbital—ELUMO; electronic energy—Ee; surface area—Sa; octanol-water partition coefficient—log P; dipole moment—DM; refractivity—R, polarizability—α) and their combinations (Rf/PSA, RM/MW, RM/VM) were tested in order to generate useful models of solutes’ skin permeability coefficient log Kp. It was established that neither RM0 nor S obtained in the conditions used in this study is a good predictor of the skin permeability coefficient. The chromatographic parameters Rf and Rf/PSA were also unsuitable for this purpose. A simple and potentially useful, purely computational model based on (N+O), log D and HD as independent variables and accounting for ca. 83% of total variability was obtained. The evaluation of parameters derived from RM (RM, RM/MW, RM/VM) as independent variables in log Kp models proved that RM/VM is the most suitable descriptor belonging to this group. In a search for a reliable log Kp model based on this descriptor two possibilities were considered: a relatively simple model based on 5 independent variables: (N+O), log D, RM/VM, ET and Eh and a more complex one, involving also Eb, MW and PSA.


Introduction
The skin is the heaviest organ in the human body. Its average surface area is ca. 2 m2 and it accounts for about 1/ 10 of the total bodyweight [1]. The skin provides a selective barrier, allowing transdermal delivery of drugs and providing protection against harmful chemicals. Topically applied drugs and other compounds enter the body via either the transepidermal pathway (diffusion across the skin layers) or the appendageal pathway (through hair follicles or sweat ducts), the latter being considered significantly less important [1]. In the transepidermal pathway the molecule permeates the skin either transcellularly or intercellularly and the preferred route depends on the solute's molecular properties-small hydrophilic molecules prefer the transcellular route and lipophilic ones favor the opposite [1]. The rate of drug permeation through skin is expressed as the flux (J)the amount of substance permeated per unit area and unit time. The flux depends on the permeabilityof the skin to the permeant (K p ) and the gradient of permeant concentration across the skin(∆c): J = K p ·∆c.
For passive diffusion, the permeability coefficient K p depends on the partition coefficient P, the diffusion coefficient D and the diffusional pathlength h: K p = P·D/h. The ability of compounds to cross the skin barrier is of great interest to pharmaceutical and medicinal chemists because the transdermal delivery of drugs is an effective alternative to systemic delivery. This ability is also important in the context of environmental toxicology, because many harmful substances enter the body through skin. Transdermal permeation of drugs may be studied using many methods, including in vitro experiments on excised human skin [2] but there are many problems connected with this methodology. Aside from ethical considerations, the compound's skin permeability may differ significantly between individuals, depending on age and race and even between skin samples taken from the same individual. Percutaneous absorption is therefore studied on other, more convenient, models, including pig, rabbit, rat, mouse or shed snake skin, cultured human skin cells, or synthetic membranes [2,3]. However, such models require much tedious experimental work, so alternative solutions have been sought. Skin permeation depends on some readily obtained physicochemical parameters of a molecule, including the well-established predictor of lipophilicity and biological activity of compounds, the partition coefficient between aqueous and organic layers log P ow [4]. However, it was demonstrated that log P ow cannot be applied as a single measure of log K p across a very wide range of chemical families, so molecular weight or volume were incorporated as additional descriptors [5,6]. Further research supported a relationship between K p and hydrogen bond donor and acceptor activity (H d and H a , respectively) [7] and a melting point M Pt [8]. Other authors stressed the importance of solvation free energy [9]. Log K p was also correlated with Abraham's solute descriptors [10][11][12]. Following QSAR studies on skin permeation, models involving log P ow and M w along with further descriptors were developed [13][14][15]. Different computational skin permeability models have been reviewed and compared [15][16][17][18][19][20][21][22], and the most interesting equations discussed there are presented in Table 1.  [10] Where: log P ow -octanol-water partition coefficient; ∆log P oct-hept -the difference between logarithms of octanolwater and heptane-water partition coefficients; M W -molecular weight; V M -molecular volume; M Pt -melting point; SsssCH-sum of E-state indices for all methyl groups; ABSQon-sum of absolute charges on nitrogen and oxygen atoms; H b -total H-bond count; A, B, S, E, V-Abraham's solvation parameters (A-hydrogen bond acidity; B-hydrogen bond basicity; S-dipolar interactions; E-excess molar refractivity; V-McGowan's characteristic volume); Neoplastic-80-antineoplastic-like property at 80% similarity; ALOGP-log P ow calculated using ALOGP algorithm; F06[C-N]-frequency of carbon-nitrogen bond at a topological distance of 06; QXXp-electrostatic interactions between electric quadrupoles of van der Waals forces; π 2 H -solute dipolarity/polarizability; Σα 2 H -solute overall hydrogen-bond acidity; β 2 H -solute overall hydrogen bond basicity.
More recent developments in the field of in silico evaluation of skin permeability include the model proposed by Chenet et al., [24], in which predictions of skin permeability may be based on certain molecular properties and topological descriptors.
The equation proposed by Potts is among the most accurate and the most frequently cited of the computational models of skin permeability coefficient [21,22], although the model developed by Mitragotri [25] is also interesting [21]-supposedly giving even more reliable predictions-but the model by Potts has the benefit of simplicity.
Chromatography is a powerful technique for acquiring measurements of physicochemical and biological properties of solutes, including the ability of compounds to cross different biological barriers [26][27][28]. It is superior to techniques based on excised human skin, animal experiments or even cell cultures because it is more reproducible and it usually involves commercially available chromatographic supports. The chromatographic techniques used to predict the skin permeability of solutes include normal-and reversedphase thin layer chromatography [29,30], immobilized artificial membrane (IAM) column chromatography [31][32][33][34], RP-18 column chromatography [33][34][35], column chromatography on a novel stationary phase based on immobilized keratin [36], biopartitioning micellar chromatography (BMC) [37][38][39], micellar electrokinetic chromatography [34], and liposome electrokinetic chromatography [40]. Available models of skin permeability parameters, based on liquid chromatographic descriptors, are presented in Table 2. Where: K p -skin permeability coefficient; R f -retention factor in thin layer chromatography; R M = log (1/R f − 1) and R M 0 are obtained by extrapolation of R M values from a series of TLC experiments conducted for different concentrations ϕ of a modifier in mobile phases to zero concentration of the modifier, according to the equation: R M = R M 0 + S ϕ; P mw -micelle-water partitioning coefficient; k-retention factor in column liquid chromatography; k IAM , k BMC , k KERATIN -retention factors in IAM, BMC and immobilized keratin chromatography, respectively; k w IAM -k IAM measured at or extrapolated to purely aqueous conditions; ∆log k w IAM -the difference between log k w IAM measured and predicted on the basis of log P ow ; V-McGowan's characteristic volume; M w -molecular weight; M Pt -melting point; K sc -human skin-water partition coefficient.
From the equations presented in Table 2, it can be concluded that skin permeability coefficient is connected with chromatographic parameters (log k or R M 0 for column and thin layer chromatography, respectively) via linear or reverse parabolic relationships. Chromatographic retention parameters are used either as sole skin permeability predictors, or they are combined with additional descriptors (log P ow , V, M w or M Pt ).
The objective of this study was to examine the relationships between the skin permeability coefficient log K p and different, RP-18 chromatography-derived descriptors and other physicochemical parameters for a large group of structurally unrelated compounds, mainly drugs.

General Considerations
The majority of compounds analyzed in this study are drugs, administered either orally or in the form of injections, and whose transdermal delivery, in certain cases, is of interest. On the other hand, some topically applied substances (e.g., cosmetic preservatives) may exhibit the unwanted ability to cross the skin barrier. The skin permeability coefficient (K p ) is an important parameter that helps in the assessment of compounds' epidermal permeability; however, the experimentally determined values of K p are available for only a few drugs, within the studied group. For this reason it was decided that the models of skin permeability involving thin layer chromatographic and calculated descriptors should be generated and validated using K p values obtained in silico with the EpiSuite software (DERMWIN v. 2 module) (log K p EPI ), recommended by the US Environmental Protection Agency [41] and, at a final stage of these investigations, tested on a sub-group of analyzed solutes whose experimental log K p values could be found (log K p exp ). The estimation methodology used by DERMWIN is based on the Equation (2), related to the Potts model: The values of log K p EPI obtained using DERMWIN are given in Table 3.  [37] Where: log K p exp -experimental values; log K p EPI -values calculated using DERMWIN software; log K p (6) to log K p (10) and log K p (12)values calculated according to Equations (6)-(10) and (12).

Thin Layer Chromatographic Parameters-Extrapolation Methodology
Reversed-phase thin layer chromatography has been used to predict physicochemical properties and bioavailability of compounds for many years [26] and the chromatographic parameter considered in these investigations is usually the R M value defined by Bate-Smith and Westall (Equation (2)) [42]: The partitioning between chromatographic supports and an aqueous mobile phase resembles that between biomembranes and the aqueous phase. The most common approach to obtain the chromatographic retention parameters for water as a mobile phase is by using a series of chromatographic experiments with mobile phases containing different concentrations ϕ of a water-miscible solvent (organic modifier). Plots of R M vs. ϕ are extrapolated to zero concentration of the modifier to furnish R M 0 and the most common method to do so is by using the linear Soczewiński-Wachmeister Equation (3) [43].
Apart from the R M 0 value, other useful chromatographic descriptors derived from the linear extrapolation of R M to purely aqueous conditions are the slope S and C 0 = −R M 0 /S.
The compounds 1 to 22 and 62 were chromatographed on the RP-18 stationary phase using methanol-pH 7.4 buffer mobile phases as described in Section 3. The values of R M were calculated and plotted against the organic modifier concentration ϕ. The chromatographic parameters R M 0 and S obtained according to Equation (3) ( Table 4) were correlated with the log K p EPI values presented in Table 3, and the resulting correlations were found to be unsatisfactory. Neither R M 0 nor S correlated with log K p EPI for the whole group of 23 compounds. 14 compounds (2, 3, 5, 6, 8, 9, 10, 11, 16, 17, 18, 20, 21, 62) gave reverse parabolic relationships with R M 0 and S ( Figure 1), but the remaining solutes did not fit any reasonable pattern.   From this, it was evident that RM 0 and S obtained in the chromatographic conditions described in Section 3 could not be used as sole predictors of log Kp. Attempts to generate a multivariate relationship between log Kp and the calculated physicochemical parameters, involving RM 0 or S, failed and a purely computational model (4) was obtained by forward stepwise regression: From this, it was evident that R M 0 and S obtained in the chromatographic conditions described in Section 3 could not be used as sole predictors of log K p . Attempts to generate a multivariate relationship between log K p and the calculated physicochemical parameters, involving R M 0 or S, failed and a purely computational model (4) The same set of dependent variables was applied to compounds 1 to 60 to furnish the Equation (5)   The Equations (4) and (5), based on the set of purely computational variables PSA, log D, E HOMO , E LUMO , E h and α, are relatively simple and logical, since they involve some properties responsible for drug absorption (lipophilicity, polar surface area and polarizability) [44,45]. However, apart from poor statistics, the differences between the parameters of the Equations (4) and (5) suggest that they lack the required universality with respect to larger groups of structurally unrelated compounds.

Single Chromatographic Run Retention Parameters-R f
The extrapolation method, although commonly used and recognized, has certain drawbacks. Several chromatographic experiments are required and the extrapolated R M 0 values depend on an organic modifier and its concentration range used to generate R M = f (ϕ) plots. Some studies, therefore, use the single chromatographic run approach in which the R f and R M values are collected using a single concentration of an organic modifier in a mobile phase [46]. The single chromatographic run approach was used to predict the lipophilicity of selected cosmetic raw materials [47] and, in separate work, yielded two chromatographic descriptors, R f and R f /PSA, used to study the blood-brain barrier permeation of solutes [48][49][50]. Apart from providing reliable lipophilicity descriptors, this methodology has additional advantages: it requires considerably less experimental work than any extrapolation or interpolation method, and none of the additional factors mentioned earlier (such as the modifier concentration range or the type of extrapolation curve) need to be taken into consideration.
Following the unsuccessful attempts to generate useful and universal models of the skin permeability coefficient using thin layer chromatographic retention parameters obtained by extrapolation (R M 0 , S), attention turned to the single RP-18 chromatographic run approach developed in the earlier study of blood-brain barrier permeation [48][49][50]. The models presented in those studies involve R f (collected as described in Section 3) and a novel parameter derived from it: R f /PSA. The latter (combined with R f and a number of calculated physicochemical parameters) is a particularly good predictor of the blood-brain barrier permeability (expressed as steady state blood-brain partition ratio, log BB). In the case of the skin permeability coefficient, however, this parameter (and R f itself) failed to be selected in the course of the forward stepwise regression. The most interesting model (6) (Figure 3) generated at this stage of our investigations did not contain any chromatographic parameters and was as follows: log K p = −1.  This model accounts for over 83% of total variability and contains the independent variables introduced in the following order: (N+O), log D, HD. The variables in the equation are strongly related to the conditions of good oral availability [44,45] and the ability to cross the blood-brain barrier [48][49][50][51]. The equation obtained at the very first stage of this regression, containing (N+O) as a sole independent variable, accounts for as much as 70% of total variability. The coefficients for (N+O) and HD are negative, which (as already observed, e.g., by Lien and Gaot [23]) suggests that excessive hydrogen bonding is an obstacle to epidermal permeability.
At this point the group of 60 studied compounds was divided into two subsets: a training set (1 to 40) and a test set (41 to 60).Equation (7) generated for the training set, and containing the same independent variables as Equation (6) The values of log K p were calculated for compounds 41 to 60 according to Equation (7) and plotted against the reference log K p EPI values presented in Table 3. The linear relationship between these two groups of log K p values improved when salicylic acid (60), whose K p EPI value seems to be overestimated compared to K p exp (due to the molecule's combined polarity and acidic properties, 60 is not a very good skin permeant), was removed as an outlier (R 2 = 0.83 and 0.89, respectively).
The model (6) was also tested on a subgroup of 9 compounds, analyzed in this study, whose log K p exp values were available (12, 16, 21, 40, 56, 60, 61, 62, 63). The resulting dependence is linear, with R 2 = 0.92, which is a much better result than that obtained for the relationship between log K p exp and log K p EPI (for n = 9, R 2 = 0.40).

Single Chromatographic Run Retention Parameters-R M
With the promising equations (6) and (7)   In Equation (8) the independent variables were introduced in the following order: The difference between this model and Equation (6) was not very significant; HD was replaced with R M /V M , and two more independent variables, responsible for a very small improvement in statistics, were introduced (E T , E h ), but (N+O) and log D still accounted for the greatest percentage of variability.
The group of 60 studied compounds was, as before, divided into a training set (1 to 40) and a test set (41 to 60). Equation (9), generated for the training set and containing the same independent variables as Equation (8) Values of log K p were calculated for compounds 41 to 60 according to Equation (9) and were plotted against the reference log K p EPI values. Just as in the case of Equation (7), the linear relationship between these two groups of log K p values improved when salicylic acid (60) was removed as an outlier (R 2 = 0.87 and 0.92, respectively).
The model (8) was also tested on a subgroup of nine compounds, analyzed in this study, whose log K p exp values were available (12, 16, 21, 40, 56, 60, 61, 62, 63). The resulting relationship was linear, with R 2 = 0.80, which is a much better result than that obtained for the relationship between log K p exp and log K p EPI (for n = 9, R 2 = 0.40), but not as good as in the case of the Equation (6).
Equations (8) and (9) account for 87 and 82% of total variability, respectively, and have the advantage of being very simple; however, an attempt was made to improve their predictive abilities by adding further independent variables. Equation (10) (Figure 5  The same set of variables as in Equation (10) was applied to the training set (compounds 1 to 40). The Equation (11) Equation (11) improved significantly after one compound (19) was removed. The difference between the K p (11) and K p EPI values for this compound is probably due to its very high energies and PSA, coincidentally leading to overestimated skin permeability, calculated according to Equation (11) The values of log K p were calculated for compounds 41 to 60 according to Equation (12) and plotted against the reference log K p EPI values. Just as in the case of Equation (9), the linear relationship between these two groups of log K p values improved when salicylic acid (60) was removed as an outlier (R 2 = 0.831 and 0.886, respectively).
The model (10) was also tested on a subgroup of nine compounds, analyzed in this study whose log K p exp values were available (12, 16, 21, 40, 56, 60, 61, 62, 63). The resulting dependence was linear, with R 2 = 0.83, which is a much better result than that obtained for the relationship between log K p exp and log K p EPI .

Chemicals
The 63 drugs analyzed during these investigations were isolated from pharmaceutical preparations, purchased from Sigma-Aldrich (St. Louis, MO, USA)or donated as free samples by Polfa-Pabianice (Pabianice, Poland). The purity of drugs isolated from pharmaceutical preparations was assessed by TLC and densitometry (Section 3). All isolated drugs gave single chromatographic spots (densitometric peaks) and were used without further purification. Drugs purchased from Sigma-Aldrich were of analytical or pharmacopeial grade. Distilled water used for chromatography was from an in-house distillation apparatus. Analytical grade acetonitrile and methanol were from Avantor Performance Materials Poland S.A. (formerly POCh S.A., Gliwice, Poland). pH 7.4 phosphate buffered saline was from Sigma-Aldrich.

Thin Layer Chromatography
Thin layer chromatography was performed on 10 × 20 cm glass-backed RP-18 F 254s TLC plates (layer thickness 0.25 mm) from Merck (Darmstadt, Germany). Before use, the plates were pre-washed with methanol-dichloromethane 1:1 (v/v) and dried overnight in ambient conditions. Solutions of compounds 1 to 63 in methanol (1 µg·µL −1 , spotting volume 1 µL), were spotted with a Hamilton microsyringe, 15 mm from the plate bottom edge, starting 10 mm from the plate edge, at 8 mm intervals. The chromatographic plates were developed in a vertical chromatographic chamber lined with filter paper and previously saturated with the mobile phase vapor for 20 min. The mobile phase consisted of acetonitrile-pH 7.4 phosphate buffered saline 70:30 (v/v) or methanol-pH 7.4 phosphate buffered saline (methanol contents from 90 to 50% v/v in 5% increments). The development distance was 95 mm from the plate bottom edge. After development, the plates were dried at room temperature and examined under UV light (254nm) and with a CD60 densitometer (Desaga, Germany, multiwavelength scan, 200-300 nm at 20 nm intervals). All chromatograms were repeated in duplicate, and the mean R f values were used in further investigations. The chromatographic data are presented in Table 4.

Calculated Molecular Descriptors
The molecular descriptors for compounds investigated during this study were calcu-  Table 5. Statistical analysis was done using Statistica v.13 or StatistiXL v. 2.

Conclusions
This work has established that R M /V M is a useful descriptor of skin permeability derived from RP-18 thin layer chromatography. In a search for reliable log K p models based on this descriptor two possibilities were considered: a relatively simple model based on 5 independent variables: (N+O), log D, R M /V M , E T and E h (Equations (8) and (9)) and a more complex one, containing also E b , M W and PSA (Equations (10)- (12)). The latter accounts for over 90% of total variability and involves all the major properties that determine the drugs' ability to cross biological barriers (lipophilicity, molecular size, ability to form hydrogen bonds). It should be mentioned here that a very simple parameter (N+O) accounts for as much as 70% of the total log K p variability; log D and R M /V M account for a further 10% and 5%, respectively.
Skin permeability is a difficult property to measure.Due to the limited availability of experimental permeability data for the solutes investigated in this study, the reference skin permeability coefficients were calculated according to a widely accepted model based on log P ow and M W (log K p EPI ). The advantages of this model are clear-it is based on easily obtained molecular properties, whose influence upon the skin permeability is well documented [21]. Of course, the reference model has also its limitations: it overestimates the results for very hydrophilic molecules [52,53], underestimates the values for non-hydrogen bonding solutes [52], and fails for extremely lipophilic compounds or solutes having a very high tendency to hydrogen bonding [18,53]. However, the group of solutes examined in this study does not include molecules of very high lipophilicity or with a very high tendency to H-bonding, and the differences between log K p EPI and log K p exp for hydrophilic solutes or non-H-bond donors are moderate. Small discrepancies between the calculated and experimental reference values (log K p EPI and log K p exp , respectively) for hydrophilic solutes or those without the H-donor sites may be a likely reason why the correlations between the experimental log K p exp values (where available) and the log K p values calculated according to the equations developed in this study ((6), (8) and (10)) are better than those between log K p exp and log K p EPI for the same subgroup of compounds. To conclude, Equations (6), (8) and (10) have been shown to be efficient tools for skin permeability predictions. Although Equation (10) provides the closest correlation, Equations (6) and (8)

Data Availability Statement:
The data presented in this study are available in this manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. H a -H-bond acceptor activity; IAM-immobilized artificial membrane; BMC-biomedical chromatography; k-retention factor in column liquid chromatography; R f -retention factor in TLC; TLC-thin layer chromatography; RP-reversed phase (chromatography); P mw -micelle-water partitioning coefficient; k IAM and k BMC -retention factors in IAM and BMC chromatography, respectively; log k w IAM -log k IAM measured at or extrapolated to purely aqueous conditions; ∆log k w IAM -the difference between log k w IAM measured and predicted on the basis of log P ow ; V-McGowan's characteristic volume; K sc -human skin-water partition coefficient; DM-total dipole moment; V M -van der Waals molar volume; S a(a) -surface area (approximate); S a(g) -surface area (grid); E HOMOenergy of the highest occupied molecular; E LUMO -energy of the lowest unoccupied molecular orbital; E T -total energy; E b -binding energy; E e -electronic energy; E h -hydration energy; Rrefractivity; α-polarizability; HD-H-bond donor count; HA-H-bond acceptor count; (N+O)-total nitrogen and oxygen atom count; Neoplastic-80-antineoplastic-like property at 80% similarity; ALOGP-log P ow calculated using ALOGP algorithm; F06[C-N]-frequency of carbon-nitrogen bond at a topological distance of 06; QXXp-electrostatic interactions between electric quadrupoles of van der Waals forces; π 2 H -solute dipolarity/polarizability; Σα 2 H -solute overall hydrogen-bond acidity; β 2 H -solute overall hydrogen bond basicity.