Journal of Chemical Theory and Computation The minimum free energy (ΔG = 0 kcal/mol) of coumarin along the bilayer normal was considered the reference state in all profiles. All free energy profiles were calculated for one lipid leaflet and have been symmetrized for the other one. Error estimates of free energy were calculated as integrated standard deviations of the mean calculated either over the bins of 100 bootstraps generated by Bayessian bootstrap analysis by the g_wham program56 in umbrella sampling or over the force distribution in 0.1 nm-spaced positions along the z-axis in constraint simulations. For comparison with experimental data (see later) free energy profiles of 7-acetoxy-4-methylcoumarin and 7-acetox- ycoumarin along the normal of a DMPC bilayer were calculated (using a DMPC bilayer structure with 128 DMPC molecules and 3,655 water molecules taken from Tieleman’s Web site,58 after replacing 20 molecules of water by 10 Na+ and 10 Cl− ions). The lipid bilayer was equilibrated for 200 ns. 7-Acetoxy- 4-methylcoumarin topology was prepared by the simulation protocol described above for coumarin (including the use of RESP partial charges). Initial structure was generated from a 20 ns unbiased free simulation. During this simulation, 7- acetoxy-4-methylcoumarin penetrated to 1.8 nm from the middle of the bilayer and was then pulled into the bilayer for 4 ns with a pulling rate of 1 nm·ns−1 and a pulling force constant of 500 kJ·mol−1·nm−2 (119.5 kcal·mol−1·nm−2). The mild force constant was chosen to avoid water artifacts, as the molecule did not penetrate into region 3 freely during the unbiased simulation and was still in touch with water molecules. A 10 ns constraint simulation was performed with the same simulation protocol as applied in the coumarin CF simulation, but the simulation near the equilibrium position (0.8−1.7 nm) was prolonged to 15 ns per simulation bin. A free energy profile of 7-acetoxycoumarin was obtained by the same protocol as for 7-acetoxy-4-methylcoumarin, except the free simulation lasted 60 ns and during this time 7- acetoxycoumarin penetrated to 0.5 nm from the center of the lipid bilayer, then the molecule was pulled further into the bilayer for 1 ns with a pulling rate of 1 nm·ns−1 and force constant of 2,000 kJ·mol−1·nm−2 (477.9 kcal·mol−1·nm−2). The molecule freely penetrated close to the bilayer center and thus was pulled with a higher force constant than in the previous case. A free energy profile was obtained by z-constraint simulation, which yielded a very large minimum energy zone, with substantial variation in mean forces, so the simulation was prolonged to 15 ns in the interval between 1.0 and 2.0 nm and more frames were added (so the distance between simulation windows was 0.05 nm in the zone between 1.0 and 2.0 nm). ■ RESULTS Unbiased MD Simulations. Five independent unbiased simulations starting from coumarin in water, 3 μs long in total (2 × 1 μs, 600 ns, 2 × 200 ns), showed a tendency for coumarin to stay at the boundary between regions 2 and 3, as coumarin was most frequently located 1.4 ± 0.1 nm from the bilayer center (see Supporting Information, Figure S1). Once a coumarin molecule entered the bilayer during the first 10 ns of the simulation, it did not leave the lipid bilayer during the rest of simulation (Figure S1). Coumarin occurred in both leaflets because it can transverse the bilayer center spontaneously. Twelve successful (and ten unsuccessful) transitions between leaflets were observed during the 3 μs of simulations (Figure S1). The transition between both leaflets took place on a 100+ ns time scale, but the transition process itself was rapid Article and lasted several nanoseconds. The unbiased simulations also identified a metastable state of coumarin in the bilayer center (Figure 2), where coumarin stayed up to 10 ns. The transition between the bilayer center and one leaflet occurred on time scales of ps up to ns (Figure S1). The bilayer center penetration barrier calculated from the partition coefficient profile (cf. Equation 3) was 2.1 kcal/mol, but the water/lipids barrier could not be calculated, as the distribution of coumarin in water was not properly sampled (Figure S1 and Figure 2). Biased Simulations. The free energy profiles of coumarin in DOPC lipid bilayer reconstructed from four types of biased (UP, CP, UF, and CF) simulations showed similar trends (Figure 2). Typically, the free energy dropped as coumarin entered region 1 (cf. Figure 1). As it moved deeper into the bilayer, the free energy decreased and the global free energy minimum was reached at the border between regions 2 and 3. When coumarin moved deeper into the bilayer center, the free energy rose. A small local minimum was located in the bilayer center (region 4). So, one global minimum at 1.35−1.53 nm (with a thermally accessible region within 1.05−1.95 nm at 310 K − by thermally accessible region we mean an area with energy barrier of RT (0.616 kcal/mol at 310 K) from the energy minimum) and one local minimum in the middle of the lipid bilayer were common features of all free energy profiles (Figure 2). The bilayer center penetration barriers (ΔGpen) obtained from the free energy profiles fitted a narrow interval, varying between 2.6−3.3 kcal/mol (Table 1, Figure 2). The water/ Table 1. Properties Extracted from the Free Energy Profiles Calculated by Four Different Simulation Protocols with 50 ns of Biased Simulation Per Windowa area within a reach of a thermal ΔGwat ΔGpen simulation position of motion at 310 protocol minimum (nm) K (nm) (kcal/mol) (kcal/mol) UP 1.53 1.15 1.95 5.7 ± 0.3 3.2 ± 0.2 UF 1.35 1.09 1.75 6.7 ± 0.1 2.6 ± 0.1 CP 1.47 1.20 1.71 5.9 ± 0.2 2.8 ± 0.1 CF 1.49 1.10 1.80 6.4 ± 0.2 3.3 ± 0.1 aUP 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. ΔGwat and ΔGpen are water/lipid and bilayer center penetration barriers, respectively. Area within reach of a thermal motion is considered to be the area surrounded by an energy barrier of RT (0.616 kcal/mol, T = 310 K). lipids barrier (ΔGwat) fitted an interval of 5.7−6.7 kcal/mol, and the values calculated with pulling initial structures were lower than those calculated in simulations with initial structures from unbiased simulations (ΔGwat values derived from UP, CP, UF, and CF simulations were 5.7 ± 0.3, 5.9 ± 0.2, 6.7 ± 0.1, and 6.4 ± 0.2 kcal/mol, respectively; Table 1 and Figure 2). The UP free energy profile also showed a very shallow local minimum at 1.95 nm with an energy barrier of 0.5 kcal/mol (Figure 2). The free energy barrier of this minimum was higher than the free energy error bar estimated by statistical bootstrap analysis (0.1 kcal/mol), but the error seemed to be under- estimated. As the depth of the shallow minimum declined with increasing duration of simulation windows (see the following paragraph “Convergence of Biased Simulations”) and no state 1204 dx.doi.org/10.1021/ct2009208 | J. Chem. Theory Comput. 2012, 8, 1200−1211