PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Mol Pharm. Author manuscript; available in PMC 2010 June 24.
Published in final edited form as:
PMCID: PMC2891763
NIHMSID: NIHMS209577

Computational Model for Predicting Chemical Substituent Effects on Passive Drug Permeability across Parallel Artificial Membranes

Abstract

Drug permeability is often a limiting step in drug action, requiring chemical optimization of a drug candidate to improve this property. Such optimization is typically performed in the context of a congeneric series, where substituents are varied to optimize the target property. Motivated by this need the present work examines the influence of chemical substituents on passive permeability (log P pass) across parallel artificial membranes (PAMPA) undertaken for three congeneric series of compounds; benzoic acids, pyridines and quinolines. PAMPA showed pyridine and quinoline to have high permeability and chemical substituents to typically reduce the permeability. On the contrary, benzoic acid showed poor permeability and chemical substituents typically increased the permeability. To quantitate these effects with respect to physical properties, models were built to explain and predict the permeability of these classes of compounds based on computed molecular descriptors. Models for the benzoic acid series in the ionized state indicated the solvent accessible surface area, cavity dispersion and the free energy of solvation in hexane as well as in water to dominate permeability. However, when the acid group is treated as neutral, the free energy of solvation in water, the fraction polar surface area, the polar surface area and difference in the free energy of solvation between hexane and water were important; these terms, among others, were also important for the neutral pyridine−quinoline series. Considering that the permeability of the benzoic acid series is about 2 orders of magnitude lower than the pyridines and quinolines and that a shift of approximately two pH units in the pK a of the acid group of benzoic acid will allow for the neutral species of the molecule to dominate under experimental conditions (pH = 6.5), these results suggest that the additional energy barrier associated with permeation of the benzoic acid series is associated with the need to protonate the acidic group, thereby forming the neutral species which may then cross the hydrophobic region of the membrane.

Keywords: PAMPA, oral absorption, bioavailability, permeability, quantum mechanics, QSAR, multiple regression, substituent effects

Introduction

Poor bioavailability reduces the potential of a drug candidate being suitable for further development. Oral bioavailability of a drug candidate is affected by a number of processes such as formulation, stability in gastro-intestinal tract, aqueous solubility, passive permeability (P pass), and interactions with various influx and efflux transporters present in gastrointestinal tract,1 among others. Hence, it is extremely difficult to determine the contribution of each of these individual processes on oral bioavailability; however, P pass plays a major role in determining oral bioavailability for most drug-like molecules.24 Therefore, additional efforts to develop approaches to rationally predict and optimize drug permeability are warranted.

Computational techniques can be used as a useful tool to rationally predict the appropriate chemical features to enhance oral absorption by increasing P pass of novel drug candidates.210 A common approach is to study a series of structurally similar drug-like molecules (i.e., a congeneric series) to determine the effect of various chemical groups on the permeability.8, 11 In this approach a congeneric series of pharmaceutically relevant compounds typically are studied and a computational model is designed to explain the effect of various chemical substituents on the absorption and/or permeability of the parent compound by analyzing the relevant molecular descriptors.9 The computational model is then used to predict suitable chemical features that enhance the permeability for that class of compounds. Benzoic acid, pyridine and quinoline (Figure 1) are very common molecular scaffolds found in many pharmaceutical compounds.1214 Yet there is a lack of knowledge to explain the effect of various chemical substituents on the permeability of these molecular scaffolds.9 This omission motivated the present study to elucidate the various substituent effects on congeneric series of benzoic acids, pyridines and quinolines. Models were developed based on experimental data obtained with parallel artificial membranes (PAMPA) along with theoretical data from both empirical and quantum mechanical computational methods.

Figure 1
Structure of (left to right) benzoic acid, pyridine and quinoline

Methods

Permeability Measurements

(a) Artificial membrane preparation and flux across PAMPA

Artificial membranes were prepared daily prior to each experiment. Parallel artificial membranes were created as described elsewhere.15 Briefly, a 96-well microtiter plate (Dynex N4500411, Chantilly, VA) served as the receiving chamber, and a 96-well microfiltration plate (Millipore/Immobilon MAIPN4510, Billerica, MA) served as the donor chamber. Artificial membrane solution (AMS) consisted of 5% (w/v) of crude egg lecithin based phosphatidyl choline (approximately 60%) (Sigma P5394, St. Louis, MO) dissolved in n-dodecane. Transport solution (pH = 6.5) consisted of 50mM sodium phosphate with 3% dimethyl sulfoxide (DMSO) (v/v) and was used to fill each receiver compartment (385µL/well). The microfiltration plate was placed directly on top of the microtiter plate. Each of the test filters was impregnated with 5µL of AMS; reference filters without lipid were impregnated with 5µL of n-dodecane. Bilayer formation in the filter pore channels was promoted by lipid contact with water interface from receiver compartment.16 Transport studies were initiated by adding 350µL of 500µM donor solution to each of the top wells. Donor solution was transport buffer with dissolved pyridine, quinoline, or benzoic acid analogs.

