Journal of Chemical Theory and Computation Article Figure 2. Density profiles of DOPC (blue dash-dotted curves) and coumarin calculated from unbiased MD simulations (3 μs in total; red dashed curves) and from the free energy profile acquired by the CF method (green dotted curve). The density profiles of coumarin obtained from both methods match each other well (upper panel − a). Free energy profiles obtained from biased simulations (lower panel − b). Both umbrella (UP − black curve and UF − red curve) and constraint (CP − dotted blue curve and CF − dotted green curve) simulations provide free energy minima positions for coumarin that overlap well with the maximum density calculated from the free simulation (cf. upper panel − a). The free energy profiles were calculated for one bilayer leaflet and were symmetrized to the other one; the density was symmetrized along the middle of the bilayer. 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. In the other approach for generating starting frames, in which free unbiased simulation (20 ns long) was applied, spontaneous coumarin penetration through the DOPC bilayer was observed. Starting frames (spaced 0.1 ± 0.02 nm apart along the z-axis with the lowest potential energy of the structure at the respective points) for umbrella and constraint simulations were chosen as described above. Simulations with these starting structures are henceforth referred to as umbrella-free (UF) and constraint-free (CF), respectively. With both sets of starting frames, umbrella sampling and constraint simulations (164 simulation windows in total) were carried out as described below for 30 ns per simulation window, except for simulation bins in the 1.0−2.5 nm z-axis region (for which simulation was prolonged for 50 ns, or up to 100 ns for CP and CF simulations, giving in the latter cases up to 6.9 μs of biased simulation in total). In umbrella sampling a harmonic potential is applied between COMs of two groups of molecules, here the drug coumarin and DOPC lipid bilayer. The distance between COMs of coumarin and DOPC was restrained by a harmonic force constant of 2,000 kJ·mol−1·nm−2 (477.9 kcal·mol−1· nm−2). The force applied on coumarin was proportional to the square of the displacement from its original position, and a free energy profile was calculated from eq 14,55 ΔG(z) = −RT ln P(z) + U(z) (1) where P(z) and U(z) are the coumarin distribution and biasing potential along the bilayer normal, respectively. The force constant and distance between simulation windows were chosen to achieve equal sampling, as the presence of regions with low sampling density increases the error of umbrella sampling. Forces acting on coumarin were analyzed, and the free energy profile was reconstructed by the weighted histogram analysis method (WHAM)20 using the g_wham program.56 We also calculated a free energy profile from constraint simulations.14,15,18 In this approach, the distance between COMs of the drug and lipid bilayer was constrained, and the constraint force was monitored. Free energy was then calculated from eq 24,16,17,57 z ΔG(z) = −∫ outside ⟨F(⃗z′)⟩t dz′ where the mean force applied on the molecule ⟨F⃗(z′)⟩t in certain bilayer depth z′ is integrated along the bilayer normal axis beginning in water until the certain bilayer depth z. Part of the free energy profile was also calculated from the partitioning displayed in an unbiased simulation using eq 315 ΔG(z) = −RT ln K(z) (2) (3) where K(z) is a partition coefficient estimated for a 0.02 nm bin in bilayer depth z, symmetrized for both leaflets. The partition coefficient is calculated from average mass density of coumarin in certain bin and by normalizing this density − the reference state is set to have K(z) = 1. 1203 dx.doi.org/10.1021/ct2009208 | J. Chem. Theory Comput. 2012, 8, 1200−1211