were used. The used methods for free energy profile calculation (umbrella and constraint simulation) are quite equivalent when applied on an error-free set of starting structures. However, if the membrane defects are present, the z-constraint simulation leads to more rapid convergence than umbrella sampling. In summary, for efficient calculation of convergent free energy profiles of druglike molecules along bilayer normals, we recommend using z-constraint biased MD simulations based on as much starting geometries acquired from unbiased MD simulations as possible, otherwise when pulling simulation is employed, the biased simulation might need far longer time to reach convergence. ■ ASSOCIATED CONTENT *S Supporting Information Details of the recommended simulation protocol, used partial charges, four supplemental figures, and one equation. Table S1: Partial charges on coumarin atoms used in the paper. Figure S1: Position of coumarin in unbiased simulations. Figure S2: Funnel shape water defect caused by pulling coumarin into the membrane. Figure S3: Effect of window spacing in constraint and umbrella simulations. Figure S4: Atom labels as used in Table S1. This material is available free of charge via the Internet at http://pubs.acs.org. AUTHOR INFORMATION Corresponding Author *E-mail: karel.berka@upol.cz (K.B.), michal.otyepka@upol.cz (M.O.). Notes The authors declare no competing financial interest. Journal of Chemical Theory and Computation This artificial minimum was most profound when short simulation times (<5 ns) for each sampling bin were applied, and it slowly disappeared when the simulation time was prolonged (Figure 4). The main reason for the slow convergence and need for longer simulation times was the slow coumarin water shell elimination in region 2. Furthermore, the presence of the artificial minimum led to underestimation of ΔGwat in simulations using the initial structure set generated by pulling simulation (CP and UP). While ΔGwat values obtained by CF and UF simulations seemed to reach convergence, they did not reach convergence in UP and CP simulations during 50 ns of simulation (or even during 100 ns of CP simulation), although both UP and CP yielded ΔGwat values close to those obtained from CF and UF simulations. In this respect, constraint biasing was more effective, as the artificial minimum was eliminated within 15 ns per bin (in CP), while there were signs of the artificial minimum in UP simulation even after 50 ns per bin (Figure 4). Even longer times may be needed in simulations of polar or charged residues, as previously shown by MacCallum et al.30 and Neale et al.,13 who found that 80 to 205 ns per bin may be required to achieve convergence in umbrella simulations with charged solutes. In contrast, for nonpolar solutes Neale et al. achieved convergence more rapidly (in some cases after 20 ns per bin)., These water artifacts seem to be present when nonequilibrated initial structures are used for biased simulation. The higher efficiency of constraint over umbrella biasing is also consistent with the recent observation by Gunsteren et al.,61 that constraint-biased simulation using force averaging is the most effective method for calculating potential of mean force with respect to a distance from a given reference point. Convergence of the free energy profiles was clearly achieved more rapidly when using starting structures acquired from unbiased (UF and CF) simulations in comparison with the pulling simulation (UP and CP), since in cases of UF and CF simulation the free energy profiles changed only marginally with increases in the length of the simulation bins (Figure 4). Therefore we recommend starting the calculation of free energy profile with an unbiased simulation for all molecules, in a case of more polar molecules a slow pulling simulation (pulling force constant <500 kJ·mol‑1·nm‑2 (119.5 kcal·mol‑1·nm‑2) and a pulling rate <1 nm·ns‑1) from the deepest position in the lipid bilayer should follow. Thus, this approach was used for comparing the calculated results with experimental data, and the calculated positions of 7-acetoxy-4-methylcoumarin and 7- acetoxycoumarin in DMPC bilayer agreed well with positions derived from NMR experiments (Figure 6). In summary, whenever possible biased simulations should start from geometries acquired from unbiased MD simulations, and constraint biasing is the recommended and quickly converging method. ■ CONCLUSION The convergence in time of free energy profiles of coumarin along a DOPC bilayer normal, calculated by both umbrella sampling and z-constraint techniques, was thoroughly analyzed. Two sets of starting structures were also considered: one based on unbiased MD simulation and the other on “pulling” coumarin along the bilayer normal. Water defects on the lipid bilayer surface were identified in the structures obtained by pulling simulation but not in structures acquired from unbiased simulation. Consequently, the free energy profiles converged more rapidly when starting frames from unbiased simulations Article ■ ■ ■ ACKNOWLEDGMENTS Support through GACR grants 303/09/1001 and P208/12/ G016 and Student Project PrF_2011_020 provided by Palacky University is gratefully acknowledged. This work was also supported by the Operational Program Research and Develop- ment for Innovations − European Regional Development Fund (CZ.1.05/2.1.00/03.0058) and European Social Fund (CZ.1.07/2.3.00/20.0017). REFERENCES (1) Orsi, M.; Essex, J. W. Passive Permeation Across Lipid Bilayers: a Literature Review. In Molecular Simulations and Biomembranes, 1st ed.; Sansom, M. S. P., Biggin, P. C., Eds.; Royal Society of Chemistry: 2010; pp 76−90. (2) Balaz, S. Modeling Kinetics of Subcellular Disposition of Chemicals. Chem. Rev. 2009, 109, 1793−1899. (3) Afri, M.; Gottlieb, H. E.; Frimer, A. A. Superoxide Organic Chemistry within the Liposomal Bilayer, part II: a Correlation between Location and Chemistry. Free Radical Biol. Med. 2002, 32, 605−618. (4) Xiang, T.-X.; Anderson, B. D. Liposomal Drug Transport: a Molecular Perspective from Molecular Dynamics Simulations in Lipid Bilayers. Adv. Drug Delivery Rev. 2006, 58, 1357−1378. (5) Berka, K.; Hendrychova,́ T.; Anzenbacher, P.; Otyepka, M. Membrane Position of Ibuprofen Agrees with Suggested Access Path Entrance to Cytochrome P450 2C9 Active Site. J. Phys. Chem. A 2011, 115, 11248−11255. (6) Alberts, B.; Johnson, A.; Lewis, J.; Raff, M.; Roberts, K.; Walter, P. Molecular Biology of the Cell, 4th ed.; Garland Science: New York, 2002. (7) Cooper, G. M. The Cell: A Molecular Approach, 2nd ed.; Sinauer Associates: Sunderland MA, 2000. 1209 dx.doi.org/10.1021/ct2009208 | J. Chem. Theory Comput. 2012, 8, 1200−1211