The work described herein demonstrates robust consistency between the in vitro PAMPA approach and the predictive capacity of MD calculations. Our predictive MD approach is calibrated using the equivalent PAMPA experimental measurements. The semiquantitative predictive power of the methodology is then rigorously tested using a custom synthesized set of structurally related compounds (compounds LLNL1–LLNL18), demonstrating the potential utility of these methods for permeability assessment of novel chemical entities.

One of the most powerful (though computer intensive)techniques to simulate the molecular process of diffusion at the atomic level is molecular dynamics (MD). (25) Coupling MD with free energy techniques provides a powerful tool to study membrane permeability in detail. (26) While computational models can often correlate well with experiment, knowledge-based QSPR methods are often less accurate (but inexpensive) compared to MD methods. (27) Although human membranes are a complicated mixture of lipids and proteins, with a judicious choice of lipid, a single bilayer can provide a good first approximation of the physicochemical properties of a cellular membrane. Indeed, phosphatidylcholine lipids are the major phospholipid within cellular membranes. (28) Previous studies have investigated the permeation of small compounds (and a limited number of drug-like molecules) through homogeneous lipid bilayers. (29-35) Thus, with increasing computing power, these methods will become an increasingly valuable tool for parameter prediction.

For prediction of passive membrane permeability, the most critical parameter in QSPR studies is lipophilicity. (20) Lipophilicity is a crucial factor governing passive membrane partitioning, (21) and calculated lipophilicity metrics are often utilized to predict drug absorption. (22) The parameter that determines the lipophilicity of a molecule is LogP (the partition coefficient of the molecule between an aqueous and lipophilic phase, usually water and octanol). LogP is a crucial factor governing passive membrane partitioning; an increase in LogP enhances permeability. (21) While the partition coefficient is used to calculate properties such as membrane permeability and water solubility, it also has importance in the prediction of biological activities, ADME, and toxicological end points. Thus, reliable prediction of LogP is of massive importance in the drug design process, as it is important to know lipophilic properties of a compound before synthesis. (23) There are several methods used to calculate the LogP parameter for a initial, rapid prediction of lipophilicity and, thus, permeability. Atomic-based LogP predictions take each atom’s contribution to the LogP of the compound and combine them additively. (24) This method is in essence a look-up table per atom and suited to smaller, less complex molecules. More advanced “hybrid” LogP methods use the LogP contribution of each atom, as well as a contribution from the neighboring atoms, and combine them with correction factors. In a fragment-based approach, data determined experimentally from different chemical moieties, or chemical “fragments” are modeled, and these contributions are again added up with correction factors. The motivation is that atomic-based approaches do not always properly capture certain molecular subtleties. Thus, fragment-based methods tend to be better for complex, larger molecules. However, the assumption is the molecule contains fragments that are similar to those from which the model was constructed.

An alternative approach to assess membrane permeability is to employ computational methods. Mostprediction approaches use quantitative structure-permeability relationship (QSPR) models. These predictive, quantitative models applied in drug design are methods that utilize statistical relationships extracted from experimental permeability measurements combined with the physiochemical properties of a set of training compounds. Thus, in QSPR studies, the permeability of the compound is an outcome that is the result of the various interactions that a compound experiences during passage through the membrane. The interactions between the compound and the membrane are assumed to be governed by the chemical structure and properties of the compound, which are termed descriptors. As such, a mathematical model of biological permeability is optimized based on a combination of a variety of descriptors for the small molecule compound. Thus, a change in structure can result in a change in permeability. QSPR models have been developed to model and predict, among others, oral bioavailability, (14) intestinal absorption, (15) Caco-2 permeability, (16) and blood-brain barrier penetration. (17) The degree of success of these models is extremely dependent on the compounds in the training set. Thus, the applicability of the approach can be limited, and transferability can be a major issue. (18) Even successful QSPR models can be limited because they do not provide any insight into the mechanisms that govern the permeability. (19)

