passive transport of drugs through membranes 
 is the main process 
 limiting their penetration into cells ( in the absence of a specific 
 active transporter ) and thus a key step in their administration to 
 the bodies of humans ( and animals ) . diffusion through membrane and 
 partitioning between water and membrane phases are the key properties 
 for this passive transport affecting kinetics and thermodynamics of 
 permeation process , respectively .
further , the equilibrium 
 position of specific drugs in target membranes also affects their 
 metabolism and transport ( both active and passive ) .
the composition 
 of biological membranes is complex and diverse , 
 varying substantially among the outer and inner leaflets of both organelles 
 and organs .
they consist of proteins and lipids , in 
 approximately equal mass proportions . while 
 proteins are responsible for active transport and signaling , lipids 
 pose the main barrier to passive membrane transport .
the most important 
 membrane for drug administration is the plasma membrane , through which 
 drugs must penetrate to reach the internal milieu of target cells . 

however , mitochondrial and endoplasmic reticulum membranes are also 
 involved in drug metabolism because they accommodate various drug - metabolizing 
 enzymes ( e.g. , cytochromes p450 and udp - glucuronosyltransferases ) .
the most abundant lipids in mammalian 
 membranes are phosphatidylcholines ( pc ) , although phosphatidylserines , 
 phophatidylethanolamines , sphingomyelins , and cholesterol are also 
 present , thus pc bilayers are commonly 
 used as simple membrane models . however , it must be remembered that 
 in vivo membranes are much more complex , so results obtained using 
 such simple models should be interpreted cautiously .
several 
 structural frameworks of lipid bilayers have been proposed , 
 including the four - region model of marrink and berendsen and others presented by neale et al . and orsi et al .
the four - region model , applied in the study presented here , describes 
 the physicochemical properties and densities of lipids in the following 
 four regions along a bilayer s normal axis ( figure 1):(i)the low headgroup density region ( hereafter region 1 ) , a polar zone with similar transport 
 conditions to water , from the point where head groups are first encountered 
 ( at minimal density ) and ending where the densities of head groups 
 and water are comparable.(ii)the highly structured high 
 headgroup density region ( region 2 ) , from the point where 
 region 1 ends to the point closer to the bilayer center where the 
 density of water decreases to below 1% and bulklike water disappears . 

strong coulombic interactions between polar groups keep polar molecules 
 in the first two regions.(iii)the high density of acyl 
 chains region ( region 3 ) is hydrophobic .
double bonds of 
 unsaturated lipids are typically localized in this region.(iv)the fourth , low 
 density of 
 acyl chains region ( region 4 ) , resides in the middle of the 
 bilayer and terminal methyl groups are primarily located in this region . 
 here ,
the two hydrophobic acyl chain regions are believed to form the main 
 barrier for most druglike molecules , which are often water - soluble .
the low headgroup density region ( hereafter region 1 ) , a polar zone with similar transport 
 conditions to water , from the point where head groups are first encountered 
 ( at minimal density ) and ending where the densities of head groups 
 and water are comparable .
the highly structured high 
 headgroup density region ( region 2 ) , from the point where 
 region 1 ends to the point closer to the bilayer center where the 
 density of water decreases to below 1% and bulklike water disappears . 

strong coulombic interactions between polar groups keep polar molecules 
 in the first two regions .
the fourth , low 
 density of 
 acyl chains region ( region 4 ) , resides in the middle of the 
 bilayer and terminal methyl groups are primarily located in this region . 
 here ,
the two hydrophobic acyl chain regions are believed to form the main 
 barrier for most druglike molecules , which are often water - soluble .
upper left panel ( a ) : density profile of dopc bilayer along the 
 normal to the lipid bilayer plane showing densities of the complete 
 system ( black ) , water ( blue ) , dopc ( red ) , phosphates ( magenta ) , cholines 
 ( cyan ) , carbonyls ( green ) , and terminal carbons ( dark blue ) .
lower 
 left panel ( b ) : free energy profile of coumarin along the dopc bilayer 
 normal calculated from constraint simulation with initial structures 
 obtained by free simulation ( cf ) .
