Journal of Chemical Theory and Computation Article Figure 3. Initial structures (at 1.9 nm from the bilayer center) obtained by pulling (a) and free simulation (b) show a difference in coumarin hydration. The structure generated by pulling simulations indicates that coumarin is pulled to the lipid bilayer with its solvation shell, which causes funnel-like bilayer deformation. Snapshots taken at 10 ns indicate that CP eliminates coumarin hydration more rapidly than UP. Both CF and UF simulations lead to similarly solvated coumarin structures. Carbons are colored in cyan, oxygens red, and hydrogens white. The olive and blue balls represent DOPC phosphate and nitrogen atoms, respectively. Waters surrounding coumarin are colored as red/white balls.UP and UF refer to umbrella simulations with initial structures obtained by pulling and free unbiased simulation, respectively; CP and CF refer to constraint simulation with pulling and free initial structures, respectively. corresponding to this minimum was observed in the unbiased simulation (Figure S1 in the Supporting Information), we considered this minimum to be an artifact and called it “the artificial minimum”. As such “artificial minima” may be due to systematic sampling errors, it was of considerable interest to determine the reasons for such behavior. The primary reason lay in the starting structures, which were generated by pulling coumarin along the bilayer normal in the UP and CP simulations. The pulling caused deformation of the lipid bilayer,13,18 leading to a funnel-shaped bilayer surface depression (see Supporting Information, Figure S2) induced by solvated coumarin (Figure 3a). In other words, the pulling procedure produced structures in which coumarin embedded in the bilayer was hydrated by a few water molecules. The defects, namely these which were deeper in lipid bilayer, were eliminated during biased MD simulations, as water was expelled from the hydrophobic bilayer interior and bilayer relaxed rapidly on hundreds of picoseconds time scale. The relaxation time grew with increasing distance from the bilayer center. In the area of the artificial minimum (∼1.7−2.0 nm), close to the bilayer surface, the water defects were eliminated on a tens of nanoseconds time scale. It is of considerable interest that relaxation occurred significantly more rapidly in the CP than in UP simulations, and thus the membrane deformation was eliminated more rapidly (Figure 3 and Figure 4). The lipid bilayer depression was not observed during a spontaneous embedment of coumarin in the lipid bilayer in the unbiased simulations, hence the starting structures for biased simulations based on the snapshots from the unbiased simulation were free of this artifact (Figure 3b). Convergence of Biased Simulations. The position of the global minimum converged more rapidly in biased simulations starting from free (unbiased) simulations (UF and CF) than in simulations starting from pulling (UP and CP) simulations (Figure 4). The free energy profile obtained from UP simulation depended strongly on the length of the simulation windows, and two energetically similar minima in region 2 (one at ∼1.5 and the other at ∼2.0 nm) were observed during the beginning of this simulation (Figure 4). After 10 ns the minimum at ∼1.5 nm became the global minimum and its position converged to 1.53 nm, while the energy barrier of the artificial minimum decreased to 0.5 kcal/mol. During the first 16 ns of UP simulation the area accessible by thermal motion (ΔGmin+RT) gradually widened from 0.90 nm after 5 ns to 1.10 nm, and the region accessible by thermal motion thereafter declined to 0.80 nm. ΔGwat gradually rose throughout the simulation, to a final value (at 50 ns) of 5.7 ± 0.3 kcal/mol, while ΔGpen dropped within the first 16 ns of simulation, slowly rose until 30 ns, and then fluctuated around a final value of 3.2 ± 0.2 kcal/mol (cf. Figure 4). The free energy profile obtained from CP simulation also displayed two minima initially, while the artificial minimum (at ∼2.0 nm) quickly vanished, and after ∼15 ns there was no sign of this minimum. The area within reach of thermal motion 1205 dx.doi.org/10.1021/ct2009208 | J. Chem. Theory Comput. 2012, 8, 1200−1211