the lipid bilayer plays critical roles in the biochemistry of cells , the most basic of which being its role in defining the shape of the cell and organelles .
the unique arrangement of hydrophilic and hydrophobic groups at the membrane exterior and interior regions , respectively , allows for a selective inhibition of foreign agents and transport of essential molecules . as a consequence ,
the functional versatility of the cell membrane is also in tune with its highly inhomogeneous composition and structure , wherein receptors , channels , peptides , and other molecules create a vast chemical mixture and diverse geometry.(5 ) the thermodynamics of adhesion and membrane insertion of biomolecules is therefore the subject of much research . in the context of protein folding , for instance ,
understanding how membrane proteins insert and simultaneously fold to assume their functional three - dimensional structure is of great interest.(6 ) this poorly understood process can involve the cooperative action of several amino acid residues that work in concert toward both membrane penetration and organization to a particular shape .
researchers are also studying how smaller biomolecules insert into the membrane and aggregate , and to what extent these processes affect the physical properties of the membrane .
several experimental methods have been employed to describe the equilibrium properties of peptide insertion , that is , to describe the initial and final states and the associated thermodynamic properties .
the most important thermodynamic parameter in this context is the free energy change , g , defined here to be the energetic change upon the transfer of the peptide from its initial aqueous environment to its final state in membrane ; g < 0 indicates a favorable insertion process .
studied the free energy of amino acid partitioning first in an octanol / water(12 ) membrane mimic and later in a synthetic membrane.(13 ) in the latter study , they derived a thermodynamic scale for the partitioning of small model peptides between membrane and water . by altering a single residue in the test peptide , the contribution of each of the 20 amino acids to the transfer free energy
for instance , in the case of the hydrophobic amino acid leucine , the free energy of transfer is 0.56  0.04 kcal mol .
also contemplated how peptide aggregation or formation of higher order molecular assemblies might occur in membranes .
concurrently , computational methods have played significant roles in elucidating the atomic - level details that are not always accessible by experiment .
the advance of computational power and the development of faster algorithms enables modeling of increasingly larger membrane systems .
for example , aliste and tieleman used md simulations to study the partitioning of the wimley and white peptides in a model membrane , and the results provided the details of the peptidemembrane atomic interactions responsible for insertion.(21 ) the partitioning of the various amino acids into a model membrane host has also been studied.(22 ) a unique benefit of computation is its ability to shed light on the free energy of peptidemembrane binding , as well as to characterize the energetic landscape or the pathways of insertion and aggregation .
one such study examined the free energy profile ( or the potential of mean force , pmf ) of the chromophore indole ring as it traveled through the membrane.(23 ) in this work , the authors explored how different parametrizations of the indole moiety can affect the free energy calculations ; they also pinpointed some key properties of the reaction coordinate , including the locations of energetic barriers . in this study , we employ a molecular dynamics ( md ) technique known as  umbrella sampling  to simulate the insertion of a model hydrophobic peptide into a membrane .
the model of choice is the small hexapeptide consisting of one tryptophan ( trp ) and five leucine residues ( wl5 ) , the same model system used in the experiments of wimley and white.(13 ) we computed the pmf as the peptide moves from the solvent and inserts into and completely traverses across the membrane . to our knowledge , this is the first attempt in computing the free energy profile of a full length peptide as it crosses a membrane .
we analyzed the dynamical and structural properties of the peptide during this process . by coupling the computed free energy profile with the experimental results of wimley and white ,
we propose a thermodynamic model for the insertion and aggregation of hydrophobic peptides in a membrane host .
each leaflet of the membrane can be divided into three regions based on the polarity of the constituent atoms : the solventheadgroup interface ( s / hg ) , where the first two solvent shells (  6  ) merge with the choline and phosphate head groups of the phospholipids ; the glycerol or headgroupcore interface ( hg / c ) , where the lipid head groups mix with the hydrophobic fatty acid chains ; and the core , the region occupied by the aliphatic lipid tails . by using these three regions as landmarks , the free energy profile of the peptide wl5 as it traverses across the membrane
see supporting information for a figure ( s1 ) demonstrating the extent of convergence of this profile ( in the allotted sampling time ) .
the profile exhibits minima at both the s / hg and the hg / c interfaces of each leaflet , the latter being the global minimum .
note that two peptide orientations ( one per leaflet ) were used to determine whether the initial orientation affects the energetics of insertion .
the dashed line in figure 1 clearly shows that the initial orientation does indeed affect the profile .
such an asymmetric profile implies that even for a rudimentary peptide such as wl5 , there may be more than one possible path of membrane insertion .
free energy profile ( symmetric - heavy blue line , asymmetric - dashed line ) of the peptide as a function of position along the z dimension ( negative / positive values correspond to the lower / upper leaflets , respectively ) .
pink marks the regions of the solvent / lipid headgroup interface ( s / hg in the lower leaflet , hg / s in the upper ) .
beige indicates the headgroup / core interface ( hg / c in the lower leaflet , c / hg in the upper ) .
the conformation of the peptide , as judged by the distribution of / angles , assumes that of a  strand , especially when the peptide is positioned somewhere in the membrane ; see figure s2 of supporting information .
most noticeably , the distribution centers on the  strand / angles when the peptide is embedded in the hg / c region of the membrane ( lower left panel of figure s2 ) .
such a conformational change may have some effect on the free energy profile of the peptide as it travels through the membrane .
yet , it is difficult to quantify this effect because we have not imposed any constraints on the / angles of the peptide . in the upper leaflet ,
the nitrogen heteroatom points toward the interfacial regions while the peptide moves through the membrane ; while in the lower leaflet , the ring nitrogen consistently points toward the core of the membrane ( see figure 2 and figure 3 ) .
this preferential orientation results in distinct and deeper energy minima at the interfacial regions of the upper leaflet while the corresponding minima in the lower leaflet are shallow and less distinct ( see figure 1 and figure 2b ) .
the location of the minima and the orientation of the indole ring in the upper leaflet indicate that the ring is in a prime position to interact with the glycerol carbonyls and other hydrogen bonding elements of the lipid head groups .
the importance and extent of hydrogen bonding between the indole ring and the membrane interfacial region has been addressed in a previous study.(26 ) here , the nitrogen heteroatom forms a significant number of hydrogen bonds with the membrane , and it does so to a greater extent where the minima occur along the reaction coordinate ( see figure 4 ) . in the upper leaflet , where the minima are deeper and more distinct ,
a larger number of hydrogen bonds are formed as compared with the minima in the lower leaflet .
( a ) density distributions of the tryptophan indole nitrogen in the upper and lower leaflets of four simulation windows .
dashed lines indicate where the peptide was constrained ( via the indole ring ) in that particular simulation .
( b ) the asymmetric free energy profile in the membrane ( dashed line ) and the same profile shifted to reflect the nitrogen heteroatom position ( dark heavy line ) . at each point , the nitrogen is distributed to the right ( more positive z ) of the indole center ; but the magnitude of that shift is not equal in all parts of the profile .
snapshots of wl5 in the upper ( above ) and lower ( below ) leaflets of the membrane . the membrane is shown as a transparent surface with the lipid head groups colored pink and the tails brown .
the peptide is shown as van der waals spheres with the indole ring highlighted in yellow .
the insets show a close - up of the indole ring orientation , with carbon , hydrogen , and nitrogen atoms colored in cyan , white , and blue , respectively .
average number of hydrogen bonds formed between the indole nitrogen and the hydrogen bonding elements of the lipids , as a function of the peptide s position along the reaction coordinate .
the cutoff criterion was set to elements within 3  and hydrogen bonding angles of 40. analysis with respect to orientation of the entire peptide reveals a tilt in the peptide axis ( see figure 5 ) , an observation common in helical transmembranepeptide systems.(27 ) as it adheres and inserts into the membrane , the peptide axis stabilizes at an angle near 70 , which is off parallel to the membrane plane ( an angle of 90 ) .
such a tilt allows the leucine residues to interact with the hydrophobic part of the membrane while maintaining the indole ring position in the interfacial region .
this appears to be the general strategy of transmembrane peptides , in that they are designed to exploit both the hydrophilic and the hydrophobic parts of the membrane during the insertion process .
the orientations of the whole peptide and the tryptophan side chain argue that the pmf in the upper leaflet represents a more plausible path of peptide insertion .
but in reality , the peptide could approach the membrane in a variety of orientations .
also , it may be the case that the peptide binds the membrane surface in a random initial orientation and then rearranges prior to further insertion and stabilization .
peptide tilt as a function of position along the z axis , with the membrane centered on 0 .
inset shows how the angle is defined by the vector fitted to contain the alpha carbons of the peptide ( light arrow ) and the vector representing the membrane normal ( dark arrow ) .
however , we do not have any a priori knowledge of the peptide orientation or information on its conformational dynamics at the solventmembrane interface .
furthermore , the constraint applied at the trp residue during the umbrella sampling does not allow for a spontaneous reorganization .
it is therefore impossible to determine the relative weights of the two profiles , and here we assume that the lower- and upper - leaflet orientations are equally probable . in other words , because the molecular dynamics simulations in each leaflet are separate and independent , one can consider the data from the upper and lower leaflets as twice the sampling of the same system .
thus , we symmetrize the free energy profile ( heavy blue line in figure 1 ) and allow this result to reflect the sampling of the entire system ( see materials and methods ) .
this symmetric profile is used for all further calculations and discussion , but we note that this symmetry is due purely to the homogeneity of the model membrane used in the simulations . in biological heterogeneous membranes ,
the peptide begins its approach toward the membrane from a region of bulk solvent , where the pmf plateaus ( 35  ) .
this plateau indicates that the peptide does not interact with the membrane and instead experiences a purely aqueous environment .
the average pmf in this region ( aqua - colored in figure 1 ) is 10.1 kcal mol .
this value will be used as the reference point from which the change in free energy ( g ) will be calculated . in other words ,
as the peptide travels closer to the membrane , at about 30  from the center , the free energy rapidly declines as the peptide encounters the solventheadgroup ( s / hg ) interface .
such a drop in the free energy is often encountered when a hydrophobic solute encounters a like environment.(28 ) inside this s / hg region , the free energy profile continues to decline until it reaches a shallow minimum near 17  , and the barrier to escape this minimum is approximately 0.5 kcal mol .
[ this barrier seems insignificant given the sampling error of this method ( see materials and methods ) , and it should be greater . nonetheless , we know from previous experiment and simulation that a minimum in the free energy profile occurs here . ] the computed free energy change as the peptide inserts into the s / hg interface is gs / hg = 8.3 kcal mol .
thus , the peptide has a tendency to localize in the s / hg region , and this can likely be attributed to the partial hydrophilic nature of the trp residue .
the peptide descends further into the headgroupcore ( hg / c ) interface and encounters a global minimum at about 7  and gains a further 2 kcal mol in free energy .
the total free energy change from the bulk solvent to the hg / c region ( ghg / c ) is 10.2 kcal mol .
a local maximum occurs in the core of the membrane ( 0  ) , and the energy barrier as measured from the hg / c region is approximately 1.0 kcal mol .
the root - mean - square fluctuation ( rmsf ) demonstrates how much the peptide deviates from its average conformation and is plotted in figure 6a as a function of the peptide location along the membrane normal .
one can conclude that the peptide is more flexible and samples a wider region of conformational space ( has a larger rmsf ) as it travels deeper into the membrane core .
this flexibility is also demonstrated by the root - mean - square deviation ( rmsd ) of the peptide structure , as seen in figure s3 of supporting information .
an alternative view of the rmsf in figure 6b further shows that the side chains of the peptide are predominantly responsible for this flexibility . in particular , the trp side chain has the largest rmsf when deep in the membrane core .
to the extent that peptide flexibility can be correlated with entropy,(29 ) one can infer that the peptide exhibits increased entropy as it inserts into the membrane .
this results in a more negative ( or favorable ) contribution to the free energy change and explains why the energy barrier in the core is not greater in the calculated profile .
therefore , a small barrier seems to suggest that the peptide can occasionally translate across the membrane center and shift positions between the two leaflets . yet , as wimley and white have demonstrated and as we shall explore further below , individual units of this model peptide aggregate in the membrane and form a supramolecular structure , a process that can counteract the tendency of a lone peptide to transfer from leaflet to leaflet .
root mean square fluctuation ( rmsf ) of the peptide by atom number and as a function of position along the z axis of the simulation box .
atom numbers 13 correspond to the acetyl group on the n - terminus ; 424 correspond to the tryptophan residue ; 2569 correspond to the leucine residues ; and 7072 correspond to the amide group on the c - terminus .
( a ) surface representation , where the red halo shows the maximum rmsf achieved by the indole ring near the core of the membrane .
the free energy profile presented here is that of a single peptide traveling through the membrane and is thus not reflective of the entire biological process .
if the final state of wl5 in the membrane is not monomeric but rather an aggregate form , then it behooves us to study this aggregation process and the pathways that lead to it .
although we have not simulated the actual aggregation , the presented work , coupled with previous simulations and experimental evidence , can shed some light on the mechanism and potential pathways of aggregation .
figure 7 shows a proposed thermodynamic cycle for the insertion and assembly of wl5 in the membrane .
in this model , the initial state of the peptide in an aqueous environment is monomeric , and the final state in membrane is an assembly of several monomers ( top left and bottom right of figure 7 , respectively ) .
the difference in the free energy between these two states as obtained by wimley and white is gexp = 5.3 kcal mol ( exp = experimental).(13 ) [ this is assuming that the transfer from monomeric aqueous to membrane aggregate form is an equilibrium process , and that the experimental result is a reflection of this . ] as suggested by grossfield et al .
, the computed free energy can be coupled to this experimental value via a correction term that accounts for the difference in the peptide and lipid concentrations between simulation and experiment.(30 ) proposed thermodynamic cycle of peptide insertion and supramolecular assembly inside a model membrane .
the yellow arrows depict a pathway in which single peptides insert first into the deep headgroupcore ( hg / c ) interface and then aggregate to form the final structure .
the red arrows mark a less likely ( yet plausible ) path , in which the peptide aggregates first in solution , and then the aggregate as a whole inserts into the membrane .
the green arrows qualitatively describe an intermediary pathway , in which individual peptides insert first into the aqueous or solventheadgroup ( s / hg ) interface , aggregate to some extent , and then insert deeper into the membrane .
each of these steps can be described by a free energy change ( g , aggr = aggregation ) .
following the notation of grossfield , gr and gmix are the correction terms for the peptide and lipid ( respectively ) concentration differences between simulation ( sim ) and experiment ( exp ) .
 is the mole fraction of the species in question ( either peptide or lipids ) .