Membrane permeability is a key metric in the drug design pipeline. A drug intended for an intracellular target, but with poor permeability will have low efficacy. As a result of these characteristics several well-characterizedandpermeability prediction methods have been developed. These methods are relatively high throughput, and are important in the early stages of the drug discovery process with some of the most common and relatively simplemethods being the parallel artificial membrane permeability assay (PAMPA), the immobilized artificial membrane (IAM) technique, and immobilized liposome chromatography (ILC). The PAMPAtechnique was developed in 1998 by Kansy et al. (4) and originally established to rapidly predict passive permeability through the gastrointestinal tract, but has since been adapted for use in other systems such as the blood brain barrier (BBB), (5) for which it demonstrated a good prediction of BBB penetration. (5, 6) Briefly, the technique involves a donor and an acceptor compartment separated by a filter supporting a liquid artificial membrane. The artificial membrane can be composed of a variety of phospholipid mixtures. The compound to be tested is placed in the donor compartment and is allowed to permeate between the donor and the acceptor compartments through the artificial membrane. As the assay can be performed in 96-well plates, high throughput screening of drug candidates is possible. Alternatively, IAMs mimic the phospholipid environment of the cellular plasma membrane by anchoring synthetic lipid analogues to silica particles in monolayer density. These particles are then used as the packing material for a high-performance liquid chromatography column. (7-10) The IAM chromatographic retention factors are used to generate predictions of membrane permeability. These systems have shown reasonable results for prediction of permeability, despite the retention time in the column not reflecting actual passage across the membrane. (11, 12) The ILC method is used to study solute-membrane interactions, where liposomes are noncovalently immobilized to gel beads. The liposomes can have varied compositions, including numerous different phospholipids. Lipids extracted from human red cells have also been used in ILC assays. (13)

Most drugs need to pass through at least one cellular membrane to reach their intended target. Although tight binding of a drug molecule to its intended target is important for potency, poor membrane permeability often translates into poor or nonexistentefficacy. Thus, a detailed understanding of the partitioning of a given species in the membrane is vitally important from a pharmacokinetics and rational drug design standpoint. In eukaryotic systems, two possible transport modes are available for a molecule to pass through a membrane: active and passive. (1) Active transport involves a transport protein that uses energy (e.g., ATP hydrolysis) to shuttle a molecule across a membrane. In contrast, passive transport, which is the most common mode of drug passage through membranes, involves diffusion of a molecule through the membrane with no outside assistance or energy input. The rate of passive diffusion across a membrane is proportional to the partition coefficient of the compound between the membrane (lipophilic environment) and the external medium (aqueous milieu), the diffusion coefficient of the compound through the membrane, and the compound’s concentration gradient across the membrane. (2) Important chemical properties of small molecules for the process of membrane binding and diffusion are lipophilicity, molecular weight, and measures of molecular polarity. (3) Thus, the development of a successful drug involves a fine balance among all these properties in a molecular scaffold, which is no trivial matter.

The LogP calculations (termed “(c)LogP”) implemented from the Chemaxon/Chemicalize server (48) are based on a pool of predefined fragments. This pool of fragments is based on the data in Viswanadhan et al., (59) with a few minor adaptations (such as redefining atomic types to accommodate electron delocalization and contributions of ionic forms, or considering the potential for hydrogen bonding to form a six membered ring between suitable donor/acceptor atoms).

The Molinspiration method for LogP prediction (miLogP) is also based on group/fragment contributions. The fragment contributions were obtained by fitting a training set of >12 000 mostly drug-like molecules with experimental LogP data. The authors state that miLogP is calculated using the methodology developed by Molinspiration as a sum of fragment-based contributions and correction factors. (57, 58)

CLogP is a proprietary method that is owned by BioByte/Pomona College. This method is a fragment-based approach, whereby experimentally determined fragments are modeled using QSPR (rather than per atom). The CLogP program (54-56) has become the benchmark, or gold-standard, method for LogP prediction in the last 35 years and is the most extensively tested fragment-based LogP calculator used in drug design.