the calculated bilayer center penetration 
 barrier , g , and water / lipids 
 barrier , g , are labeled .
the 
 free energy profile was calculated for one bilayer leaflet and was 
 symmetrized to the other one , the densities of the system were symmetrized 
 along the middle of the bilayer .
the vertical bins labeled by numbers 
 denote four bilayer regions : 1  low density of head groups 
 ( 2.22.9 nm ) , 2  high density of head groups ( 1.452.2 
 nm ) , 3  high density of acyl chains ( 0.51.45 nm ) , 
 and 4  low density of acyl chains ( 00.5 nm ) .
right 
 panel ( c ) : structure of dopc bilayer , together with snapshots of coumarin 
 initial structures .
molecular dynamics ( md ) simulations can be used 
 to estimate the 
 equilibrium position of a drug in a lipid bilayer , its partition coefficients 
 and diffusion coefficients simultaneously at subpicosecond and atomic 
 resolution . the partitioning and mean position can be well described 
 by a free energy ( g ) profile along the normal 
 to the lipid bilayer , also known as a potential of mean force 
 ( pmf ) . in principle
, such free energy profiles can be calculated from 
 partitioning values obtained by long unbiased simulation
. however , 
 this approach only provides reliable free energy profiles when all 
 states along the g profile are thoroughly 
 sampled .
this is challenging for unbiased md simulations , because 
 they usually do not sample adequately at the available simulation 
 time scales .
the sampling problem is based on the dependence of probability 
 of drug crossing a membrane on the energy barrier for this phenomenon . 

if a barrier for a drug crossing a membrane is higher than 10 
 kcal / mol , the statistical probability of spontaneous membrane crossing 
 is very low within typical time scales ( hundreds of nanoseconds ) accessible 
 by unbiased atomistic md simulations .
therefore during unbiased md 
 simulations , the polar molecules do not usually enter freely the deeper 
 parts of bilayer and the nonpolar molecules do not sample enough of 
 the area of bulk water .
great advances have been made in this field in recent years , and numerous 
 methods for obtaining free energy profiles have been developed , including 
 umbrella sampling , z - constraint method , metadynamics , adaptive biasing force , particle insertion , and others .
however , although these techniques undoubtedly enhance sampling , 
 all of them have drawbacks for estimating free energy profiles along 
 bilayer normals .
for example , in an analysis of interactions between 
 charged and neutral forms of ibuprofen and aspirin with a dipalmitoylphosphatidylcholine 
 ( dppc ) bilayer , boggara and krishnamoorti noted that the drugs ( especially charged forms ) caused deformations 
 of the lipid bilayer in z - constraint simulations .
similarly , maccallum 
 et al . observed  water defects  
 at the water
 lipid interface in umbrella simulations applied 
 for calculating free energy profiles of amino acids along the normal 
 of a dioleoylphosphatidylcholine ( dopc ) bilayer .
such deformation 
 of lipid bilayer in simulations was recently analyzed by neale et 
 al .
neale and co - workers stressed that this systematic sampling 
 error complicated the convergence of free energy profiles of druglike 
 molecules along lipid bilayer normals , especially for charged molecules .
as this systematic error may also originate during the generation 
 of starting structure sets used for biased md simulations , this issue 
 has been addressed in several previous studies using various strategies . 

neale et al . used an inflategro procedure , in which ( briefly ) a pre - equilibrated bilayer 
 is expanded , a molecule of interest is inserted , and the bilayer is 
 compressed and re - equilibrated .
boggara and krishnamoorti inserted molecules into a lipid bilayer manually 
 using the vmd visualization program . other 
 approaches for generating starting structures have included growing 
 the molecule inside the bilayer from zero size , pulling simulation , snapshots 
 from unbiased simulation , and estimation 
 of reaction coordinates using metadynamics .
however , no methodological analyses of this problem have been previously 
 published . here ,
based on an examination of the embedding of 
 a nonpolar druglike 
 molecule ( coumarin ; 1,2-benzopyrone ) in dopc bilayer , we show that 
 systematic sampling error is difficult to avoid , but it can be reduced 
 by using appropriate biased simulation and initial structure set generation 
 method .