this entire correction term , gcorr , is calculated to be + 3.0 kcal mol ( see materials and methods for further information ) and must be added to the computational results . in terms of insertion and aggregation , figure 7 shows two paths on opposite ends of the spectrum . in the yellow path , each peptide inserts individually into the deeper membrane interface .
the aggregation step then follows : 

 equation 5 is the corrected version of eq 4 , where the correction term has been added to computational result ghg / c = 10.2 kcal mol .
if we can set gyellowpath to be the experimental value determined by wimley and white ( 5.3 kcal mol ) , the aggregation term is then calculated to be + 1.9 kcal mol .
the red path in figure 7 illustrates the opposite approach , in which individual peptides in aqueous solution first assemble , and then the whole aggregate inserts into the membrane .
the first leg of this pathway is plausible and has a negative free energy change , because hydrophobic peptides in an aqueous environment are energetically more stable in aggregate instead of monomeric form.(31 ) yet , the stability and energetic gain of the aggregate in solution could deter the adhesion and insertion of the aggregate into the membrane , thus making the second leg of the red path less likely .
furthermore , if the orientation of each monomeric peptide can affect the thermodynamics of the insertion process ( as discussed above ) , it may become energetically costly for the peptide units of the complex to reorganize and assume the correct orientation .
others have also shown that transmembrane proteins do not fold entirely or form complex structures before inserting into a membrane.(32 ) instead , it is more likely that a transmembrane peptide adsorbs on the membrane surface , aggregates to some extent , and then inserts into the membrane . as the insertion event proceeds , it is likely that the membrane continues to shape the peptide aggregate structure until it achieves a final stable form in the membrane host . note that membranes exhibit the mechanical properties to do just this.(33 ) for our model system , this is depicted in the green path , in which the monomeric peptide first adheres to the s / hg interface , aggregates , and then inserts : 

 equation 7 represents the corrected version , with the correction term added again to the computational result , gs / hg = 8.3 kcal mol .
