Journal of Chemical Theory and Computation reasonably well (1.46 and 1.31 nm, respectively, see Figure 6 down). In summary, the MD results for both coumarin derivatives agreed reasonably well with the experimental results, notably 7- acetoxy-4-methylcoumarin was located deeper in the bilayer than 7-acetoxycoumarin, and the carbon atoms’ positions calculated from simulations matched those deduced from experiments. ■ DISCUSSION Coumarin Preferentially Stays in Bilayer Regions 2 and 3 in Unbiased Simulations. During the five independent unbiased simulations (3 μs long in total) coumarin preferred the lipid bilayer phase rather than the aqueous phase, because it quickly (within <10 ns) entered the lipid bilayer and remained there for the rest of the simulation time (Figure S1 in the Supporting Information). The preferentially occupied position was at 1.4 ± 0.1 nm, at the border of regions 2 and 3, and the molecule was oriented mainly with its oxygens pointing toward the water phase (data not shown). In addition, coumarin penetrated the lipid bilayer spontaneously, i.e., moved from one leaflet to the other, remaining at the preferentially occupied positions in both leaflets for several hundreds of nanoseconds between brief (a few ns) visits to the lipid bilayer center (Figure S1 in the Supporting Information). Similarly, the transition movements were quite rapid, generally occurring within several nanoseconds. Some key penetration properties were identified from the unbiased simulations, namely positions of local and global free energy minima and qualitative estimates of the height of energy barriers (Figure 2). The water/lipids barrier, ΔGwat, seemed to be higher than the bilayer center penetration barrier, ΔGpen. The number of penetration events allowed us to roughly estimate the absolute value of the bilayer center penetration barrier, at 2.1 kcal/mol (eq 3). However, the estimated ΔGpen value should be interpreted with care, due to the limited sampling as only a small number of transitions between the minima were observed, and ΔGwat could not be calculated as the water phase was not sampled adequately (Figure 2). Nevertheless, the spontaneous embedding of coumarin in DOPC bilayer strongly indicates that this is a barrierless process and that coumarin prefers the bilayer phase, in accordance with expectations based on coumarin’s logPoct/wat value. Free Energy Profiles Obtained by Biased and Unbiased Simulation Agree. The final free energy profiles obtained by all simulation protocols (UP, CP, UF, and CF) were in accord (Figure 2), but those obtained from the unbiased simulations provided more accurate information. The global energy minimum was found at 1.44 ± 0.09 nm, while a local energy minimum was localized in the membrane center. The presence of a local energy minimum in the lipid bilayer center agrees with previous findings presented by Bemporad et al.,23 of such local minima for some other small solutes, e.g. water and acetamide. In our case, the bilayer center penetration barrier (ΔGpen) of coumarin spanned 2.6−3.3 kcal/mol (Table 1, Figure 2), close to the ΔGpen estimated roughly from the unbiased simulation (2.1 kcal/mol). In contrast, the water/ lipids barrier (ΔGwat) varied significantly with time and method used (see below). The estimated ΔGwat from CF simulation (which is taken as reference, as constraint biasing eliminates possible artificial errors quicker) was 6.4 ± 0.2 kcal/mol (Figure 2). The free energy profiles therefore confirmed that Article coumarin more readily penetrates the bilayer than escapes to the water phase. Partial Charges − The Force Field Issue. Neither accurate unbiased simulation, nor accurate free energy profile calculation, is possible without a careful choice of force field. The Berger force field (using Mulliken partial charges calculated at HF/6-31G* level (in gas phase)50) used for lipids39,44 was tested and shown to provide area per lipid and volume per lipid values that correspond well with experimental values.44 Further, as coumarin parameters were not available in standard data sets for lipid simulations, they had to be acquired separately. Generally, atom types and corresponding parame- ters can be adopted for a nonlipid molecule from the standard data sets quite safely, but the set of partial charges had to be carefully considered, as it may introduce a serious systematic sampling error in lipid bilayer-guest molecule simulations. We addressed this issue by using three sets of partial charges (Figure 5, Table S1 in the Supporting Information): one generated by the PRODRG server, one assigned by Mulliken population analysis, and the third generated by applying the RESP procedure in B3LYP/cc-pVDZ calculations of electro- static potential in gas phase. Generally, increasing the dipole moment of a molecule (by use of PRODRG or Mulliken partial charges) resulted in a lower ΔGwat and higher ΔGpen, in accordance with expectations, given the higher polarity of coumarin bearing PRODRG or Mulliken charges in comparison with RESP partial charges (Figure 5). With only these profiles it would be difficult to decide which partial charges provided more reliable results. However, PRODRG charges led to overestimation of the dipole moment of coumarin and (as mentioned above) partial charges assigned by the PRODRG server lead to unrealistically strong partitioning in water in cyclohexane/water systems as found by Lemkul et al.40 The latter finding agrees with the trends observed in our lipid bilayer simulations. In summary, RESP charges seem to provide more accurate models for simulations of lipid bilayer-guest molecule systems than PRODRG charges (although whether the RESP charges should ideally be based on gas phase or solvent-polarized ESP, and if they can be robustly combined with the Berger force field for lipids, remains to be determined). Convergence of Free Energy Profiles − The Artificial Minimum Issue. We have shown here that the convergence of free energy profiles was significantly influenced by the generation of initial structures when followed by the biasing method. The biased simulations starting from the pulling simulations (UP and CP) suffered from bilayer deformation induced by pulling coumarin from the water phase toward the bilayer center (Figure 3). Similar bilayer deformations have been repeatedly previously observed13,18,30 and identified as a systematic sampling artifact in biased lipid bilayer simulations. For example, Neale et al.13 observed bilayer deformations when a charged molecule was embedded in the bilayer. We observed a funnel-shape bilayer surface depression (Figure S2 in the Supporting Information), caused by water hydrating the polar parts of coumarin penetrating the lipid bilayer more deeply and thereby exacerbating bilayer deformation during the pulling simulations. The bilayer deformation caused an artificial minimum (∼2.0 nm) in the free energy profiles in region 2 (Figure 2), whereas in unbiased simulation coumarin never stayed longer in this position, and its behavior showed no sign of reaching a local energy minimum. 1208 dx.doi.org/10.1021/ct2009208 | J. Chem. Theory Comput. 2012, 8, 1200−1211