Transport of each pyridine, quinoline, and benzoic acid across PAMPA membranes was performed using n=6. Transport was allowed to occur for a total of 6 hours. At 6 hr, the two plates were separated and donor and receiver samples were sampled for subsequent analysis. In parallel studies, the permeabilities of reference compounds 14C-mannitol and 3H-metoprolol were measured (n=6) within each 96 well plate. 14C-mannitol served to assess membrane integrity and 3H-metoprolol served as a reference for moderately high permeability. 14C-mannitol and 3H-metoprolol assessed plate-to-plate variation across a broad permeability range.

(b) Analytical

Pyridine, quinoline, and benzoic acid samples were analyzed by high-pressure liquid chromatography (HPLC) using a Hitachi L-7000 HPLC module, Hitachi L-7200 auto sampler, and Hitachi L-7400 UV detector. Chromatograms were collected and integrated using Hitachi D7000 Chromatography Data Station Software (ver 3.1). Briefly, 100 µL injections were made onto a reverse phase C18 column (Suppelco Discovery Column 37170-01, Bellefonte, PA) at a 1.1 mL/min flow rate. Mobile phase for the benzoic acid series consisted of 70:30 10 mM potassium phosphate buffer (pH=2.5):acetonitrile. Mobile phase for the pyridine and quinoline series was 65:35 10 mM potassium phosphate buffer (pH=7.0):acetonitrile. The wavelength of absorption used for the benzoic acid, pyridine, and quinoline series was 240 nm, 254 nm, and 254 nm, respectively.

14C-mannitol and 3H-metoprolol were quantified by adding 300 µL of sample to 3 mL of scintillation cocktail (Econosafe RPI111174, Mount Prospect, IL). The samples were subjected to liquid scintillation counting using a Beckman LS-5801 Scintillation Counter at 1 minute scans per sample.

Permeability data was calculated using:

P=VrdCrdtACd
(1)

where Vr is the volume of the receiver compartment, dCr is the change in concentration of the receiver compartment, dt is the change in time, A is the area of the filter membrane corrected for porosity, and Cd is the change in concentration of the donor compartment.

Computational Methods

Congeneric series of benzoic acids, pyridines and quinolines were modeled by varying the chemical substituent at 2−, 3− and 4− positions on the aromatic ring structures (Figure 2). These compounds were built using the program CHARMM17 with the all-atom CHARMM general force field (CGenFF, A.D.M., Work in progress).

Figure 2
Different configuration of benzoic acid based on the orientation of carboxyl group with respect to the substituent (R2, R3 or R4).

Each molecule was subjected to 1000 steps of steepest descent (SD) and 500 steps of Newton−Raphson (NRAP) energy minimization in the gas phase to a gradient of 10−4 kcal/mol Å. CHARMM-minimized structures were then subjected to quantum mechanical (QM) geometry optimization in the gas phase at the RIMP2/6–31G* level. QM calculations were performed with the ab initio package Q-CHEM.18 QM optimized structures were then used to calculate a variety of molecular descriptors for model development. The free energy of solvation of the compounds in water (ΔG-wat) and n-hexane (ΔG-hex) were calculated using the semiempirical quantum chemical program AMSOL6.8 by invoking the AM1-SM5.42R continuum solvent model,19 as this model was developed to avoid the need for further geometry optimization in the presence of the solvent term. The free energy of solvation of each compound in a solvent medium was calculated by taking the difference between the heat of formation plus solvation energy of the solvated system in its relaxed electronic state and that for the gas-phase system. Free energy of transfer from water to n-hexane (ΔG-h/w = ΔG-hex – ΔG-wat) was also calculated. The same program was also used to evaluate the free energy for creating a solute cavity inside the solvent combined with solute-solvent dispersion interactions (CDS) in water, and n-hexane. Solvent accessible surface area (SASA), polar surface area (PSA) and fractional polar surface are (fPSA) were calculated using the Lee and Richards method20 as implemented in CHARMM. The radius of the solvent molecule was taken as 1.4 Å, which approximates the radius of a water molecule.21 PSA for each compound was calculated by adding the contribution of polar atoms (N, O, F, Cl, Br and hydrogen atoms covalently attached to any of these atoms). For calculation of PSA halogens were included as they could potentially take part in polar interaction with other polar atoms.2224 fPSA was calculated by taking the ratio of accessible solvent area formed by the polar atoms to the total solvent accessible area. Other descriptors including molecular weight (MW), octanol/water partition coefficient (logP(o/w)), hydrogen bond acceptor (H_acc) and hydrogen bond donor (H_don) were calculated using the MOE program.25 Octanol/water distributioin coefficient (logD(o/w)) was estimated at experimental pH (6.5) using the following equations, as discussed by Scherrer and Howard.26

logDacid=logP+log(11+10pHpKa)logDbase=logP+log(11+10pKapH)
(2)

pKa values for the compounds were taken from Scifinder Scholar database.27 Calculated results for the above descriptors for all species studied are presented in Table S1 to Table S3 of the Supplementary Information.

Model Development

