Article pubs.acs.org/JCTC Convergence of Free Energy Profile of Coumarin in Lipid Bilayer Market́a Paloncyóva ,́ † Karel Berka,*,† and Michal Otyepka*,† †Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacky University, tr. 17 listopadu 12, 771 46, Olomouc, Czech Republic *S Supporting Information ABSTRACT: Atomistic molecular dynamics (MD) simulations of druglike molecules embedded in lipid bilayers are of considerable interest as models for drug penetration and positioning in biological membranes. Here we analyze partitioning of coumarin in dioleoylphosphatidylcholine (DOPC) bilayer, based on both multiple, unbiased 3 μs MD simulations (total length) and free energy profiles along the bilayer normal calculated by biased MD simulations (∼7 μs in total). The convergences in time of free energy profiles calculated by both umbrella sampling and z-constraint techniques are thoroughly analyzed. Two sets of starting structures are also considered, one from unbiased MD simulation and the other from “pulling” coumarin along the bilayer normal. The structures obtained by pulling simulation contain water defects on the lipid bilayer surface, while those acquired from unbiased simulation have no membrane defects. The free energy profiles converge more rapidly when starting frames from unbiased simulations are used. In addition, z-constraint simulation leads to more rapid convergence than umbrella sampling, due to quicker relaxation of membrane defects. Furthermore, we show that the choice of RESP, PRODRG, or Mulliken charges considerably affects the resulting free energy profile of our model drug along the bilayer normal. We recommend using z- constraint biased MD simulations based on starting geometries acquired from unbiased MD simulations for efficient calculation of convergent free energy profiles of druglike molecules along bilayer normals. The calculation of free energy profile should start with an unbiased simulation, though the polar molecules might need a slow pulling afterward. Results obtained with the recommended simulation protocol agree well with available experimental data for two coumarin derivatives. ■ INTRODUCTION Passive transport of drugs through membranes is the main process limiting their penetration into cells (in the absence of a specific active transporter) and thus a key step in their administration to the bodies of humans (and animals). Diffusion through membrane and partitioning between water and membrane phases are the key properties for this passive transport affecting kinetics and thermodynamics of permeation process,1,2 respectively. Further, the equilibrium position of specific drugs in target membranes also affects their metabolism and transport (both active and passive).3−5 The composition of biological membranes is complex and diverse, varying substantially among the outer and inner leaflets of both organelles and organs.6,7 They consist of proteins and lipids, in approximately equal mass proportions.8 While proteins are responsible for active transport and signaling, lipids pose the main barrier to passive membrane transport. The most important membrane for drug administration is the plasma membrane, through which drugs must penetrate to reach the internal milieu of target cells. However, mitochondrial and endoplasmic reticulum membranes are also involved in drug metabolism because they accommodate various drug- metabolizing enzymes (e.g., cytochromes P450 and UDP- glucuronosyltransferases).5,9,10 The most abundant lipids in mammalian membranes are phosphatidylcholines (PC), although phosphatidylserines, phophatidylethanolamines, sphingomyelins, and cholesterol are also present,11 thus PC bilayers are commonly used as simple membrane models. However, it must be remembered that in vivo membranes are much more complex, so results obtained using such simple models should be interpreted cautiously. Several structural frameworks of lipid bilayers have been proposed, including the four-region model of Marrink and Berendsen12 and others presented by Neale et al.13 and Orsi et al.14 The four-region model, applied in the study presented here, describes the physicochemical properties and densities of lipids in the following four regions along a bilayer’s normal axis (Figure 1): (i) The low headgroup density region (hereafter region 1), a polar zone with similar transport conditions to water, from the point where head groups are first encountered (at minimal density) and ending where the densities of head groups and water are comparable. (ii) The highly structured high headgroup density region (region 2), from the point where region 1 ends to the point closer to the bilayer center where the density of water decreases to below 1% and bulklike water disappears. Strong Coulombic interactions between polar groups keep polar molecules in the first two regions.15 (iii) The high density of acyl chains region (region 3) is hydrophobic. Double bonds of unsaturated lipids are typically localized in this region. (iv) The fourth, low density of acyl chains region (region 4), resides in the middle of the bilayer and terminal methyl groups are primarily located in this region. Here, movement of all molecules is faster due to its low density. The two hydrophobic acyl chain regions are Received: December 22, 2011 Published: February 24, 2012 © 2012 American Chemical Society 1200 dx.doi.org/10.1021/ct2009208 | J. Chem. Theory Comput. 2012, 8, 1200−1211