The concentration data was used to calculate effective permeability (cm/s) as follows:whereis the filter insert area (0.3 cm),is the volume of the donor well (0.3 mL),is the volume of the acceptor/receiver well (0.2 mL),) is the concentration in the acceptor well at time, and= incubation time (18 000 s).was calculated according to the following equation:where) is the concentration in donor well at time t, and the other values are the same as previously defined.

Compounds in solution were quantified using the Acquity ultra performance liquid chromatography (UPLC) system (Waters) with the Acquity BEH C18 column (1.7 μm × 2.1 mm × 50 mm) at a pressure of 7500 PSI. Here, 10 μL were injected for each sample. Detection was based on retention time and UV absorbance. Specific detection protocols were determined, where necessary for each compound ( Supporting Information ).

PAMPA was applied to screen compounds for passive diffusion. This study employed the Gentest Precoated PAMPA Plate System (Corning Discovery Labware). The system is composed of a 96-well plate/insert system. Two fluid-filled compartments (donor well and receiver well) are separated by a polyvinylidene fluoride (PVDF) filter plate precoated with a phospholipid-oil-phospholipid trilayer primarily consisting of DOPC phospholipids. (53) Compounds of interest were dissolved in Hanks Balanced Salt Solution (HBSS) at a concentration of 100 μM. Solutions were added to donor wells at a volume of 0.3 mL. The filter plate, containing acceptor wells filled with 0.2 mL blank HBSS, was placed on top of the donor plate. The entire system was incubated for 5 h at 25 °C. Following incubation, solution was removed from donor and acceptor wells. Compounds were then reserved for analysis via ultra performance liquid chromatography (UPLC). Certain compounds exhibited very low traversal to the acceptor well, such that material could not be detected via UPLC upon initial analysis. Such compounds were concentrated from donor and acceptor solution fractions by extraction into methanol, followed by vacuum concentration and resuspension in the minimal amount of appropriate UPLC solvent.

In order to determine the error of the, the data in the 50 diffusion profiles were combined with the data in the 50 PMF profiles using eq 3 to produce 50 similar, but different,s. The standard deviation for this set of 50s was taken as the error.

Resampling of the PMF curves was calculated using the bootstrapping method as implemented within the g_wham (52) analysis program of GROMACS. Again, 50 separate PMF profiles were generated, from which the mean and standard deviations calculated. These standard deviations of the PMF were used as an additional quality control metric. Any PMF profile that exhibited a standard deviation of >0.5 kcal molwas subjected to extended simulation and supplementary sampling in order to reduce this deviation. The PMF histogram overlap was also checked to ensure correct sampling.

As a metric for confidence in our results, the error of P eff was determined. In order to calculate this overall error of P eff , we first had to separately resample both the PMF and diffusion profiles. The sensitivity analysis used a random resampling of the position-specific diffusion coefficients. As described previously, for each compound, position-specific diffusion coefficients were calculated as well as the standard deviation of each position-specific coefficient. These calculations were based on independent subsamples of the trajectory from each umbrella-sampling simulation. Since the mean and variation of each calculated coefficient is an independent calculation, the correlation of the variations is zero. The resampled diffusion coefficients were thus based on independent draws from a set of normal distributions, where each such distribution has a mean and standard deviation corresponding to the calculated mean and standard deviation of each position-specific diffusion coefficient. We used 50 random draws for each coefficient to produce a set of 50 diffusion profiles that were distributed around the calculated data according to the normal distribution of that data.