Single variable linear regression analysis was employed between the log P pass and individual molecular descriptors for each series of congeneric molecules separately. Pearson correlation coefficient (r 2) and p-values were obtained for each parameter to estimate the significance of that descriptor on P pass. Microsoft Excel 2002 was used for regression analysis. The congeneric series of pyridines and quinolines were combined together to form a single group while benzoic series was analyzed alone; efforts to group the three series of compounds and perform correlation analysis were unsuccessful as judged by individual r 2 values of 0.42 or less. Next the molecular descriptors with r 2 > 0.01 were subjected to multivariable regression analysis separately. First, molecular descriptors from each series were taken in all possible combinations of groups of two and all groups were subjected to multivariable regression analysis. No combination of descriptors having internal correlation greater than 0.8 was considered for multivariable regression analysis as determined from the correlation matrix of the descriptors. The correlation matrices for different series of congeners are presented in Table S4 to Table S6 of the Supplementary Information. Groups of descriptors having r 2 greater than 0.6 and p-values < 0.05 of the independent variables and intercept values were retained for further analysis. The same method was applied increasing the number of descriptors in combinations until all p-values of a particular number of descriptors in combination became greater than 0.05. This method was applied on both series of compounds separately and each of the selected group of descriptors formed candidate model for the respective congeneric series.

Multivariable r 2 values of a number of candidate models were significantly high requiring additional criteria for model selection. In order to rank these candidate models AIC (Akaike Information Criterion) analysis28 was applied. AIC with a second order correction for small sample size is given by

AICC=nlnRSSn+2k+2k(k+1)nk1
(3)

where n = sample size, k = number of variable including the intercept, RSS = residual sum of squares.

For larger sample sizes the last term of the eq 2 vanishes to yield the original AIC equation. The model having the smallest AIC value among all the candidate models is designated as the best model for the given data. The AIC weight factor (Wi) is used to measure the probability of the model having the smallest AIC value to be the best in the entire candidate model set. Wi is calculated using the following equation:

Wi=exp(Δi/2)i=1Nexp(Δi/2)
(4)

where N = number of candidate models and Δi = AICi – min(AIC).

The robustness of the predictive power of these two models was tested by leave-one-out cross validation analysis. In this method one of the observations was kept as validation data with the rest of the data used to form the training set. The permeability of the test compound was then predicted using the model based on the training set compounds. This procedure was repeated for all other compounds until each of them served once as a test compound. The predictive power of the model was then assessed by calculating the cross-validated r 2 or Q 2 using the following equation:

Q2=1(ypredyobs)2(yobsymean)2
(5)

Results and Discussion

The goal of the present study was the development of models to quantitatively relate experimental passive permeability to calculated molecular descriptors for drug-like molecules. Such models would allow for a better understanding of rate limiting steps to passive diffusion for the series of compounds studied as well as have the potential to be used a predictive models. For the study three series of compounds based on the benzoic acids, pyridines and quinolines were selected. There were fifteen different chemical groups used to form the congeneric series of benzoic acid, e.g., −Cl, −Br, −F, −OH, −NH2, −COCH3, −CH3, −CH2CH3, −OCH2CH3, −CH2CH2CH3, −OCH2CH2CH3, −CH2CH2CH2CH3, −OCH2CH2CH2CH3, −CH2CH(CH3)2 and −H. The series of pyridine was made of ten different chemical substituents including −COOH, −NH2, −F, −Cl, −Br, −CH3, −CH2CH3, −COCH3, −OCH3, −CHO and −H. Only three substituents (-CHO, -Br and -H) were used for the quinolines. The mono-substituted congeneric series of the parent compounds were obtained by placing the substituting groups at 2-, 3- or 4-position on the aromatic ring.

Passive permeability of the congeneric compounds was measured across parallel artificial membrane to study the effect of the presence as well as the location of chemical substituent on their cellular permeability. The experimental values of P pass of the compounds are listed in Table 1, Table 2 and Table 3 for the benzoic acids, pyridines and quinolines, respectively. Unsubstituted benzoic acid showed poor bilayer permeability (24.5 × 10−8 cm/s) as opposed to pyridine (1880 × 10−8 cm/s) and quinoline (2080 × 10−8 cm/s). Most of the chemical substituents had enhancing effect on P pass of the benzoic acid while the similar substituents had the opposite effect in most of the cases for pyridine and quinoline.

Table 1
Benzoic Acid Derivatives and Their Mean Permeability Values
Table 2
Pyridine Derivatives and Their Mean Permeability Values
Table 3
Quinoline Derivatives and Their Mean Permeability Values

Experimental Results for the Benzoic Acid Series

All bulky alkyl and alkoxy groups (e.g. −OCH2CH3, −CH2CH2CH3, −OCH2CH2CH3, −CH2CH2CH2CH3, −OCH2CH2CH2CH3, −CH2CH(CH3)2 ) enhanced the permeability of the benzoic acid significantly. In case of smaller alkyl groups (e.g. −CH3) the substitution location played the major role; 2- and 3- position had very small reducing and enhancing effects, respectively, on permeability while the 4-position showed a significant increasing effect. All the halogen substituents decreased the permeability of benzoic acid except for 4-F, 4-Cl, 3-Br and 4-Br. Polar groups like −NO2, −NH2 and −OH reduced the permeability. Acetyl substitution reduced the permeability at the 2- and 3- positions while at the 4-position it had negligible effect.

Experimental Results for the Pyridine–Quinoline Series