coumarin naturally occurs in diverse plants , including tonka 
 beans ( dipteryx odorata ) , vanilla grass ( anthoxanthum odoratum ) , sweet woodruff ( galium odoratum ) , sweet clover ( meliotus l. ) , sweet grass ( hierochloe odorata ) , and cassia cinnamon ( cinnamomum 
 aromaticum ) .
it is absorbed by humans both orally from food 
 and through the skin from perfumes .
further , it is a valuable test 
 substance because it is a small , planar , rather rigid , nonpolar ( logpoct / wat 1.39 ) , biologically significant 
 druglike molecule , as its skeleton can be recognized in many drugs 
 ( e.g. , the anticoagulant warfarin and antispasmodic / insecticide hymecromone ) and other biologically active compounds ( e.g. , 
 scopoletin )
. therefore , coumarin is an ideal model for assessing the 
 quality of various methods for calculating g profiles of small low - polar druglike molecules along normals of 
 lipid bilayers .
dopc bilayer has been previously used as a model of 
 endoplasmic reticulum membrane , in which coumarin is metabolized by 
 membrane - anchored cytochrome p450 2a6 .
the 
 main aim of the study was to identify the mean position of 
 coumarin in dopc bilayer from calculations of free energy profiles 
 using different biased md simulations and different sets of initial 
 structures .
we also discuss convergence , advantages and disadvantages 
 of z - constraint , and umbrella sampling methods using starting structures 
 obtained by pulling and unbiased simulations .
we focus on the systematic 
 bias caused by choice of the initial structure set and the possibilities 
 of avoiding this bias .
the effect of choice of partial charges is 
 also analyzed and results of biased and unbiased md simulations ( 3 
 s in total ) are compared .
finally , a robust simulation protocol 
 for obtaining a convergent free energy profile along a bilayer normal 
 is suggested and tested against available experimental data for two 
 coumarin derivatives embedded in dimyristoylphosphatidylcholine ( dmpc ) 
 bilayer .
the structure and topology of coumarin ( 1,2-benzopyrone ; 
 cas number 
 91 - 64 - 5 ) was generated by the prodrg2 beta server using the gromos 53a6 forcefield . however
, partial charges assigned by prodrg have been found to lead 
 to unrealistic partitioning between water and cyclohexane phases .
as the partial charges used can introduce another 
 systematic error into free energy calculations , we addressed this 
 problem by also using mulliken partial charges and restrained fit 
 of electrostatic potential ( resp ) partial charges .
the resp partial 
 charges were successfully adopted by the second generation of amber 
 family force fields .
the electrostatic potential ( esp ) and esp partial 
 charges were calculated by applying b3lyp / cc - pvdz method to coumarin 
 geometry optimized at the same level of theory in gaussian 03 .
mulliken partial charges , that were adopted by berger lipid 
 force field , were calculated at the hf/6 - 31 g * 
 level in gas phase .
hereafter , all mentioned coumarin charges are 
 resp charges , except those explicitly named as prodrg or mulliken 
 charges .
the lipid bilayer , as prepared and equilibrated by 
 siu et al . , contained 128 dopc molecules , 
 64 in each leaflet , 
 with a structure generated by the lipidbook server .
the bilayer was oriented perpendicularly to the z - axis of the simulation box and equilibrated for another 
 10 ns by a free md simulation .
water and salt ( nacl ) were added to 
 give a physiological concentration , of 0.154 m , of salt in the aqueous 
 phase ( excluding the lipid bilayer from the volume calculation ) .
the 
 equilibrated box contained 5,188 molecules of flexible simple point 
 charge ( spc ) water , 19 na and 19 cl ions .
the equilibrated surface area 
 per lipid was 0.638 nm , and the start of the z - axis was set in the middle of the bilayer .
the gromacs 4.0.7 
 package and united 
 atom berger lipid force field were used 
 for md simulations .
the latter reduces the number of atoms in simulations , 
 as it merges nonaromatic and nonpolar hydrogens with their carbons . 

