Journal of Chemical Theory and Computation Article Figure 1. Upper left panel (a): density profile of DOPC bilayer along the normal to the lipid bilayer plane showing densities of the complete system (black), water (blue), DOPC (red), phosphates (magenta), cholines (cyan), carbonyls (green), and terminal carbons (dark blue). Lower left panel (b): free energy profile of coumarin along the DOPC bilayer normal calculated from constraint simulation with initial structures obtained by free simulation (CF). The calculated bilayer center penetration barrier, ΔGpen, and water/lipids barrier, ΔGwat, are labeled. The free energy profile was calculated for one bilayer leaflet and was symmetrized to the other one, the densities of the system were symmetrized along the middle of the bilayer. The vertical bins labeled by numbers denote four bilayer regions: 1 − low density of head groups (2.2−2.9 nm), 2 − high density of head groups (1.45−2.2 nm), 3 − high density of acyl chains (0.5−1.45 nm), and 4 − low density of acyl chains (0−0.5 nm). Right panel (c): structure of DOPC bilayer, together with snapshots of coumarin initial structures. Carbons are colored in cyan, oxygens red, and hydrogens white. The olive and blue balls represent DOPC phosphate and nitrogen atoms. believed to form the main barrier for most druglike molecules, which are often water-soluble.4,16 Molecular dynamics (MD) simulations can be used to estimate the equilibrium position of a drug in a lipid bilayer, its partition coefficients and diffusion coefficients simultaneously at subpicosecond and atomic resolution.4,15,17,18 The partition- ing and mean position can be well described by a free energy (ΔG) profile along the normal to the lipid bilayer,5,13,18,19 also known as a potential of mean force (PMF). In principle, such free energy profiles can be calculated from partitioning values obtained by long unbiased simulation. However, this approach only provides reliable free energy profiles when all states along the ΔG profile are thoroughly sampled. This is challenging for unbiased MD simulations, because they usually do not sample adequately at the available simulation time scales. The sampling problem is based on the dependence of probability of drug crossing a membrane on the energy barrier for this phenomenon. The probability of membrane crossing decreases exponentially with the energy barrier (Eq. S1 in the Supporting Information). If a barrier for a drug crossing a membrane is higher than ∼10 kcal/mol, the statistical probability of spontaneous membrane crossing is very low within typical time scales (hundreds of nanoseconds) accessible by unbiased atomistic MD simulations. Therefore during unbiased MD simulations, the polar molecules do not usually enter freely the deeper parts of bilayer and the nonpolar molecules do not sample enough of the area of bulk water. Free energy profiles can also be calculated by biased MD simulations. Great advances have been made in this field in recent years, and numerous methods for obtaining free energy profiles have been developed, including umbrella sampling,18,19 z-constraint method,12,14,15,17,18,22,23 metadynamics,24,25 adap- tive biasing force,26,27 particle insertion,22 and others.28,29 However, although these techniques undoubtedly enhance sampling, all of them have drawbacks for estimating free energy profiles along bilayer normals. For example, in an analysis of interactions between charged and neutral forms of ibuprofen and aspirin with a dipalmitoylphosphatidylcholine (DPPC) bilayer, Boggara and Krishnamoorti18 noted that the drugs (especially charged forms) caused deformations of the lipid bilayer in z-constraint simulations. Similarly, MacCallum et al.30 observed “water defects” at the water−lipid interface in umbrella simulations applied for calculating free energy profiles of amino acids along the normal of a dioleoylphosphatidylcho- line (DOPC) bilayer. Such deformation of lipid bilayer in simulations was recently analyzed by Neale et al., who identified it as a systematic sampling error of considerable interest. Neale and co-workers stressed that this systematic sampling error complicated the convergence of free energy profiles of druglike molecules along lipid bilayer normals, especially for charged molecules. 1201 dx.doi.org/10.1021/ct2009208 | J. Chem. Theory Comput. 2012, 8, 1200−1211