Except for the 2- and 3-CHO substitutions all other substituents had a reducing effect on permeability for the quinoline compounds. 2-Br, 4-CH3 and 4-CH2CH3 were the only substituents leading to an increase in permeability for the pyridine compounds. Alternatively, the amino and acid substituents, which are expected to be charged at the experimental pH of 6.5, led to the largest decrease in permeability. Notably, these compounds had permeabilities similar to those observed in the benzoic acid series, indicating the impact of charge on membrane crossing.

Computational Model

Logarithms of the permeability values were used for model development since according to the Arrhenius transition state model barrier crossing has a logarithmic relationship with respect to the rate constant of the chemical process.29 No compound having zero or nonmeasurable P pass was included as mathematically the logarithm of zero approximates to negative infinity. Various molecular descriptors were correlated individually and in selected groups with the membrane permeability (log P pass). These physicochemical descriptors were chosen to elucidate explicit solvent−solute interactions as it may be anticipated that such terms contribute to the rate limiting step for cellular permeability of drugs. The calculations were performed on two separate groups of compounds: one consisting of the congeners of pyridines and quinolines, and the other one had the benzoic acid congeners.

The ionization state of the compounds was determined based on the experimental pH (6.5). The carboxyl group of benzoic acid (pK a ~ 4.19)30 would typically be ionized at pH 6.5. Accordingly, initial model development was based on the ionized species. In addition, model development for the benzoic acids was also performed with using the neutral state. At the experimental pH the amine group on pyridine is partially protonated as its pK a ranges from 7 to 9 depending on the position on pyridine. Experimental data show that amine substitution reduces P pass of pyridine indicating the probable polar character of the amine substituted compounds (Table 2). In order to identify the correct ionization state both protonated and neutral amines were studied, and protonated amino substituents showed better correlation with the log P pass for pyridine congeners. Hence, protonated amine substituents were selected for the final pyridine series. 2-NH2 was the only amine substitution in the benzoic acid series having nonzero permeability, and it is only 0.4% ionized (pK a ~ 4.9) at the experimental pH. Congeneric series of pyridines contained three −COOH substituents at the 2-, 3- and 4-positions and the pK a value of them ranges from 0.14 to 2.2. At the experimental pH these −COOH groups will be ionized. The reducing effect on permeability of 3- and 4-COOH reflects the anionic state of the substituents. On the contrary 2-COOH had very minimal effect on the permeability of the parent pyridine compound, indicating that at 2-position carboxylate group may retain its neutral form. Based on these assumptions deprotonated carboxylate groups were used for 3- and 4-substitutions while neutral −COOH was used for the substitution at 2-position. Quinoline series did not include any carboxylate substitution.

Four different configurations of protonated benzoic acid were studied to account for the asymmetry of the substituted benzene ring due to the various substituents on the ring as well as for the E and Z configurations of the proton on the acid group. Figure 2 shows the four initial protonated structures of the benzoic acid, and R1, R2 and R3 represent the substituents at 2-, 3- and 4-positions, respectively. Substituents were placed at the assigned positions for each of the four protonated configurations of benzoic acid to form four sets of benzoic acid congeneric series. Each structure was then energy minimized, and the free energy of solvation was calculated for each of them. Out of the four different configurations for each benzoic acid congener the configuration with the most favorable free energy of aqueous solvation was selected for further study.

Initial model development for the benzoic acid series initially assumed the acid group to be in its ionized state. Linear correlation analysis showed the highest correlation to occur with SASA followed by the cavity dispersion term in hexane and polar solvent accessible area terms being the next best predictors (Table 4). Multivariable regression was then performed yielding seven different combinations of descriptors with r 2 values greater than 0.6 and p-values less than 0.05 (data not shown). Of these the AICc showed only one to be important with 99% probability of being the best model, respectively. The top two models included the logD(o/w), H_don, SASA, CDS-hex and CDS-wat (Table 5). In the first model, the coefficients of logD(o/w) and H_don are positive and negative respectively indicating increase in polar character disfavors permeability. For the second model the coefficient of the SASA term are positive, indicating increased solvent accessible surface to increase permeability; the coefficient with CDS-hex was negative indicating that the decreases in this term lead to increased permeability; and the positive coefficient of CDS-wat indicated increase this term enhances permeability. In both models there is a tendency for terms that favor nonpolar solvation to increase permeability while those that favor aqueous solubility tend to lower permeability. These trends suggested that partitioning of the benzoic acids into the nonpolar region of the lipid bilayer is impacting the permeability. As the charge on the acidic group would be highly unfavorable, this results suggested that the neutral state of the benzoic acids may be more relevant for passive diffusion across membranes. Accordingly, model development was also performed with the neutral forms of the benzoic acids.

Table 4
Results from Single Variable Linear Regression of Molecular Descriptors for the Congeneric Series of Benzoic Acids with the Acid Group Assumed To Be Ionized.
Table 5
Calculation of Modified Akaike Information Criterion (AICC) and Weight Factor (Wi) for Various Combinations of Molecular Descriptors for the Congeneric Series of Benzoic Acids with the Acid Group Assumed To Be Ionizeda.