this simplification likely results in higher diffusion coefficients 
 than those observed in all - atom model simulations .
berger lipid force field uses 
 the mulliken partial charges calculated at the hf/6 - 31 g * level ( in 
 gas phase ) .
simulations were taken with 
 2-fs integration time steps under periodic boundary conditions in 
 all directions , with particle - mesh ewald ( pme ) electrostatics , a van der waals cutoff at 1 nm , bond constraints 
 determined by the lincs algorithm , v - rescale 
 temperature coupling to 310 k , and berendsen 
 anisotropic pressure coupling to 1 bar 
 with 10 ps time constant and compressibility of 4.5  10 bar .
a coumarin molecule 
 was placed at the top of the simulation box , 
 a 0.5 ns md simulation was executed to pre - equilibrate the system , 
 then five independent md simulations with a total time of 3 s 
 were generated . from the pre - equilibrated simulation
two sets of starting 
 frames for biased md simulations were generated : one by pulling coumarin 
 to the bilayer center and the other from unbiased md simulation .
the 
 first set of starting frames for biased md simulations was obtained 
 by pulling the center - of - mass ( com ) of coumarin against that of the 
 lipid bilayer ( in its center ) .
coumarin was pulled along the bilayer 
 normal ( the z - axis ) for 6 ns using a pulling force 
 constant of 10,000 kjmolnm ( 2,390 kcalmolnm ) and pulling rate of 1 nmns .
pulling applies 
 a harmonic potential on molecule and moves the center of this potential 
 with a given pull rate .
starting positions were collected as snapshots 
 from the pulling simulation , spaced 0.1  0.02 nm apart along 
 the z - axis from the area of bulk water ( 4 nm from 
 the bilayer center ) to the middle of the bilayer . from the structures 
 at one distance bin , the structure with the lowest potential energy 
 was chosen as the starting frame for biased md simulations at a given 
 distance from the center of the lipid bilayer .
hereafter , constraint 
 and umbrella simulations with initial structures generated by pulling 
 simulations are referred to as constraint - pulling ( cp ) and umbrella - pulling 
 ( up ) , respectively . in the other approach for generating starting 
 frames , in which 
 free unbiased simulation ( 20 ns long )
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.02.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 kjmolnm ( 477.9 kcalmolnm ) .
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 1(4,55)1where 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 ) using the g_wham program .
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 2(4,16,17,57)2where 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 3(15)3where 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 .
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 program 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-acetoxycoumarin 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 , 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 nmns and a pulling 
 force constant of 500  kjmolnm ( 119.5 kcalmolnm ) .
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.81.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 
 nmns and force constant of 2,000 kjmolnm ( 477.9 kcalmolnm ) .
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 ) .
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 ) .
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 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 ) .
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 
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 
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 .
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.351.53 nm ( with a thermally 
 accessible region within 1.051.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 ( g ) obtained from the free energy profiles fitted a narrow 
 interval , varying between 2.63.3 kcal / mol ( table 1 , figure 2 ) . the water / lipids barrier 
 ( g ) fitted an interval of 5.76.7 
 kcal / mol , and the values calculated with pulling initial structures 
 were lower than those calculated in simulations with initial structures 
 from unbiased simulations ( g 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 underestimated . 

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 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  . 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 .
g and g 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 ) . 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 , 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.72.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 ) .
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 .
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 .
convergence of free energy profiles , positions of energy 
 minima 
 and energy barriers .
the simulation windows within 1.02.5 
 nm have been simulated for 50 and 100 ns in case of up and uf simulations 
 of cp and cf simulations , respectively ; the rest of each profile is 
 calculated from 30 ns of simulation .
free energy profiles calculated 
 from short simulation times ( < 5 ns ) are biased by high error , because 
 of small data set and nonequilibrium starting structures .
free energy 
 profiles obtained by up and cp simulations ( left ) show slow elimination 
 of an artificial minimum ( 2.0 nm ) and deepening of the global 
 minimum ( 1.5 nm ) .
