Journal of Chemical Theory and Computation As this systematic error may also originate during the generation of starting structure sets used for biased MD simulations, this issue has been addressed in several previous studies using various strategies. Neale et al.13 used an inflategro procedure,31 in which (briefly) a pre-equilibrated bilayer is expanded, a molecule of interest is inserted, and the bilayer is compressed and re-equilibrated. Boggara and Krishnamoorti18 inserted molecules into a lipid bilayer manually using the VMD visualization program.32 Other approaches for generating starting structures have included growing the molecule inside the bilayer from zero size,14 pulling simulation,5 snapshots from unbiased simulation,33 and estimation of reaction coordinates using metadynamics.24 However, no methodological analyses of this problem have been previously published. Here, based on an examination of the embedding of a nonpolar druglike molecule (coumarin; 1,2-benzopyrone) in DOPC bilayer, we show that systematic sampling error is difficult to avoid, but it can be reduced by using appropriate biased simulation and initial structure set generation method. Coumarin naturally occurs in diverse plants, including tonka beans (Dipteryx odorata), vanilla grass (Anthoxanthum odor- atum), sweet woodruff (Galium odoratum), sweet clover (Meliotus L.), sweet grass (Hierochloe odorata), and cassia cinnamon (Cinnamomum aromaticum). It is absorbed by humans both orally from food and through the skin from perfumes. Further, it is a valuable test substance because it is a small, planar, rather rigid, nonpolar (logPoct/wat 1.3934), biologically significant druglike molecule, as its skeleton can be recognized in many drugs (e.g., the anticoagulant warfarin and antispasmodic/insecticide hymecromone35) and other biologically active compounds (e.g., scopoletin). Therefore, coumarin is an ideal model for assessing the quality of various methods for calculating ΔG profiles of small low-polar druglike molecules along normals of lipid bilayers. DOPC bilayer has been previously used as a model of endoplasmic reticulum membrane, in which coumarin is metabolized by membrane- anchored Cytochrome P450 2A6.36,37 The main aim of the study was to identify the mean position of coumarin in DOPC bilayer from calculations of free energy profiles using different biased MD simulations and different sets of initial structures. We also discuss convergence, advantages and disadvantages of z-constraint, and umbrella sampling methods using starting structures obtained by pulling and unbiased simulations. We focus on the systematic bias caused by choice of the initial structure set and the possibilities of avoiding this bias. The effect of choice of partial charges is also analyzed and results of biased and unbiased MD simulations (3 μs in total) are compared. Finally, a robust simulation protocol for obtaining a convergent free energy profile along a bilayer normal is suggested and tested against available experimental data for two coumarin derivatives embedded in dimyristoyl- phosphatidylcholine (DMPC) bilayer. ■ METHODS The structure and topology of coumarin (1,2-benzopyrone; CAS number 91-64-5) was generated by the PRODRG2 Beta server38 using the GROMOS 53a6 forcefield.39 However, partial charges assigned by PRODRG have been found to lead to unrealistic partitioning between water and cyclohexane phases.40 As the partial charges used can introduce another systematic error into free energy calculations, we addressed this problem by also using Mulliken partial charges and restrained fit of electrostatic potential (RESP) partial charges. The RESP Article partial charges were successfully adopted by the second generation of AMBER family force fields. The electrostatic potential (ESP) and ESP partial charges were calculated by applying B3LYP/cc-pVDZ method to coumarin geometry optimized at the same level of theory in Gaussian 03.41 RESP fit42 was implemented by Antechamber from the AMBER 11 software package.43 Mulliken partial charges, that were adopted by Berger lipid force field,44 were calculated at the HF/6-31G* level in gas phase. Hereafter, all mentioned coumarin charges are RESP charges, except those explicitly named as PRODRG or Mulliken charges. The lipid bilayer, as prepared and equilibrated by Siu et al.,45 contained 128 DOPC molecules, 64 in each leaflet, with a structure generated by the Lipidbook server.46 The bilayer was oriented perpendicularly to the z-axis of the simulation box and equilibrated for another 10 ns by a free MD simulation. Water and salt (NaCl) were added to give a physiological concentration, of 0.154 M, of salt in the aqueous phase (excluding the lipid bilayer from the volume calculation). The equilibrated box contained 5,188 molecules of Flexible Simple Point Charge (SPC) water,47 19 Na+ and 19 Cl− ions. The equilibrated surface area per lipid was 0.638 nm2, and the start of the z-axis was set in the middle of the bilayer. The GROMACS 4.0.7 package48 and united atom Berger lipid force field44 were used for MD simulations. The latter reduces the number of atoms in simulations, as it merges nonaromatic and nonpolar hydrogens with their carbons. This simplification likely results in higher diffusion coefficients than those observed in all-atom model simulations.49 Berger lipid force field44 uses the Mulliken partial charges calculated at the HF/6-31G* level (in gas phase).50 Simulations were taken with 2-fs integration time steps under periodic boundary conditions in all directions, with particle-mesh Ewald (PME) electro- statics,51 a van der Waals cutoff at 1 nm, bond constraints determined by the LINCS algorithm,52 V-rescale temperature coupling53 to 310 K, and Berendsen anisotropic pressure coupling54 to 1 bar with 10 ps time constant and compressibility of 4.5 × 10−5 bar−1. A coumarin molecule was placed at the top of the simulation box, a 0.5 ns MD simulation was executed to pre-equilibrate the system, then five independent MD simulations with a total time of 3 μs were generated. From the pre-equilibrated simulation two sets of starting frames for biased MD simulations were generated: one by pulling coumarin to the bilayer center and the other from unbiased MD simulation. The first set of starting frames for biased MD simulations was obtained by pulling the center-of-mass (COM) of coumarin against that of the lipid bilayer (in its center). Coumarin was pulled along the bilayer normal (the z-axis) for 6 ns using a pulling force constant of 10,000 kJ·mol−1·nm−2 (2,390 kcal·mol−1·nm−2) and pulling rate of 1 nm·ns−1. Pulling applies a harmonic potential on molecule and moves the center of this potential with a given pull rate. Starting positions were collected as snapshots from the pulling simulation, spaced 0.1 ± 0.02 nm apart along the z-axis from the area of bulk water (4 nm from the bilayer center) to the middle of the bilayer. From the structures at one distance bin, the structure with the lowest potential energy was chosen as the starting frame for biased MD simulations at a given distance from the center of the lipid bilayer. Hereafter, constraint and umbrella simulations with initial structures generated by pulling simulations are referred to as constraint-pulling (CP) and umbrella-pulling (UP), respectively. 1202 dx.doi.org/10.1021/ct2009208 | J. Chem. Theory Comput. 2012, 8, 1200−1211