Ranking of the individual descriptors for the neutral benzoic acid series are shown in Table 6 along with the corresponding coefficients. The highest correlation occurs with logD(o/w), followed by the logP(o/w), ΔG-hex, SASA and the CDS-hex. The molecular descriptors were then combined in groups of two and more for multivariable regression analysis. The neutral benzoic acid series yielded eight different combinations of descriptors having r 2 greater than 0.6 and p-values less than 0.05 (data not shown). Among them two models were identified by AICC analysis to be important with 65% and 35% probability of being the best model, respectively. Both of the top two models included logD(o/w) with positive coefficients. First model also consisted of ΔG-h/w with negative coefficient. The positive coefficients of logD(o/w) indicated the greater the polar character the lower the permeability; consistent with this is the negative coefficients of ΔG-h/w. As ΔG-h/w is defined as the difference between ΔG-hex and ΔG-wat [(ΔG-hex) – (ΔG-wat)], the tendency of a molecule to partition into aqueous solution versus hexane (i.e. ΔG-h/w is positive) leads to a decrease in permeability. The second model also contains H_don. Negative coefficient of H_don was consistent with other terms as it suggested that increased polar character would inhibit permeability.

Table 6
Results from Single Variable Linear Regression of Molecular Descriptors for the Congeneric Series of Benzoic Acids with the Acid Group Assumed To Be Neutral

For the pyridine−quinoline series the correlation coefficients for the individual descriptor are shown in Table 8. Three terms show high correlation ΔG-h/w, ΔG-wat, and ΔG-hex. Multiple regression analysis for the combined series of pyridines and quinolines yielded fifteen candidate models (data not shown). Final models were again selected based on the AICC indicators, with 34%, 34%, 13% and 12% probability of being the best model, respectively (Table 9). ΔG-h/w, ΔG-wat, and ΔG-hex occurred at last once in each of the models, with other terms being the CDS-wat, SASA and MW. Coefficients with ΔG-wat and CDS-wat are all positive, indicating that as the free energy of solvation in aqueous solution becomes more favorable, the permeability is decreased. The negative coefficient with ΔG-h/w indicates that favorable aqueous to lipid transfer leads to improved permeability. The only term not directly indicating increased aqueous solvation to lower permeability is the positive coefficient with ΔG-hex. It may be interpreted as too much lipid solubility can trap the molecule inside the bilayer preventing it from entering the intracellular medium. Thus a balance between polar and non-polar character should be maintained so that the molecules can enter as well as leave the lipid layer. In the second model MW appeared with a positive coefficient suggesting that bulkier molecules would be favorable for permeability. To further verify the quality of the selected models for all three systems, Q 2 values for all the candidate models were calculated. The values are included in Table 5, Table 7 and Table 9 for the ionized benzoic acid, the neutral benzoic acids and the pyridine−quinoline series. In all the selected models for benzoic acid series with neutral acid group and pyridine−quinoline series, significantly high Q 2 (>0.7) are obtained indicating the high predictive power and the robustness of the models. Figure 3 shows the relation between experimental and predicted values of log P pass as obtained from the best models for the benzoic acids with a neutral acid group and the pyridine−quinoline series, respectively.

Figure 3
Experimental and predicted values of logPpass for the congeneric series of (A) benzoic acid (with neutral acid group) and (B) pyridine-quinoline as obtained from the respective best model
Table 7
Calculation of Modified Akaike Information Criterion (AICC) and Weight Factor (Wi) for Top Two Combinations of Molecular Descriptors for the Congeneric Series of Benzoic Acids with the Acid Group Assumed To Be Neutrala
Table 8
Results from Single Variable Linear Regression of Molecular Descriptors for the Congeneric Series of Pyridines and Quinolines
Table 9
Calculation of Modified Akaike Information Criterion (AICC) and Weight Factor (Wi) for Various Combinations of Molecular Descriptors for the Congeneric Series of Pyridines and quinolinesa

Mechanism of Permeability in PAMPA Using the Developed Models

The process of solute transport across lipid bilayer via passive diffusion is characterized by random movement of solute molecules, typically along a concentration gradient. In this process the solute molecules cross the bilayer from the extracellular region to the intracellular region. This process can be divided into the following steps:5, 31 (1) diffusion of the solute molecules across the aqueous boundary layer (ABL); (2) desolvation from the aqueous medium; (3) insertion of the desolvated solute molecules into the lipid layer, i.e., solute partitioning from aqueous to the hydrophobic lipid environment; (4) diffusion across the lipid bilayer towards the intracellular side of the bilayer; (5) solute partitioning from the membrane bilayer into intracellular aqueous medium.

The rate of passive diffusion or P pass may be limited by any of the above steps depending on the nature of the solute molecule crossing the bilayer. For example, step 2 could be the limiting step for a hydrophilic molecule. A lipophilic molecule should have greater ease to enter the hydrophobic region of the bilayer compared to a hydrophilic molecule; however, it cannot be too lipophilic as it will have a high probability of being trapped inside the lipid bilayer due to high affinity for the lipid interior and an unfavorable lipid−water partitioning energy (e.g., ΔG-h/w). Various physicochemical molecular descriptors were considered in this study to explain the process of passive permeation of monosubstituted benzoic acids, pyridines and quinolines. ΔG-wat is the solvation free energy of a solute molecule in water and an unfavorable value of ΔG-wat may indicate easier desolvation from the aqueous medium (step 2). Alternatively, a favorable value of ΔG-hex and a negative value of ΔG-h/w may indicate favorable partitioning from aqueous to the hydrophobic lipid phase (step 3).