free energy profiles obtained from uf and 
 cf simulation are consistent in coumarin positioning and have deeper 
 energy minima than those obtained from up and cp simulations .
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 .
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 .
g gradually rose throughout the simulation , to a final value 
 ( at 50 ns ) of 5.7  0.3 kcal / mol , while g 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 gradually narrowed 
 from 1.38 nm after 5 ns of simulation to 0.51 nm after 40 ns and thereafter 
 remained constant .
g gradually 
 rose throughout the simulation , to 5.9  0.2 kcal / mol , while 
 g grew during the first 11 ns 
 of simulation , until 20 ns of simulation it gradually declined to 
 2.8  0.1 kcal / mol and then fluctuated around this value . in 
 the last 50 ns of the simulation
prolonged to 100 ns the position 
 of the energy minimum remained constant ; g fluctuated around 2.8  0.1 kcal / mol and g continued to rise , to 6.2  0.2 kcal / mol .
the position of the minimum in the uf free energy profile was almost 
 constant ( within 1.291.35 nm ) during the whole simulation 
 time , with the area accessible by thermal motion slowly widening from 
 0.44 to 0.66 nm .
g slowly decreased 
 during the first 19 ns of simulation , then very slowly increased , 
 and after 30 ns g converged 
 to 6.7  0.1 kcal / mol , while g became convergent after 20 ns , fluctuating within 2.6  0.1 
 kcal / mol .
the cf free energy profile showed a minimum position 
 within 1.291.49 
 nm , and the area of thermal motion slowly widened from 0.41 to 0.70 
 nm .
g fluctuated around 6.47.0 
 kcal / mol during the whole simulation time , while g decreased during the first 10 ns of simulation and 
 then fluctuated around 2.93.3 kcal / mol .
the prolonged simulation 
 to 100 ns showed similar trends  a free energy minimum at 
 1.29 nm , thermal motion within 0.6 nm , a constant g value of 7.0 kcal / mol after 80 ns , and g already convergent with a final value of 
 3.1  0.1 kcal / mol .
as assignment of 
 partial charges might introduce another systematic sampling error 
 into free energy calculations , we carried out 10 ns long cf simulations 
 with prodrg and mulliken charges ( assigned partial charges are listed 
 in supporting information table s1 ) , to 
 assess the extent to which the partial charges affected the free energy 
 profiles .
coumarin with partial charges assigned by prodrg bore a 
 dipole moment of 9.5 d , assignment of mulliken partial charges led 
 to 6.0 d , and resp partial charges resulted in a dipole moment of 
 4.9 d ( figure 5 ) .
the dipole moment based on 
 resp charges was close to that of coumarin in the gas phase calculated 
 by the hybrid dft method ( b3lyp / cc - pvdz ) of 4.6 d , mulliken partial 
 charges represent a compromise between the dipole moment in water 
 ( represented by continuum dielectrics with r = 78.39 ) 
 and heptane ( r = 1.92 ) , which we calculated by the 
 cpcm / b3lyp / cc - pvdz method and that resulted in 6.7 and 5.4 d , respectively . 
 considering these values , the dipole moment stemming from prodrg charges 
 seemed to be unreliably overestimated , which could systematically 
 bias free energy profiles based on prodrg charges .
the global minimum 
 of the g profile of coumarin bearing resp 
 partial charges was located at 1.29 nm ( cf with a 10 ns sampling window ) , 
 energy minimum of coumarin bearing mulliken partial charges was localized 
 at 1.20 nm , and the minimum for prodrg - charged coumarin was shifted 
 toward the bilayer / water interface , at 1.62 nm .
the global free energy 
 minimum for resp - charged coumarin was also considerably deeper than 
 for mulliken - charged or prodrg - charged coumarin ( g : 7.5 , 5.6 , and 3.3 kcal / mol , respectively ) , and the 
 bilayer center penetration barriers of the systems also differed ( g : 3.1 , 4.6 , and 10.1  kcal / mol , respectively ) 
 ( figure 5 ) . as expected , the energy cost of 
 bilayer center penetration grows with the increasing dipole moment , 
 and reversely the energy barrier between lipid bilayer and water decreases 
 with the growing dipole moment .