The membrane permeability of a small molecule was calculated using the same methodology as previously published. (36) Briefly, information from the PMF 1D energy landscape was combined with diffusion coefficients within the membrane to calculate the membrane permeability rates. The position-specific diffusion coefficients were calculated from the MD simulation data using the method of Hummer. (50) For each umbrella sampling simulation, the compound positional variance and autocovariance as a function of lag time were calculated. These calculations were repeated for a number of subsamples of each trajectory. The resulting autocovariance curves decay roughly exponentially with increasing lag time. The characteristic time, τ, of each autocovariance decay was estimated by making a least-squares fit to the log of the autocovariance data. The diffusion coefficient for each subsample was then calculated as:where ⟨⟩ is the average of theaxis position of the compound center of mass during the simulation and) is the variance of the compound center of mass (the autocovariance at “lag zero”), and the average over the subsamples was calculated. The standard deviation of the diffusion coefficient over the subsamples was also calculated for use in a sensitivity analysis, described below. The effective membrane permeability,, was calculated from the effective resistivity,, (as derived by Marrink and Berendsen (51) ),where) is the resistivity of every “slice” of the membrane at position “”,where β is the inverse of the Boltzmann constant times the temperature, and Δ) is the free energy from the PMF calculations. The integral is over the width of the membrane.

The potential of mean force (PMF) free-energy profiles for the partitioning of compounds was calculated using umbrella sampling simulations. A single harmonic restraint with a force constant of 1000 kJ molnmwas applied to the-axis (normal to the bilayer) distance between the center of mass of the compound and the center of mass of the DOPC bilayer. One hundred separate simulations were performed, with the compound harmonically restrained in increments of 0.1 nm along the-axis direction. These 100 simulations completely span the 10 nm distance from bulk water, through the entire membrane and out into bulk water again. Each simulation window was initially run for ∼50 ns, for a total of ∼5 μs of simulation time for a complete set. As the position of a molecule within the membrane can shift its pvalue, most compounds were simulated in both their neutral forms and most physiologically relevant charge species, resulting in ∼10 μs of simulation data for each compound. A total of 27 different compounds were simulated for a combined ∼250 μs of simulation time. On the basis of our previous experience, the initial 20 ns of each simulation was discarded as equilibration time and as such only the final ∼30 ns of each simulation was used for all subsequent analyses. The weighted histogram analysis method, (49) as implemented within GROMACS, was used to calculate the PMF for each compound, and all free-energy profiles were normalized to zero in bulk water.