Based on the physical properties that contribute to the final models for the three series of compounds suggestions as to the rate limiting steps can be made. To facilitate this analysis the linear r 2 values and coefficients from the linear regression in Table 4, Table 6 and Table 9 are combined in Table 10. The ordering of the physical properties in the Table is based on the r 2 values for the ionized benzoic acid series. LogD(o/w) and logP(o/w) were the top two descriptors occurred in both acidic and neutral benzoic acids. The positive coefficients of these terms suggested that greater hydrophobicity favored permeability. The coefficient with SASA is positive for all three congeneric series, with high correlation occurring with both the neutral and acidic benzoic acid series. The higher correlation of SASA in the benzoic acids is due to the larger hydrophobic moieties that are not present in the pyridine−quinoline series, leading to some significant changes in SASA. However, correlation between SASA and CDS-hex for the acidic and neutral series yielded correlation coefficients of –0.7 and –0.74, respectively (Table S4 and Table S5). Thus, it appears that the SASA contribution is actually related to hydrophobicity of the larger functional groups.

Table 10
Linear correlation coefficients and regression coefficients of the molecular descriptors for two types of benzoic acid series (acid group assumed to be ionized and neutral) and pyridine-quinoline series. Molecular descriptors are ranked in the descending ...

Of the remaining physical properties there is a trend where terms that may be related to aqueous solvation have coefficients that indicate increased aqueous solubility to decrease permeability. Examples include the positive coefficients with CDS-wat and ΔG-wat. These results indicate that the role of desolvation is contributing to the rate limiting step to permeability, corresponding to steps 2 and 3 in the scenario presented above. Consistent with this is the ionized benzoic acid series, where the dominant terms are the PSA, CDS-hex and ΔG-hex (Table 5). The coefficients of PSA, CDS-hex and ΔG-hex are all negative, such that increase polar surface area disfavor permeability while increased solvation in hexane favors permeability. Interestingly, the ΔG-h/w term does not make a significant contribution to the ionized benzoic acid series, though it does make a contribution to both the neutral benzoic acid and pyridine−quinoline series based on both the linear regression and multiple regression analyses. While speculative, these similarities along with the higher correlation coefficient for the neutral versus the ionized benzoic acid series suggest that similar physical events are dictating the rate limiting step for the neutral benzoic acids and the pyridine−quinoline series. Based on the scenario presented above this would involve step 2 (desolvation from aqueous solution) and 3 (insertion into the hydrophobic core of the lipid bilayer). The present data certainly does not allow for these two contributions to be separated and, indeed, the aqueous desolvation and membrane insertion may occur simultaneously.

With the benzoic acids, both the ionized and neutral states were used for model development. The higher r 2 values for both linear and multiple regression for the neutral series indicate this to be the more appropriate ionization state for modeling permeability and that the neutral state is that which is involved in the rate limiting step to membrane permeation. Notably, when the neutral species is used, many of the physical descriptors that have predictive abilities are common to the pyridine−quinoline series, suggesting that the barrier to permeation is similar. Why then are the permeabilities of the benzoic acids approximately 2 orders of magnitude lower than that of pyridine and quinoline (Table 1, Table 2 and Table 3)? From the Arrhenius equation it can be shown that the change in free energy associated with a two-order change in the magnitude of the permeability is approximately 2.7 kcal/mol. What is causing the barrier to be higher by 2.7 kcal/mol in the benzoic acid series? The pK a of unsubstituted benzoic acid is 4.2, which is significantly lower than the pH at which the experiments were performed, 6.5. Thus, for a significant population of the neutral species of benzoic acid to be achieved there would have to be a pK a shift of a least 2 units. The amount of energy required for such a shift may be determined from the following equation, as discussed in detail by Chen and MacKerell33.

ΔpKa=12.3RTΔΔG
(6)

where ΔΔG = change in the change in Gibbs free energy. As may be seen the amount of energy required to shift the pKa by two units is also 2.7 kcal/mol, due to the natural log relationship of the changes in equilibrium to energy based on the van’t Hoff equation. Similar observations have been reported by other research groups explaining the process of permeation of ionizable species through lipid bilayer via passive diffusion.3436

A schematic energy diagram of bilayer permeation by the benzoic acids is presented in Figure 4. In aqueous solution, the acid group is ionized, consistent with the pK a of 4.19.30 Upon reaching the aqueous−lipid interface and moving into the membrane (steps 2 and 3 above) the ionization state changes to the neutral form, requiring additional energy, thereby leading to a larger barrier and concomitantly lower permeabilities. This step is indicated as step 3’ in Figure 4. Consistent with this hypothesis is the experimental data showing that polar substituents that lower the pKa of substituted benzoic acid, causing a decrease in the percent of the neutral species as pH 6.5, lead to decreased permeability; on the contrary, nonpolar substituents that increase the pKa of the acid moiety enhance the permeability.37 For example, polar substituents like 2-OH and 2-Br decreased the pK a of benzoic acid (pKa=4.19)30 to 3.01 and 2.9,27, 30 respectively and in turn, reduced the permeability of benzoic acid from 24.5 × 10−8 cm/s to 1.95 × 10−8 cm/s and 4.96 × 10−8 cm/s, respectively (see Table 1). On the other hand, nonpolar substituents like 4-CH2CH3 and 4-CH2CH2CH3 increased the permeability to 218 × 10−8 cm/s and 652 × 10−8 respectively (see Table 1) consistent with increased the pK a values of 4.35 and 4.37,27 respectively. Further supporting this hypothesis are the permeabilities of the charged pyridine analogs, which have permeabilities approaching 2 orders of magnitude below that of many of the neutral compounds. Thus, it is suggested that the generally decreased permeability of the benzoic acids as compared to the pyridines and quinolines is associated with the need for the acidic group to be protonated, a process that requires energy therefore effectively raising the barrier to permeation. This process is suggested to also contribute to the decreased permeability of the charged pyridine analogs.