left panel : free energy profiles calculated 
 for coumarin with prodrg 
 ( red dotted curve ) , mulliken ( blue dashed curve ) , and resp charges 
 ( black curve ) by constraint simulation ( cf ) with initial structures 
 obtained by free simulation using 10 ns windows .
coumarin with prodrg 
 partial charges is shifted to the outer part of the lipid bilayer . 

the bilayer center penetration barriers grow and the water / lipids 
 barriers decrease with increasing dipole moment .
the right panel shows 
 that the partial charges ( mapped on the vdw surface ) calculated by 
 resp ( upper part ) and mulliken population analysis ( middle ) are spread 
 along the whole molecule , while partial charges assigned by prodrg 
 ( lower part ) are localized close to coumarin oxygens . to our knowledge , 

the precise positioning of bare coumarin in a dopc bilayer has not 
 yet been studied experimentally ; therefore , we compared the results 
 of our theoretical calculations to data obtained in experiments with 
 coumarin derivatives .
depths of several coumarin derivatives in dimyristoylphosphatidylcholine 
 ( dmpc ) bilayer have been studied in nmr investigations , in which chemical shifts of c - labeled derivatives were used to assess the polarity of 
 the surroundings of c atoms and hence estimate their 
 depth in the lipid bilayer .
results of the cited experiments indicate 
 that the mean position of 7-acetoxy-4-methylcoumarin is at the border 
 of regions 2 and 3 in dmpc lipid bilayer ( 0.7 nm from the dmpc choline 
 nitrogens , corresponding to 1.2 nm from the center of the bilayer ) ; c - labeled carbons of the derivative ( c2 and c4 , see figure 6 ) appeared to be located 0.72 and 0.70 nm from the 
 choline nitrogens , corresponding to 1.18 and 1.20 nm from the bilayer 
 center , respectively .
another coumarin derivative , 7-acetoxycoumarin , 
 was apparently located closer to the bilayer interface in region 2 , 
 with its c - labeled ( c2 and c4 ) carbons 0.44 and 0.59 
 nm from the choline nitrogen , corresponding to 1.46 and 1.31 nm from 
 the bilayer center , respectively .
we recalculated the experimental 
 positions ( originally expressed as distances from the bilayer surface ) 
 as distances from the bilayer center to facilitate direct comparison 
 with results of this study . in this recalculation ,
the distance between 
 the dmpc bilayer surface and center was set at 1.9 nm : the mean distance 
 between the bilayer center and maximum density of nitrogens ( regarded 
 as the membrane surface in the cited nmr experiments ) in corresponding 
 md simulations .
free energy profiles and structures of coumarin derivatives 
 ( 7-acetoxy-4-methylcoumarin , 
 upper panel ( a ) , and 7-acetoxycoumarin , lower panel ( b ) ) along a dmpc 
 lipid bilayer normal calculated from constraint simulation with initial 
 structures obtained by free simulation ( cf ) .
nmr - observed positions 
 of c - labeled carbons ( c2 and c4 ) are displayed as red 
 and green circles , respectively , and the positions of c2 and c4 carbons 
 calculated from simulation are depicted as green and red curves , respectively . 

the positions of marked carbons of both coumarin derivatives are in 
 good agreement with the positions observed by nmr . for the comparison we employed the most effective 
 simulation protocol 
 of those considered here to calculate the free energy profiles of 
 the coumarin derivatives described above , briefly comprising unbiased 
 simulation followed by constraint simulation with 10 ns simulation 
 bins ( see methods for details ) .
the profile 
 obtained for 7-acetoxy-4-methylcoumarin indicated the free energy 
 minimum position of its com to be 1.0 nm from the center of the dmpc 
 bilayer ( figure 6 ) , with a thermally accessible 
 region between 0.8 and 1.3 nm .
the thermally accessible region estimated 
 from the simulation ( 0.81.3 nm ) matched that acquired from 
 nmr experiments , where 7-acetoxy-4-methylcoumarin was located 1.2 
  0.1 nm from the bilayer center .