The simulation setup and protocol was the same as our previous published investigation. (36) In brief, each simulation contained a single copy of the compound harmonically restrained at specified locations along the-axis in a solvated DOPC bilayer system. The complete system contained 5124 water molecules, 72 DOPC molecules, and one small compound. A typical system contained a total of ∼19 300 atoms. This system is comparable to similar studies. We acknowledge that larger membrane systems may be required to investigate the permeability of larger molecules. The simulations were run using GROMACS 4.5.5, (37) with the Berger et al. force field used for the DOPC molecules (38) as adapted by the Tieleman group ( http://wcm.ucalgary.ca/tieleman/downloads ), the GROMOS 53A6 force field used for the small compounds, (39) and the SPC model used to represent the water, (40) with the simulation parameters the same as Carpenter et al. (36) The topologies for the small molecule compounds were generated using the Automated Topology Builder server ( https://atb.uq.edu.au ). (41-43) The topology builder combines a knowledge-based approach with quantum mechanical calculations to produce parameters that are compatible and consistent with a given force field. During the initial stage, the molecule was optimized at the HF/STO-3G level of theory. The molecule was then reoptimized at the B3LYP/6-31G* level of theory (44-46) in implicit water. The charges were initially estimated by fitting the electrostatic potential using a Kollman–Singh scheme. (47) The Hessian matrix was then calculated and used to estimate the force constants for both the bonds and angles. All compounds were modeled in their neutral forms and other physiologically relevant charged states. The ChemAxon/Chemicalize server (48) was used to calculate the pvalues for the ionizable atoms of each compound, and thus determine the most prevalent charge species at physiological pH of 7.4. The parameters for the calibration data set are freely available from the Automated Topology Builder server.

The prediction set of compounds used in this study was synthesized via an amide forming reaction between an amine and an ethyl ester species. The library of generated materials can be assembled from an ethyl ester and a series of amines that serve as the point structural diversity origin. Since the carboxylic acid component (ethyl ester) is not an activated species, the amine was used in slight excess to it (1.2 equiv) and the mixture was heated to 70 °C in an overnight basis in ethanol. In some instances, the products precipitate out of the solution (purity >99% after ethanol washes during filtration) while in some cases flash column chromatography (3:7 → 7:3 EtOAc:hexanes) is needed for their purification. Purity of the compounds was assessed by 1 H NMR (DMSO- d6 ) and in all cases found to be >98%.

In this study, two sets of compounds were characterized, the calibration data set and the prediction data set. The calibration data set was selected to cover a wide range of experimentally known permeability values and is comprised of the following compounds: progesterone, chlorpromazine, promazine, atropine, diazepam, theophylline, 2-PAM, Hi-6, and MMB4 ( Figure 1 ). The calibration compounds were commercially obtained (Sigma-Aldrich). The prediction data set (compounds LLNL1–LLNL18) is a structurally related set of compounds that were designed and synthesized in order to provide a set of compounds with a narrow range of permeabilities with which to rigorously test our methodology and predictive capabilities.

The LogP prediction tools performed more poorly than our PMF-based methodology in both the correlation against the control compounds and the semiquantitative prediction. While the correlation coefficient for our PMF-based methodology is 0.97, the range for the LogP tools is only 0.44–0.75. The LogP tools have difficulty with the unusual structures of our negative controls that are known to be impermeable. Furthermore, all of the LogP tools have a significant number of false negative permeability predictions, with no method able to correctly identify more than two (of the eight) compounds that exhibit any form of experimentally categorized permeability. In fact, two of the methods (SLogP and miLogP) fail to identify any of these compounds that have categorized permeability. Indeed, regardless of calibration, each of the LogP calculation tools assign a higher LogP (greater permeability) to some compounds measured to be impermeable than those measured to have a medium (or even high) permeability. In comparison, our PMF-based methodology has no such occurrences and correctly identifies all of the compounds with experimentally measured categorized permeability. Indeed, even the four impermeable compounds that our method incorrectly categorize as low permeability compounds still all have lower predicted permeability that the four compounds correctly identified to have low permeability. Thus, a “manual recalibration” of the impermeable/low permeability threshold cutoffs would allow our data to be categorized such that all 18 compounds are predicted correctly. This is impossible to achieve with any of the LogP data sets. Additionally, the compounds miscategorized by the PMF method are only a single cutoff (predicts to be “low permeability” but measures as “impermeable”) removed from their correct region, while many of the LogP predictions are misplaced by two, or even three regions (predicts to be “impermeable” but measures as “high permeability”). As can be expected, the fragment-based approaches perform marginally better than the atomic/hybrid method. CLogP is the top performing LogP tool (albeit only slightly, with only a ∼ 60% success rate and only predicting 2/8 permeable compounds). This is again anticipated as CLogP is the current gold-standard tool for LogP prediction.

The semiquantitative permeability prediction results generated from our PMF-based methodology are compared with the equivalent calibration and predictions using several commonly used LogP calculation tools ( Table 3 ). LogP prediction tools are used in lieu of QSPR models, which are primarily built upon LogP values. These LogP calculation tools are the atomic/hybrid methodology of MOE’s SLogP, and the fragment-based methods of miLogP (from Molinspiriation), (c)LogP (from Chemicalize), and CLogP (from Biobyte/Pomona College). The details of these techniques are described in the Methods section. The LogP calculation tools are subjected to the same analytic approach as for theresults. The various LogP calculations are first calibrated against themeasurements for the control compounds to define LogP cutoffs between the four different permeability categories. These defined thresholds for each of the LogP tools are used to categorize our 18 custom compounds ( SI , Figures S1 and S2).

Using these cutoffs, we attempt to semiquantitatively predict the permeability of 18 of our synthesized set of structurally related compounds. This prediction is significantly more difficult than correlating the control compounds, as the range of permeabilities of these compounds is much smaller. Despite this, however, we were able to correctly predict the permeability category of 78% (14 of 18) of the novel compounds ( Figure 4 ). The four compounds that are assigned into the wrong permeability category are computationally predicted to have a “low permeability”, while experimentally they are classified as impermeable ( Table 2 ). Thus, we have obtained a few “false positive” results. However, we have successfully identified (and categorized) all eight of the compounds with some level of permeability (low, medium, or high). From a drug design point of view, a false positive is a better outcome than obtaining false negatives that actually do possess permeability, but are predicted to be impermeable.

The linear correlation between our predictedpermeability rates and those experimentally measuredis extremely good (= 0.97). Using these measured experimental permeability rates, as well as permeability data from literature, we are able to divide our data into four groups; high permeability (green), medium permeability (yellow), low permeability (orange), and impermeable (red). The boundaries between these regions are defined by using themidpoints between different compounds of known permeability categories (identified from literature). The boundary between the “impermeable” and “low permeability” groups is chosen such that the size of the “low permeability” region is uniformly distributed about thevalue for the “low permeability” control compound (theophylline). These cutoff values of the experimentally measuredare determined such that impermeable compounds have a logof ←6.14, low permeability compounds have a logof > −6.14 and < −5.66, medium permeability compounds have a logof > −5.66 and < −5.33, and high permeability compounds have a logof > −5.33. By fitting against the linear regression line, we can use thesecutoff values to define the equivalent cutoff values for our predicted log Figure 3 and Table 1 ). Thus, the cutoff values for our predictions are impermeable compounds, log< −2.05; low permeabilty compounds, log> −2.05 and <0.15; medium permeabilty compounds, log> 0.15 and <1.62; high permeability compounds, log> 1.62. These values are in agreement with our previous studies, where the boundary between impermeable/poorly permeable and highly permeable is defined as a logof ∼0. (36) By defining what thecutoffs and thresholds are for these different permeability categories, we are able to compute thes for novel compounds and use those values to predict the permeability categories for those compounds.