Figure 4
Schematic diagram of the steps associated with permeation across a lipid barrier of benzoic acid. The steps as shown are (1) diffusion across the aqueous boundary layer, (2) aqueous desolvation, (3) insertion into the hydrophobic interior of the lipid ...

Conclusion

Developed models relating computed physical properties to permeability for the two sets of congeneric series indicated the importance of a proper balance in the polar and nonpolar character of a molecule to be able to possess a high passive permeability. The various steps associated with the process of passive diffusion were greatly reflected in the selected models. Since a molecule needs to pass through an aqueous layer as well as the hydrophobic interior of a lipid layer at different stages of the process, it is evident that either absolute polar or absolute nonpolar nature would hinder its permeability. Instead an appropriate equilibrium of these two opposing characters would help the molecule to efficiently diffuse from the extracellular region to intracellular region through the lipid bilayer. The process of passive permeability also depends on the effective pK a of the ionizable chemical moieties present in the molecules. The congeneric series of benzoic acid presented an interesting example of the importance of pK a for this process, where the present data suggested a model where in order to cross the bilayer membrane the pK a of the acid moiety of the benzoic acids must shift to obtain a significant population of the neutral species at pH 6.5. This process requires additional energy input, leading to the decreased permeation rates for the benzoic acids as well as for pyridines with charged chemical substituents.

Supplementary Material

supp data

Acknowledgements

NIH grant (DK067530) is acknowledged for financial support and the University of Maryland Computer-Aided Drug Design Center, the Advanced Biomedical Computing Center and the Pittsburgh Supercomputing Center are acknowledged for computational support.

Reference