in addition , the positions 
 of its c2 and c4 carbons calculated from simulations ( 1.25 and 1.21 
 nm , respectively ) agreed well with those estimated from experiments 
 ( 1.18 and 1.20 nm , respectively ) , although the c4 carbon seems to 
 flip - flop between two positions in the bilayer ( the other at 0.76 
 nm , with ca .
the free energy minimum position of the com of 7-acetoxycoumarin 
 in the simulation was located 1.2 nm from the bilayer center , with 
 a thermally accessible region between 0.75 and 1.35 nm .
the simulated 
 distances of the c2 and c4 carbons from the bilayer center at this 
 point ( 1.34 and 1.19 nm , respectively ) again matched those obtained 
 from the nmr data 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 .
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 , g , seemed to be higher 
 than the bilayer center penetration barrier , g .
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 
 g value should be interpreted 
 with care , due to the limited sampling as only a small number of transitions 
 between the minima were observed , and g 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 .
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 . , of such local minima for some 
 other small solutes , e.g. water and acetamide . in our case ,
the bilayer 
 center penetration barrier ( g ) of coumarin spanned 2.63.3 kcal / mol ( table 1 , figure 2 ) , close to the g estimated roughly from the unbiased simulation 
 ( 2.1 kcal / mol ) .
in contrast , the water / lipids barrier ( g ) varied significantly with time and method 
 used ( see below ) . the estimated g 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 coumarin more readily penetrates the bilayer 
 than escapes to the water phase .
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 - 31 g * level 
 ( in gas phase ) ) used for lipids was tested and shown to provide area per lipid and volume per lipid 
 values that correspond well with experimental values .
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 parameters 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 electrostatic potential in gas phase .
generally , increasing the 
 dipole moment of a molecule ( by use of prodrg or mulliken partial 
 charges ) resulted in a lower g and higher g , 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 .
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 ) .
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 observed and identified as a systematic 
 sampling artifact in biased lipid bilayer simulations .
for example , 
 neale et al . 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 .
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 g in simulations using the initial structure set generated by pulling 
 simulation ( cp and up ) . while g 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 g 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 . and neale et 
 al . , who found that 80 to 205 ns per 
 bin may be required to achieve convergence in umbrella simulations 
 with charged solutes .
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 . , 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 kjmolnm ( 119.5 kcalmolnm ) and a pulling rate < 1 nmns ) 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 .
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 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 .
 atomistic molecular dynamics ( md ) simulations of druglike 
 molecules 
 embedded in lipid bilayers are of considerable interest as models 
 for drug penetration and positioning in biological membranes . here 
 
 we analyze partitioning of coumarin in dioleoylphosphatidylcholine 
 ( dopc ) bilayer , based on both multiple , unbiased 3 s md simulations 
 ( total length ) and free energy profiles along the bilayer normal calculated 
 by biased md simulations ( 7 s in total ) . 
 the convergences 
 in time of free energy profiles calculated by both umbrella sampling 
 and z - constraint techniques are thoroughly analyzed . 
 two sets of starting 
 structures are also considered , one from unbiased md simulation and 
 the other from  pulling  coumarin along the bilayer 
 normal . 
 the structures obtained by pulling simulation contain water 
 defects on the lipid bilayer 
 surface , while those acquired from unbiased simulation have no membrane 
 defects . 
 the free energy profiles converge more rapidly when starting 
 frames from unbiased simulations are used . 
 in addition , z - constraint 
 simulation leads to more rapid convergence than umbrella sampling , 
 due to quicker relaxation of membrane defects . 
 furthermore , we show 
 that the choice of resp , prodrg , or mulliken charges considerably 
 affects the resulting free energy profile of our model drug along 
 the bilayer normal . 
 we recommend using z - constraint biased md simulations 
 based on starting geometries acquired from unbiased md simulations 
 for efficient calculation of convergent free energy profiles of druglike 
 molecules along bilayer normals . 
 the calculation of free energy profile 
 should start with an unbiased simulation , though the polar molecules 
 might need a slow pulling afterward . 
 results obtained with the recommended 
 simulation protocol agree well with available experimental data for 
 two coumarin derivatives . 
Introduction
Methods
Results
Discussion
Conclusion