Figure 3. Comparison between the experimentally measured P eff (using PAMPA, x -axis) and the computationally predicted P eff (using the PMF, y -axis). The plots are divided into regions of differing permeability. A red region is defined as impermeable, an orange region is defined as having a low permeability, a yellow region is defined as having a medium permeability, and a green region is defined as having a high permeability.

The PMF values for the calibration compounds (along with diffusion data also extracted from the MD simulations) are used to determine(effective permeability) of the compounds (termed “”). Thesevalues are then compared to the equivalentmeasured from running the same compounds throughPAMPA assays (termed “”), which also measure passive permeability ( Figure 3 ).

Figure 2. PMF free energy curves for the passage of the control compounds from water ( x -axis >3 nm) to the center of the bilayer ( x -axis <1 nm). As the membranes are symmetrical, the absolute distances of the compound from the bilayer center are used generate these curves. Generally, a negative free energy at the bilayer center correlates with a more permeable compound.

The potential of mean force (PMF) free energy curves for these compounds is generated from the umbrella sampling MD simulations ( Figure 2 ), and in keeping with previous results, (36) it shows that (as a general rule) if a compound has a higher free energy in the middle of the bilayer, it will be less permeable.

We are focused on improving the predictive capability of our membrane permeability model based on our previous work. (36) In order to do this, we continue to run simulations and carry out experimental assay validation on our data set of “calibration” compounds ( Figure 1 ). This set of compounds is comprised of molecules that fit a specific set of criteria: (1) a significant amount of previous experimental work is available to compare against our results, and (2) the compounds must also cover a wide range of permeability values to evaluate the performance of our model across this range.

Our computational methodology is first calibrated against the well-established PAMPA protocol using a set of diverse, highly studied compounds. The predictive power of the approach is then tested with a synthesized set of structurally similar compounds and compared against the equivalent predictive power of readily available LogP calculators.

Discussion ARTICLE SECTIONS Jump To