1. Egan WJ, Lauri G. Prediction of intestinal permeability. Adv Drug Deliv Rev. 2002;54(3):273–289. [PubMed]
2. Wessel MD, Jurs PC, Tolan JW, Muskal SM. Prediction of human intestinal absorption of drug compounds from molecular structure. J Chem Inf Comput Sci. 1998;38(4):726–735. [PubMed]
3. Hou T, Wang J, Zhang W, Wang W, Xu X. Recent advances in computational prediction of drug absorption and permeability in drug discovery. Current Medicinal Chemistry. 2006;13:2653–2667. [PubMed]
4. Walters WP, Murcko MA. Prediction of 'drug-likeness'. Adv Drug Deliv Rev. 2002;54(3):255–271. [PubMed]
5. Chen IJ, Taneja R, Yin D, Seo PR, Young D, MacKerell AD, Jr, Polli JE. Chemical substituent effect on pyridine permeability and mechanistic insight from computational molecular descriptors. Mol Pharm. 2006;3(6):745–755. [PMC free article] [PubMed]
6. Barbosa F, Horvath D. Molecular similarity and property similarity. Curr Top Med Chem. 2004;4(6):589–600. [PubMed]
7. Leach AR, Hann MM. The in silico world of virtual libraries. Drug Discov Today. 2000;5(8):326–336. [PubMed]
8. Stenberg P, Norinder U, Luthman K, Artursson P. Experimental and computational screening models for the prediction of intestinal drug absorption. J Med Chem. 2001;44(12):1927–1937. [PubMed]
9. Subramanian G, Kitchen DB. Computational approaches for modeling human intestinal absorption and permeability. J Mol Model. 2006;12(5):577–589. [PMC free article] [PubMed]
10. Clark DE, Pickett SD. Computational methods for the prediction of ‘drug-likeness’ DDT. 2000;5(2):49–57. [PubMed]
11. Stenberg P, Luthman K, Artursson P. Virtual screening of intestinal drug permeability. J Control Release. 2000;65(1–2):231–243. [PubMed]
12. Pillai AD, Rathod PD, P XF, Patel M, Nivsarkar M, Vasu KK, Padh H, Sudarsanam V. Novel drug designing approach for dual inhibitors as anti-inflammatory agents: implication of pyridine template. Biochem Biophys Res Commun. 2003;301(1):183–186. [PubMed]
13. Nazare M, Matter H, Klingler O, Al-Obeidi F, Schreuder H, Zoller G, Czech J, Lorenz M, Dudda A, Peyman A, Nestler HP, Urmann M, Bauer A, Laux V, Wehner V, Will DW. Novel factor Xa inhibitors based on a benzoic acid scaffold and incorporating a neutral P1 ligand. Bioorg Med Chem Lett. 2004;14(11):2801–2805. [PubMed]
14. Gill A. Protein kinases in drug discovery and development. Drug Discov Today. 2004;9(1):16–17. [PubMed]
15. Kansy M, Senner F, Gubernator K. Physicochemical high throughput screening: parallel artificial membrane permeation assay in the description of passive absorption processes. J Med Chem. 1998;41(7):1007–1010. [PubMed]
16. Thompson M, Krull UJ, Worsfold PJ. The structure and electrochemical properties of a polymer-supported lipid biosensor. Analytica Chimica Acta. 1980;117:133–145.
17. Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983;4:187–217.
18. Shao Y, Fusti-Molnar L, Jung Y, Kussmann J, Ochsenfeld C, Brown ST, Gilbert ATB, Slipchenko LV, Levchenko SV, O'Neill DP, DiStatio RA, Jr, Lochan RC, Wang T, Beran GJO, Besley NA, Herbert JM, Lin CY, Van Voorhis T, Chien SH, Sodt A, Steele RP, Rassolov VA, Maslen PE, Korambath PP, Adamson RD, Austin B, Baker J, Byrd EFC, Dachsel H, Doerksen RJ, Dreuw A, Dunietz BD, Dutoi AD, Furlani TR, Gwaltney SR, Heyden A, Hirata S, Hsu C-P, Kedziora G, Khalliulin RZ, Klunzinger P, Lee AM, Lee MS, Liang W, Lotan I, Nair N, Peters B, Proynov EI, Pieniazek PA, Rhee YM, Ritchie J, Rosta E, Sherrill CD, Simmonett AC, Subotnik JE, Woodcock HL, III, Zhang W, Bell AT, Chakraborty AK, Chipman DM, Keil FJ, Warshel A, Hehre WJ, Schaefer HF, III, Kong J, Krylov AI, Gill PMW, Head-Gordon M. Q-CHEM, Version 3.1. Pittsburgh, PA: Q-CHEM, Inc; 2007.
19. Hawkins GD, Geisen DJ, Lynch GC, Chambers CC, Rossi I, Storer JW, Li J, Thompson JD, Winget P, Rinaldi D, Liotard DA, Cramer CJ, Truhlar DG. AMSOL Version 6.8. Minneapolis, MN: University of Minnesota;
20. Lee B, Richards FM. The interpretation of protein structures: estimation of static accessibility. J Mol Biol. 1971;55(3):379–400. [PubMed]
21. Chothia C. Structural invariants in protein folding. Nature. 1975;254:304–308. [PubMed]
22. Auffinger P, Hays FA, Westhof E, Ho PS. Halogen bonds in biological molecules. Proc Natl Acad Sci U S A. 2004;101(48):16789–16794. [PubMed]
23. Berski S, Ciunik Z, Drabent K, Latajka Z, Panek J. Dominant role of C-Br---N halogen bond in molecular self-organization. Crystallographic and quantum-chemical study of schiff-base-containing triazoles. Journal of Physical Chemistry B. 2004;108:12327–12332.
24. Panigrahi SK, Desiraju GR. Strong and weak hydrogen bonds in the protein-ligand interface. Proteins: Structure, Function, and Bioinformatics. 2007;67(1):128–141. [PubMed]
25. MOE (Molecular Operating Environment), Version 6.0. 1010 Sherbrooke Street West, Montreal, Canada H3A2R7: Chemical Computing Group; http://www.chemcomp.com)
26. Scherrer RA, Howard M. Use of distribution coefficients in quantitative structure-activity relationships. Journal of Medicinal Chemistry. 1977;20(1):53–58. [PubMed]
27. SciFinder Scholar 2007 Database. Washington DC, USA: American Chemical Society;
28. Akaike H. Likeliood of a Model and Information Criteria. Journal of Econometrics. 1981;16:3–14.
29. IUPAC Compendium of Chemical Terminology. 2nd Edition
30. CRC Handbook of Chemistry and Physics. 84th Edition. CRC PRESS; 2003–2004.
31. Walter A, Gutknecht J. Permeability of small nonelectrolytes through lipid bilayer membranes. J Membr Biol. 1986;90(3):207–217. [PubMed]
32. Feller SE, Venable RM, Pastor RW. Computer Simulation of a DPPC Phospholipid Bilayer: Structural Changes as a Function of Molecular Surface Area. Langmuir. 1997;13(24):6555–6561.
33. Chen IJ, MacKerell AD. Computation of the influence of chemical substitution on the pKa of pyridine using semiempirical and ab initio methods. Theor Chem Acc. 2000;103:483–494.
34. Hamilton JA. Fatty acid transport: difficult or easy? J Lipid Res. 1998;39(3):467–481. [PubMed]
35. Kessel A, Musafia B, Ben-Tal N. Continuum solvent model studies of the interactions of an anticonvulsant drug with a lipid bilayer. Biophys J. 2001;80(6):2536–2545. [PubMed]
36. Ulander J, Haymet AD. Permeation across hydrated DPPC lipid bilayers: simulation of the titrable amphiphilic drug valproic acid. Biophys J. 2003;85(6):3475–3484. [PubMed]
37. Hollingsworth CA, Seybold PG, Hadad CM. Substituent effects on the electronic structure and pKa of benzoic acid. Internation Journal of Quantum Chemistry. 2002;90(4–5):1396–1403.