the last two terms on the left - hand side of equation 7 are the unknowns , and we can not calculate the exact contribution from gaggr(s / hg ) or gaggr(insert ) . again , if we equate ggreenpath to the experimental value of 5.3 kcal mol , the sum of these unknown terms can be calculated , gaggr(s / hg ) + gaggr(insert ) = 0 kcal mol .
thus , in the yellow path , the aggregation free energy change is positive and indicative of an unfavorable process . in other words ,
if the monomeric form of the peptide were to insert deep into the hg / c region of the membrane , aggregation would not be likely . in the green path
however , the collective aggregation term is 0 , indicating an equilibrium ( and more likely ) process .
we have shown the following relation : 

 where the braces group the unknown terms of each path , as described above .
the results suggest that the green path in figure 7 , where the peptide first adheres to the solventheadgroup interface , aggregates , and then inserts , more closely represents the in vitro mechanism of insertion and organization of wl5 .
furthermore , following the previous discussion of the red path in which we stated that a formed aggregate would have difficulty inserting into the membrane , we can speculate that gaggr(insert ) > 0 and thus gaggr(s / hg ) < 0 .
in other words , aggregation occurs immediately upon peptide binding to the membrane , in the solventheadgroup region .
the structure of the solvated peptidemembrane system was obtained from a previous study , in which we explored the properties of the system in a 75 ns md simulation.(25 ) snapshots at the end of this simulation , at which point , the peptides were embedded within the headgroup region of the membrane , were used as starting configurations to initiate 50 different all - atom simulations each representing a window in an umbrella sampling scheme .
two different snapshots of the peptide were employed ( one for each leaflet of the membrane ) ; in other words , the 25 simulation windows within a leaflet were initiated with the same initial peptide conformation .
these two snapshots are taken to be the equilibrated conformation of the peptide when embedded in the solventheadgroup interfacial region of the membrane . in total
then , the 50 simulation windows consisted of a wl5 peptide ( acetylated and amidated at its n- and c - termini , respectively ) , a preconstructed bilayer of 128 dimyristoylphosphatidylcholine ( dmpc ) lipids,(20 ) and 6631 water molecules , resulting in a 25  000 atom system assembled in a simulation box of 60  60  95  . the membrane was positioned at the center of the box , leaving approximately 23  for the solvent on either side .
each simulation corresponds to a window of width 1.9  such that the position of wl5 was different in each simulation . in the first window ,
wl5 was constrained ( along the z axis ) in the bulk solvent at the upper half of the box and well away from the membrane surface . in the next ,
wl5 was moved by 1.9  closer to the membrane and constrained at that position . in subsequent simulations ,
wl5 was moved closer to the center of the simulation box each time stepped along the z axis by the increment of 1.9  .
the constrained molecule thus spans the entire reaction coordinate ( in the z direction ) of approximately 95  . for organizational purposes ,
the simulation box is divided into two halves along the z axis , representing the  upper  ( positive z ) and  lower  ( negative z ) leaflets with respect to the membrane center at z = 0 . among the 25 simulation windows in each leaflet , 12 had the peptide constrained in solvent and 13 were carried out with the peptide constrained at various z locations in the membrane . the initial configuration of the peptide differs between the two leaflets but is the same among the 25 simulation windows of each leaflet .
the reaction coordinate was defined by the separation along the z axis between the centers - of - mass of the indole chromophore , the tryptophan side chain , and the membrane .
wl5 was constrained by a harmonic potential of force constant 500 kj mol nm ( 1.20 kcal mol  ) .
first , unfavorable contacts were relieved by two cycles of 5000 steps steepest descent followed by 5000 steps of conjugate - gradient energy minimizations , each with the peptide held fixed and set free .
second , except when the peptide is in bulk solvent , an annealing step was conducted to equilibrate the lipid tails around the peptide . with
the peptide and lipid head groups restrained , the system was gradually heated ( over 1 ns ) from 310 k ( the temperature of the initial snapshot ) to 410 k , then gradually cooled back down to 310 k. third , adding pressure coupling , a restrained md was performed on all 50 windows for 1 ns , with the restraint applied to nonsolvent molecules .
fourth , with only the indole ring of tryptophan constrained to its selected position , production runs commenced .
those with the peptide located outside of the membrane were sampled for 5 ns , while those with the peptide somewhere in the membrane were sampled for 25 ns . in total , including all pre- and post production runs , approximately 850 ns of all - atom md simulations have been carried out .
a time step of 2 fs was used while coordinates / velocities were recorded every 500 steps ( 1 ps )
full electrostatics were calculated using the particle mesh ewald ( pme ) method , with coulombic and van der waals cut - offs of 0.9 and 1.4 nm , respectively .
a sixth - order spline interpolation was used for pme along with a tolerance of 1  10 .
nearest neighbor lists were updated every 10 steps using the grid method and periodic boundary conditions ( in xyz ) were employed with a cutoff of 0.9 nm .
berendsen temperature coupling was used , with a reference temperature of 310 k and a coupling of 0.1 ps .
a semi - isotropic pressure coupling scheme , in which the xy dimensions are coupled together while the z direction is allowed to fluctuate independently , was used with a reference pressure of 1.0 bar , a pressure coupling of 0.5 ps , and a compressibility of 4.5  10 bar .
the simulations were performed with the gromacs md package.(34 ) the ffgmx force field was used for the peptide and the spc water model was used as the solvent .
the lipid parameters of berger et al . were employed for the dmpc lipids.(35 ) the analysis of the simulation trajectories was completed using the various gromacs tools . for further details regarding these tools and specific procedures in running gromacs md simulations
; see the gromacs manual at http://www.gromacs.org/. visualization and rendering were done by the visual molecular dynamics ( vmd ) application.(36 ) given the construction of the simulation windows , one might expect the resulting profile in figure 1 to be perfectly symmetric about the membrane center .
in other words , the profile in one leaflet should be a mirror image of the other provided that the peptide exhibits the same orientation while penetrating through both leaflets of a homogeneous membrane . as mentioned in the text , the free energy profile of figure 1 is not exactly symmetric , and this was attributed mainly to the starting conformation of the peptide in the simulation windows of each leaflet . yet , even if the peptide conformation was the same in each leaflet , it is doubtful that one could achieve great symmetry , simply due to convergence limitations inherent in md simulations .
in fact , the asymmetric profile may be used as a convergence criterion with respect to certain free energy calculations .
for instance , the values of the computed ghg / c in the upper and lower leaflets are within 1 kcal mol of each other while the values of gs / hg differ more significantly .
this could signify that the simulation of deep insertion ( from bulk solvent into the deep interface ) produces a more converged result than insertion into the solventheadgroup ( s / hg ) interface .
one could reason that this is the case because the s / hg interface is more chemically diverse and dynamic , thus requiring more sampling time for energetic calculations . in figure 1 ,
about the center of the profile ( at 0  ) , the average was taken between the two points on either side of the profile .
for instance , if the free energy was 1.8 kcal mol at + 10  ( upper leaflet ) and 1.5 kcal mol at 10  ( lower leaflet ) , the free energy in the symmetric plot at 10  is taken to be 1.65 kcal mol . the asymmetric free energy profile in figure 1
the differences at mirror sites along the reaction coordinate ( for instance , at 7  ) reveal a sampling error on the order of 0.51.0 kcal mol . as reported by wimley and white ,
the experimental concentration of the wl5 peptide is ( at most ) 100 m.(13 ) using this concentration , and that of water to be 55.5 m ( or 55.5  10 m ) , we can calculate the mole fraction of the peptide in the experiment to be exp , peptide = 1.80  10 . in each simulation
, we have one peptide , in either the upper or lower leaflet of the system . because of periodic boundary conditions , the peptide is solvated in the total solvent of 6631 water molecules .
thus , the mole fraction of the peptide in the simulation is sim , peptide = 1.51  10 . substituting these mole fractions along with the simulation temperature ( 310 k ) into the expression for gr
, we arrive at + 2.73 kcal mol , which corrects for the differences in peptide concentration between simulation and experiment . to calculate the mole fractions of the lipids , a volumetric ratio is computed in a fashion similar to grossfield et al .
taking the average thickness of a dmpc membrane to be 36  and the cross - sectional area to be 60  60  ( per simulation setup ) , we compute the volume occupied by the membrane to be 1.30  10(5 )  .
when the ratio of the membrane volume to the box volume is computed , we obtain sim , lipids = 0.38 ; we also use exp , lipids = 0.24 , for an experimental vesicle concentration of 4 mm.(29 ) these mole fractions result in gmix = + 0.28 kcal mol , and this corrects for the differences in lipid concentration between simulation and experiment .
the combined correction term is then the sum , gr + gmix = gcorr  + 3.0 kcal mol .
 a variety of biomolecules mediate physiological processes by inserting and reorganizing in cell membranes , and the thermodynamic forces responsible for their partitioning are of great interest . 
 recent experiments provided valuable data on the free energy changes associated with the transfer of individual amino acids from water to membrane . 
 however , a complete picture of the pathways and the associated changes in energy of peptide insertion into a membrane remains elusive . to this end , computational techniques supplement the experimental data with atomic - level details and shed light on the energetics of insertion . here , we employed the technique of umbrella sampling in a total 850 ns of all - atom molecular dynamics simulations to study the free energy profile and the pathway of insertion of a model hexapeptide consisting of a tryptophan and five leucines ( wl5 ) . 
 the computed free energy profile of the peptide as it travels from bulk solvent through the membrane core exhibits two minima : a local minimum at the watermembrane interface or the headgroup region and a global minimum at the hydrophobichydrophilic interface close to the lipid glycerol region . a rather small barrier of roughly 1 kcal mol1 exists at the membrane core , which is explained by the enhanced flexibility of the peptide when deeply inserted . 
 combining our results with those in the literature , we present a thermodynamic model for peptide insertion and aggregation which involves peptide aggregation upon contact with the membrane at the solventlipid headgroup interface . 
Introduction
Results and Discussion
Materials and Methods