in vivo efficacy and bioavailability of candidate drug compounds. A range of physiological membranes are relevant to such concerns. Permeability across gastrointestinal membranes is critical to bioavailability Assessing permeability is of critical importance to understanding and predictingefficacy and bioavailability of candidate drug compounds. A range of physiological membranes are relevant to such concerns. Permeability across gastrointestinal membranes is critical to bioavailability (62) following oral administration. Transdermal permeation must be evaluated for topical administration and management of skin pathologies. (63) Traversal across the BBB represents a particularly confounding challenge, and is critical to delivery of compounds requiring central nervous system access for treatment of infection, neurologic malignancies, and neurodegenerative disease, (64, 65) among other maladies. Computational prediction of permeability represents a need that has been historically challenging and is often limited by the requirement for a robust training set, which ultimately restricts applicability of the resultant model. This study sought to address this need by refining and prospectively validating an MD-based model, demonstrating effective prediction of permeability for candidate compounds across a physiological lipid membrane.

The results of our study show that the developed computational model can predict the PAMPA-defined permeability category of a compound with greater accuracy than compared LogP-based methods. This is demonstrated both when examining well-characterized calibration compounds exhibiting a wide range of permeation capacity, as well when prospectively studying a set of novel compounds with a more narrow dynamic range. The methodology is extremely successful in characterizing the nine selected calibration compounds (R2 = 0.97), and is able to correctly categorize ∼80% of our internal “prediction set” of 18 compounds.

P eff when the compound has a lower permeability versus a higher permeability (if there is an energy barrier at the bilayer center rather than an energy well).P eff . As expected, the correlation with the calibration compounds is higher than with the prediction set of compounds because the calibration compounds are structurally diverse and their experimental permeabilities are spread over ∼5 log units, while the prediction set are structurally similar and only spread over <2 log units. Calculations based on this prediction set are thus a challenging task, especially as the properties of many of the compounds are so similar. Furthermore, several of the experimental permeabilities of the prediction set are measured at the lower end of the sensitivity threshold of the PAMPA methodology, resulting in higher relative uncertainty in this region due to the inherent challenge of experimentally measuring very low levels of diffusion precisely ( Figure 4 ). Amplified relative uncertainty at the low end of the permeability spectrum is correspondingly observed in computational predictions; we have previously shown that equivalent errors in the PMF profile have a larger effect on the calculatedwhen the compound has a lower permeability versus a higher permeability (if there is an energy barrier at the bilayer center rather than an energy well). (36) Thus, any error in calculating the PMF for these poorly permeable compounds will be reflected in larger error in the

These results are a significant improvement in accuracy over predictions made using current high throughput LogP calculation techniques. Indeed, some of the LogP-based methods fail completely when attempting to categorize the permeability of the compounds. Depending on which LogP tool is used, the same compound could be characterized into three different regions. Furthermore, our PMF-based predicted permeabilities are not dependent on training sets, but only the physical properties recapitulated in the MD simulations. Our first-principles method is especially useful for a novel set of compounds for which there is no QSPR or experimental data available, and we do not run the risk of overfitting our model to a narrow spectrum of physicochemical parameters.

Limitations in permeability and, ultimately, bioavailability, of candidate drug compounds can substantially slow and derail the development process. Ideally, computational model systems would serve as a first pass screen for permeation. Existing methods such as QSPR and LogP are not computationally intensive, but are restrictive in their input capacity, and were shown in our study to produce a large proportion of false negatives, thereby running the risk of ruling out potentially promising leads. Although MD-based simulations are more computationally intensive than LogP calculations (a few days compared to a few minutes), the result is a more accurate, actionable model that does not discard valuable compounds. Further, the computational time required for the method described in our study is still less than the time required to synthesize, purify, and experimentally characterize the same novel compound.

From a drug design perspective, this predictive capability would facilitate compound evaluation by ruling out impermeable candidates with a demonstrably low false-negative rate. The resultant data would provide a more in depth characterization of the permeability of hits identified though high throughput virtual screening or analogues of existing lead compounds. Application of the described predictive methods would facilitate faster iterative improvement, allowing for selective retention of compounds with access to the intended physiological target. Such a model could represent an important component of a more comprehensive design pipeline, capable of evaluating leads independent of structure or novelty of the given compound.